glmer.run.R 8.19 KiB
###########################
###########################
### FUNCTION TO RUN GLMER ESTIMATION
library(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) }