Commit 6c1cc9a2 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

reformat code

parent 87e61f43
......@@ -4,7 +4,7 @@ D3 := output/processed
D4 := figs/test.format.tree
D5 := figs/test.traits
D6 := figs/test.CWM
sites:= Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS US
sites:= US Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS
D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) ))
D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) ))
D2tree := $(addsuffix /tree.csv,$(addprefix $(D2)/, $(sites) ))
......
#######################################
#######################################
###### EXPLORE DATA SET BEFORE ANALYSIS
rm(list = ls())
source("R/utils/mem.R")
source("R/analysis/lmer.run.nolog.all.R")
filedir <- "output/processed"
traits <- c('Wood.density', "SLA", 'Leaf.N', 'Max.height')
for (t in traits){
print(t)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE)
print(1)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, cat.TF=TRUE)
print(2)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, fname = 'data.all.csv')
print(3)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, cat.TF=TRUE, fname = 'data.all.csv')
print(4)
}
#!/usr/bin/env Rscript
source("R/analysis/glmer.output-fun.R")
files <- list.files("output/glmer", recursive=TRUE, full.names =TRUE, pattern = "glmer.results.rds")
out <- lapply(files, summarise.glmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/glmer.list.out.rds')
names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$glmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
saveRDS(DF.results,file='output/glmer.DF.rds')
......@@ -5,25 +5,20 @@ library(lme4)
run.models.for.set.all.traits <- function(model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),
type.filling, ...){
for(trait in traits)
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, type.filling=type.filling, ...)
}
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait,
type.filling, ...){
type.filling, fname = 'data.all.no.log.csv',
cat.TF = FALSE, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait,
type.filling=type.filling, ...))
type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait,
type.filling, ...){
type.filling, fname, cat.TF, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling,...))
try(fun.model(model.file, trait, type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
......@@ -32,20 +27,34 @@ run.model.for.set.one.trait <- function(model.file,fun.model, trait,
#=====================================================================
model.files.glmer.Tf.1 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.R")
c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.R")
c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
model.files.glmer.Tf.3 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.MAT.MAP.R")
model.files.glmer.Tf.4 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.MAT.MAP.R")
fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
glmer.output <- glmer(formula=formula,data=df.lmer,family=binomial)
glmer.output <- glmer(formula=formula,data=df.lmer,
family=binomial(link = 'cloglog'))
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
......@@ -54,12 +63,14 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA,
type.filling) {
type.filling, fname , cat.TF) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
path.out <- output.dir.glmer(model$name, trait,
if(fname == 'data.all.no.log.csv') dir <- 'all.no.log'
if(fname == 'data.all.csv') dir <- 'all.log'
path.out <- output.dir.glmer(dir, model$name, trait,
type.filling=type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
......@@ -67,8 +78,10 @@ run.glmer <- function (model.file, trait,
'all',"trait",
trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
min.obs, sample.size,
type.filling=type.filling,
fname = fname,
cat.TF = cat.TF) # return a DF
cat("Ok data with Nobs",nrow(df.glmer),"\n")
#= Run model
fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
......@@ -76,8 +89,25 @@ run.glmer <- function (model.file, trait,
#========================================================================
output.dir.glmer <- function (model, trait, type.filling) {
file.path("output/glmer", 'all',trait,type.filling,model)
output.dir.glmer <- function (dir, model, trait, type.filling) {
file.path("output/glmer", dir,trait,type.filling,model)
}
fun.load.data.all <- function(base.dir,fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
intern = TRUE))
data.tree.tot <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE,
nrows = nrows,
colClasses = classes)
return(data.tree.tot)
}
......@@ -86,37 +116,58 @@ output.dir.glmer <- function (model, trait, type.filling) {
#============================================================
load.and.prepare.data.for.glmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
fname = 'data.all.no.log.csv',
base.dir = "output/processed/",
cat.TF = FALSE){
### load data
data.all.sample <- read.csv(file=file.path(base.dir, "data.all.csv"),
stringsAsFactors=FALSE,nrows=1000)
classes <- sapply(data.all.sample,class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',
intern=TRUE))
data.tree.tot <- read.csv(file=file.path(base.dir, "data.all.csv"),
stringsAsFactors=FALSE,
nrows=nrows,
colClasses=classes)
fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- readRDS(file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
return(df)
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,
load.and.save.data.for.glmer <- function(trait,
min.obs = 10, sample.size = NA, type.filling = 'species',
fname = 'data.all.no.log.csv',
base.dir = "output/processed/", cat.TF = FALSE){
### load data
data.tree.tot <- fun.load.data.all(base.dir,fname)
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling, cat.TF = cat.TF)
saveRDS(df,file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs,
min.perc.neigh=0.80){
min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) )
(!is.na(data.tree[["D"]])) &
(!is.na(data.tree[['year']])))
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh &
!is.na(data.tree[[perc.neigh]]))
# select obs abov min perc.neigh species
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
min.perc.neigh.sp &
!is.na(data.tree[[perc.neigh.sp]]))
# select obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs & !is.na(data.tree[[perc.neigh.gs]]))
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(
......@@ -128,10 +179,18 @@ fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}
fun.standardized.var <- function(x){
return(x/sd(x))
}
### get variables for lmer
fun.get.the.variables.for.glmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,
abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
logyear <- log(data.tree[["year"]])
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
......@@ -145,17 +204,82 @@ sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(dead=dead,
logD=logD,
logyear=logyear,
MAT = MAT,
MAP = MAP,
species.id=species.id,
tree.id=tree.id,
plot.id=plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.diff= fun.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs= fun.standardized.var(sumTnTfBn.abs),
Tf= fun.standardized.var(data.tree[[tf]]),
sumTnBn=fun.standardized.var(sumTnBn),
sumTfBn=fun.standardized.var(sumTfBn),
sumBn=fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
}
### get variables for lmer
fun.get.the.variables.for.glmer.no.tree.id.cat <- function(data.tree,BATOT,CWM.tn,
abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
logyear <- log(data.tree[["year"]])
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])
#get the three cwm per cat
vec.CWM.tn <- paste(CWM.tn, c('A_EV', 'A_D', 'C'), sep = '.')
vec.abs.CWM.tntf <- paste(abs.CWM.tntf, c('A_EV', 'A_D', 'C'), sep = '.')
sumTnTfBn.abs.A_EV<- data.tree[[vec.abs.CWM.tntf[1]]]
sumTnTfBn.abs.A_D <- data.tree[[vec.abs.CWM.tntf[2]]]
sumTnTfBn.abs.C <- data.tree[[vec.abs.CWM.tntf[3]]]
sumTnBn.A_EV<- data.tree[[vec.CWM.tn[1]]]
sumTnBn.A_D <- data.tree[[vec.CWM.tn[2]]]
sumTnBn.C <- data.tree[[vec.CWM.tn[3]]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *sumTfBn
sumBn <- data.tree[[BATOT]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
return(data.frame(dead=dead,
logD=logD,
logyear=logyear,
MAT = MAT,
MAP = MAP,
species.id=species.id,
tree.id=tree.id,
plot.id=plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.diff= fun.center.and.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs= fun.center.and.standardized.var(sumTnTfBn.abs),
Tf= fun.center.and.standardized.var(data.tree[[tf]]),
sumTnBn=fun.center.and.standardized.var(sumTnBn),
sumTfBn=fun.center.and.standardized.var(sumTfBn),
sumBn=fun.center.and.standardized.var(sumBn)))
sumTnTfBn.abs.A_EV= fun.standardized.var(sumTnTfBn.abs.A_EV),
sumTnTfBn.abs.A_D= fun.standardized.var(sumTnTfBn.abs.A_D),
sumTnTfBn.abs.C= fun.standardized.var(sumTnTfBn.abs.C),
Tf.A_EV= fun.standardized.var(Tf.A_EV),
Tf.A_D= fun.standardized.var(Tf.A_D),
Tf.C= fun.standardized.var(Tf.C),
sumTnBn.A_EV=fun.standardized.var(sumTnBn.A_EV),
sumTnBn.A_D=fun.standardized.var(sumTnBn.A_D),
sumTnBn.C=fun.standardized.var(sumTnBn.C),
sumTfBn.A_EV=fun.standardized.var(sumTfBn.A_EV),
sumTfBn.A_D=fun.standardized.var(sumTfBn.A_D),
sumTfBn.C=fun.standardized.var(sumTfBn.C),
sumBn=fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
}
......@@ -164,7 +288,7 @@ return(data.frame(dead=dead,
# that will be used in a simple linear model
#============================================================
fun.data.for.glmer <- function(data.tree,trait,min.obs=10,
type.filling='species') {
type.filling='species', cat.TF = FALSE) {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA",
"Wood.density","Max.height"))
stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
......@@ -174,11 +298,31 @@ CWM.tn <- paste(trait,"CWM",'fill',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".")
data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,
perc.neigh,BATOT,min.obs)
perc.neigh.sp <- paste(trait, "perc", 'species', sep = ".")
perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
data.tree <- fun.select.data.for.analysis(data.tree,
abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs)
#= DATA LIST FOR JAGS
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
if(!cat.TF) {
print('no cat')
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,
BATOT,
CWM.tn,
abs.CWM.tntf,
tf)
}
if(cat.TF) {
print('cat')
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id.cat(data.tree,
BATOT,
CWM.tn,
abs.CWM.tntf,
tf)
}
return(glmer.data)
}
......
###########################
###########################
### FUNCTION TO RUN GLMER ESTIMATION
library(lme4)
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait,
type.filling, fname = 'data.all.no.log.csv',
cat.TF = FALSE, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait,
type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait,
type.filling, fname, cat.TF, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
#=====================================================================
# Run glmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.glmer.Tf.1 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.R")
model.files.glmer.Tf.3 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.MAT.MAP.R")
model.files.glmer.Tf.4 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.MAT.MAP.R")
fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
glmer.output <- glmer(formula=formula,data=df.lmer,
family=binomial(link = 'cloglog'))
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.global.rds"))
}
run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA,
type.filling, fname , cat.TF) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
if(fname == 'data.all.no.log.csv') dir <- 'all.no.log'
if(fname == 'data.all.csv') dir <- 'all.log'
path.out <- output.dir.glmer(dir, model$name, trait,
type.filling=type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run glmer for model",model.file," for set",
'all',"trait",
trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait,
min.obs, sample.size,
type.filling=type.filling,
fname = fname,
cat.TF = cat.TF) # return a DF
cat("Ok data with Nobs",nrow(df.glmer),"\n")
#= Run model
fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
}
#========================================================================
output.dir.glmer <- function (dir, model, trait, type.filling) {
file.path("output/glmer", dir,trait,type.filling,model)
}
fun.load.data.all <- function(base.dir,fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
intern = TRUE))
data.tree.tot <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE,
nrows = nrows,
colClasses = classes)
return(data.tree.tot)
}
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait,
min.obs, sample.size, type.filling,
fname = 'data.all.no.log.csv',
base.dir = "output/processed/",
cat.TF = FALSE){
### load data
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- readRDS(file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
return(df)
}
load.and.save.data.for.glmer <- function(trait,
min.obs = 10, sample.size = NA, type.filling = 'species',
fname = 'data.all.no.log.csv',
base.dir = "output/processed/", cat.TF = FALSE){
### load data
data.tree.tot <- fun.load.data.all(base.dir,fname)
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling, cat.TF = cat.TF)
saveRDS(df,file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs,
min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) &
(!is.na(data.tree[['year']])))
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
# select obs abov min perc.neigh species
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
min.perc.neigh.sp &
!is.na(data.tree[[perc.neigh.sp]]))
# select obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs & !is.na(data.tree[[perc.neigh.gs]]))
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(
data.tree[["sp"]]))>min.obs])