An error occurred while loading the file. Please try again.
-
Georges Kunstler authored01d2e58c
######################################################### FUNCTION TO RUN GLMER ESTIMATIONlibrary(lme4)get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){ sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]}run.models.for.set.all.traits <- function(set,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, set, type.filling=type.filling, ...)}run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){ for (m in model.files) try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling, ...))}run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){ fun.model <- match.fun(fun.model) for (e in ecoregions) try(fun.model(model.file, trait, set, e, type.filling=type.filling,...))}#=====================================================================# Run glmer() (in package lme4) for one ecoregion, one trait, one model#=====================================================================model.files.glmer.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.R", "R/analysis/model.glmer/model.glmer.LOGLIN.R.R", "R/analysis/model.glmer/model.glmer.LOGLIN.ER.R")model.files.glmer.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.R", "R/analysis/model.glmer/model.glmer.LOGLIN.AD.R", "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.R")model.files.glmer.Tf.1 <- 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.ER.Tf.R")model.files.glmer.Tf.2 <- 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.HD.Tf.R", "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")fun.test.if.multi.census <- function(data){return("tree.id" %in% names(data))}fun.call.glmer.and.save <- function(formula,df.lmer,path.out){ Start <- Sys.time() glmer.output <- glmer(formula=formula,data=df.lmer,family=binomial) end <- Sys.time() print(end - Start) print(summary(glmer.output)) saveRDS(glmer.output,file=file.path(path.out, "glmer.results.rds"))}run.glmer <- function (model.file, trait, set, ecoregion, min.obs=10, sample.size=NA, type.filling) { require(lme4) source(model.file, local = TRUE) model <- load.model() #= Path for output path.out <- output.dir.glmer(model$name, trait, set, ecoregion,type.filling=type.filling)7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run glmer for model",model.file," for set",
set,"ecoregion",ecoregion,"trait",
trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait, set, ecoregion,
min.obs, sample.size,
type.filling=type.filling) # 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 (model, trait, set, ecoregion,type.filling) {
file.path("output/glmer", set,ecoregion,trait,type.filling,model)
}
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait, set, ecoregion,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.csv"), stringsAsFactors = FALSE)
fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) )
## 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 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])
# 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]]))
return(data.tree)
}
fun.center.and.standardized.var <- function(x){
return((x-mean(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){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(dead=dead,
logD=logD,
species.id=species.id,
tree.id=tree.id,
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
sumBn=sumBn))
}
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"]]))
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(dead=dead,
logD=logD,
species.id=species.id,
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
sumBn=sumBn))
}
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
fun.data.for.glmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
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 ")
# get vars names
CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
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)
#= DATA LIST FOR JAGS
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
return(glmer.data)
}