lmer.run.nolog.all.R 9.29 KB
Newer Older
Georges Kunstler's avatar
Georges Kunstler committed
1
2
3
4
5
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA for all in one big model
library(lme4)

6
7
8
run.models.for.set.all.traits  <- function(model.file, fun.model,  traits =
         c("SLA", "Wood.density", "Max.height", "Leaf.N", "Seed.mass"),
                                           type.filling, ...){
Georges Kunstler's avatar
Georges Kunstler committed
9
  for(trait in traits)
10
11
   run.multiple.model.for.set.one.trait(model.file, fun.model, trait,
                                        type.filling = type.filling, ...)
Georges Kunstler's avatar
Georges Kunstler committed
12
13
}

14
15
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
                                                 type.filling, ...){
Georges Kunstler's avatar
Georges Kunstler committed
16
  for (m in model.files)
17
18
    (run.model.for.set.one.trait (m, fun.model, trait,
                                  type.filling = type.filling, ...))
Georges Kunstler's avatar
Georges Kunstler committed
19
20
21
}


22
23
run.model.for.set.one.trait <- function(model.file, fun.model, trait,
                                        type.filling, ...){
Georges Kunstler's avatar
Georges Kunstler committed
24
    fun.model <- match.fun(fun.model)
25
    try(fun.model(model.file, trait,  type.filling = type.filling, ...))
Georges Kunstler's avatar
Georges Kunstler committed
26
27
28
29
30
31
32
33
34
35
}


#=====================================================================
# 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")
36
37
model.files.lmer.Tf.2 <-
    c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.Tf.R",
Georges Kunstler's avatar
Georges Kunstler committed
38
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.Tf.R",
39
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.R",
Georges Kunstler's avatar
Georges Kunstler committed
40
41
42
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.HD.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.R")

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
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) ) )
Georges Kunstler's avatar
Georges Kunstler committed
59
   print(summary(lmer.output))
60
  saveRDS(lmer.output, file = file.path(path.out, "results.nolog.all.rds"))
Georges Kunstler's avatar
Georges Kunstler committed
61
62
}

Georges Kunstler's avatar
Georges Kunstler committed
63
run.lmer <- function (model.file, trait,
64
                      min.obs = 10, sample.size = NA,
Georges Kunstler's avatar
Georges Kunstler committed
65
66
67
68
69
70
                      type.filling) {
    require(lme4)
    source(model.file, local = TRUE)
    model <- load.model()
    #= Path for output
    path.out <- output.dir.lmer(model$name, trait, 'all',
71
                                type.filling = type.filling)
Georges Kunstler's avatar
Georges Kunstler committed
72
    print(path.out)
73
74
75
    dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
    cat("run lmer for model", model.file, " for trait", 
         trait, "\n")
Georges Kunstler's avatar
Georges Kunstler committed
76
      df.lmer <- load.and.prepare.data.for.lmer(trait,
Georges Kunstler's avatar
Georges Kunstler committed
77
                                           min.obs, sample.size,
78
79
80
                                         type.filling = type.filling)
                                        # return a DF
     cat("Ok data with Nobs", nrow(df.lmer),
Georges Kunstler's avatar
Georges Kunstler committed
81
82
         "\n")
        #= Run model
83
84
        fun.call.lmer.and.save(formula = model$lmer.formula.tree.id,
                               df.lmer, path.out)
Georges Kunstler's avatar
Georges Kunstler committed
85
86
87
88
89
}


#========================================================================

90
91
output.dir.lmer <- function (model, trait, set, type.filling) {
  file.path("output/lmer", set, trait, type.filling, model)
Georges Kunstler's avatar
Georges Kunstler committed
92
93
94
95
96
97
}


#============================================================
# Function to prepare data for lmer
#============================================================
Georges Kunstler's avatar
Georges Kunstler committed
98
99
load.and.prepare.data.for.lmer <- function(trait,
                                  min.obs, sample.size, type.filling,
Georges Kunstler's avatar
Georges Kunstler committed
100
101
                                  base.dir = "output/processed/"){
    ### load data
Georges Kunstler's avatar
Georges Kunstler committed
102
##     library(data.table)
Georges Kunstler's avatar
Georges Kunstler committed
103

104
## data.tree.tot <- fread(file.path(base.dir, "data.all.csv"),
Georges Kunstler's avatar
Georges Kunstler committed
105
##                               stringsAsFactors = FALSE)
106
107
108
data.all.sample <- read.csv(file = file.path(base.dir, "data.all.csv"),
                       stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
Georges Kunstler's avatar
Georges Kunstler committed
109
classes[classes=='integer'] <- "numeric"
110
111
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',
                           intern = TRUE))
Georges Kunstler's avatar
Georges Kunstler committed
112

113
114
115
116
data.tree.tot <- read.csv(file = file.path(base.dir, "data.all.csv"),
                     stringsAsFactors = FALSE,
                     nrows = nrows,
                     colClasses = classes)
Georges Kunstler's avatar
Georges Kunstler committed
117

118
    fun.data.for.lmer(data.tree.tot, trait, type.filling = type.filling)
Georges Kunstler's avatar
Georges Kunstler committed
119
120
}

121
122
123
124
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){
Georges Kunstler's avatar
Georges Kunstler committed
125
## remove tree with NA
126
data.tree <- subset(data.tree, subset = (!is.na(data.tree[["BA.G"]])) &
Georges Kunstler's avatar
Georges Kunstler committed
127
128
                                     (!is.na(data.tree[["D"]])) )
## remove tree with zero
129
data.tree <- subset(data.tree, subset = data.tree[["BA.G"]]>min.BA.G &
130
                    data.tree[["BA.G"]]<max.BA.G &
Georges Kunstler's avatar
Georges Kunstler committed
131
                    data.tree[["D"]]>9.9)
Georges Kunstler's avatar
Georges Kunstler committed
132
133
134
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
                    !is.na(data.tree[[BATOT]]))
135
data.tree <- data.tree[select.temp, ]
Georges Kunstler's avatar
Georges Kunstler committed
136
# select  obs abov min perc.neigh
137
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh]] > min.perc.neigh &
138
                    !is.na(data.tree[[perc.neigh]]))
Georges Kunstler's avatar
Georges Kunstler committed
139
# select  species with minimum obs
140
data.tree <- subset(data.tree, subset = data.tree[["sp"]] %in%
141
142
                    names(table(factor(data.tree[["sp"]])))[
                                   table(factor(data.tree[["sp"]]))>min.obs])
Georges Kunstler's avatar
Georges Kunstler committed
143
144
145
146
147
148
149
150
return(data.tree)
}

fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}

### get variables for lmer
151
fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn,
Georges Kunstler's avatar
Georges Kunstler committed
152
                                                   abs.CWM.tntf, tf,
153
                                                   min.BA.G = 60){
Georges Kunstler's avatar
Georges Kunstler committed
154
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
155
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
156
157
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
158
159
160
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
Georges Kunstler's avatar
Georges Kunstler committed
161

162
163
164
165
166
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])

 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
167
168
169
                   data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
170
sumTnTfBn.diff <- sumTnBn-sumTfBn
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
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)))
Georges Kunstler's avatar
Georges Kunstler committed
187
188
}

Georges Kunstler's avatar
Georges Kunstler committed
189

Georges Kunstler's avatar
Georges Kunstler committed
190
191
192
193
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
194
195
196
197
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"))
Georges Kunstler's avatar
Georges Kunstler committed
198
199
    stop("need trait to be in SLA Leaf.N Seed.mass
         SLA Wood.density or  Max.height ")
Georges Kunstler's avatar
Georges Kunstler committed
200
# get vars names
201
202
203
CWM.tn <- paste(trait, "CWM", 'fill', sep = ".")
abs.CWM.tntf <- paste(trait, "abs.CWM", 'fill', sep = ".")
tf <- paste(trait, "focal", sep = ".")
Georges Kunstler's avatar
Georges Kunstler committed
204
BATOT <- "BATOT"
205
perc.neigh <- paste(trait, "perc", type.filling, sep = ".")
Georges Kunstler's avatar
Georges Kunstler committed
206
207
208
209
210
211
212
213
214
215
216
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)
Georges Kunstler's avatar
Georges Kunstler committed
217
218
219
220
221
    return(lmer.data)
}