########################### ########################### ### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA for all in one big model 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, ...){ for (m in model.files) (run.model.for.set.one.trait (m, fun.model, trait, type.filling = type.filling, ...)) } run.model.for.set.one.trait <- function(model.file, fun.model, trait, type.filling, ...){ fun.model <- match.fun(fun.model) try(fun.model(model.file, trait, type.filling = type.filling, ...)) } #===================================================================== # Run lmer() (in package lme4) for one ecoregion, one trait, one model #===================================================================== model.files.lmer.Tf.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.R.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.Tf.R") model.files.lmer.Tf.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.HD.Tf.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.R") model.files.lmer.Tf.3 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.R.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.Tf.MAT.MAP.R") model.files.lmer.Tf.4 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.HD.Tf.MAT.MAP.R", "R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.MAT.MAP.R") fun.call.lmer.and.save <- function(formula, df.lmer, path.out){ lmer.output <- lmer(formula = formula, data = df.lmer, REML = FALSE, control = lmerControl(optCtrl = list(maxfun = 40000) ) ) print(summary(lmer.output)) saveRDS(lmer.output, file = file.path(path.out, "results.nolog.all.rds")) } run.lmer <- function (model.file, trait, 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.lmer(model$name, trait, 'all', type.filling = type.filling) print(path.out) dir.create(path.out, recursive = TRUE, showWarnings = FALSE) cat("run lmer for model", model.file, " for trait", trait, "\n") df.lmer <- load.and.prepare.data.for.lmer(trait, min.obs, sample.size, type.filling = type.filling) # return a DF cat("Ok data with Nobs", nrow(df.lmer), "\n") #= Run model fun.call.lmer.and.save(formula = model$lmer.formula.tree.id, df.lmer, path.out) } #======================================================================== output.dir.lmer <- function (model, trait, set, type.filling) { file.path("output/lmer", set, trait, type.filling, model) } #============================================================ # Function to prepare data for lmer #============================================================ load.and.prepare.data.for.lmer <- function(trait, min.obs, sample.size, type.filling, base.dir = "output/processed/"){ ### load data ## library(data.table) ## data.tree.tot <- fread(file.path(base.dir, "data.all.csv"), ## stringsAsFactors = FALSE) 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.lmer(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.8, min.BA.G = -60, max.BA.G = 500){ ## remove tree with NA data.tree <- subset(data.tree, subset = (!is.na(data.tree[["BA.G"]])) & (!is.na(data.tree[["D"]])) ) ## remove tree with zero data.tree <- subset(data.tree, subset = data.tree[["BA.G"]]>min.BA.G & data.tree[["BA.G"]]9.9) ## select species with no missing traits select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) & !is.na(data.tree[[BATOT]])) data.tree <- data.tree[select.temp, ] # 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 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]) 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.lmer.tree.id <- function(data.tree, BATOT, CWM.tn, abs.CWM.tntf, tf, min.BA.G = 60){ logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G)) logD <- fun.center.and.standardized.var(log(data.tree[["D"]])) 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']]) #= 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(logG = logG, logD = logD, MAT = MAT, MAP = MAP, species.id = species.id, tree.id = tree.id, set.id = set.id, ecocode.id = ecocode.id, plot.id = plot.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))) } #============================================================ # Function to prepare data for lmer with new CWM data # that will be used in a simple linear model #============================================================ fun.data.for.lmer <- 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', 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) #= DATA LIST FOR LMER lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree, BATOT, CWM.tn, abs.CWM.tntf, tf) return(lmer.data) }