lmer.run.nolog.all.R

###########################
###########################
### 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"]]<max.BA.G &
                    data.tree[["D"]]>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)
}