lmer.run.nolog.all.R 13.3 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
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
Georges Kunstler's avatar
Georges Kunstler committed
15
                                                 type.filling, cat.TF, ...){
Georges Kunstler's avatar
Georges Kunstler committed
16
  for (m in model.files)
17
    (run.model.for.set.one.trait (m, fun.model, trait,
Georges Kunstler's avatar
Georges Kunstler committed
18
                                  type.filling = type.filling, cat.TF = cat.TF, ...))
Georges Kunstler's avatar
Georges Kunstler committed
19
20
21
}


22
run.model.for.set.one.trait <- function(model.file, fun.model, trait,
Georges Kunstler's avatar
Georges Kunstler committed
23
                                        type.filling, cat.TF, ...){
Georges Kunstler's avatar
Georges Kunstler committed
24
    fun.model <- match.fun(fun.model)
Georges Kunstler's avatar
Georges Kunstler committed
25
    try(fun.model(model.file, trait,  type.filling = type.filling,cat.TF = cat.TF, ...))
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
                      type.filling, cat.TF, fname = 'data.all.no.std.csv') {
Georges Kunstler's avatar
Georges Kunstler committed
66
67
68
69
    require(lme4)
    source(model.file, local = TRUE)
    model <- load.model()
    #= Path for output
Georges Kunstler's avatar
Georges Kunstler committed
70
71
72
    if(fname == 'data.all.no.std.csv') dir <- 'all.no.log'
    if(fname == 'data.all.csv') dir <- 'all.log'
    path.out <- output.dir.lmer(model$name, trait, dir,
73
                                type.filling = type.filling)
Georges Kunstler's avatar
Georges Kunstler committed
74
    print(path.out)
75
    dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
76
    cat("run lmer for model", model.file, " for trait",
77
         trait, "\n")
Georges Kunstler's avatar
Georges Kunstler committed
78
  df.lmer <- load.and.prepare.data.for.lmer(trait,
Georges Kunstler's avatar
Georges Kunstler committed
79
                                           min.obs, sample.size,
Georges Kunstler's avatar
Georges Kunstler committed
80
81
                                         type.filling = type.filling,
                                                fname = fname, cat.TF = cat.TF)
82
83
                                        # return a DF
     cat("Ok data with Nobs", nrow(df.lmer),
Georges Kunstler's avatar
Georges Kunstler committed
84
85
         "\n")
        #= Run model
86
87
        fun.call.lmer.and.save(formula = model$lmer.formula.tree.id,
                               df.lmer, path.out)
Georges Kunstler's avatar
Georges Kunstler committed
88
89
 rm(df.lmer)
 gc()
Georges Kunstler's avatar
Georges Kunstler committed
90
91
92
93
94
}


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

95
96
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
97
98
}

Georges Kunstler's avatar
Georges Kunstler committed
99
100
fun.load.data.all <- function(base.dir,fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
101
102
                       stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
Georges Kunstler's avatar
Georges Kunstler committed
103
classes[classes=='integer'] <- "numeric"
Georges Kunstler's avatar
Georges Kunstler committed
104
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
105
                           intern = TRUE))
Georges Kunstler's avatar
Georges Kunstler committed
106

Georges Kunstler's avatar
Georges Kunstler committed
107
data.tree.tot <- read.csv(file = file.path(base.dir, fname),
108
109
110
                     stringsAsFactors = FALSE,
                     nrows = nrows,
                     colClasses = classes)
Georges Kunstler's avatar
Georges Kunstler committed
111
112
return(data.tree.tot)
}
Georges Kunstler's avatar
Georges Kunstler committed
113

Georges Kunstler's avatar
Georges Kunstler committed
114
115
116
117
118
119
120
121
122
123
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait,
                                           min.obs, sample.size, type.filling,
                                           cat.TF,
                                           fname = 'data.all.no.std.csv',
                                           base.dir = "output/processed"){
 data.tree.tot <- fun.load.data.all(base.dir,fname)
 fun.data.for.lmer(data.tree.tot, trait, type.filling = type.filling, cat.TF = cat.TF)
Georges Kunstler's avatar
Georges Kunstler committed
124
125
}

126
127
fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp,
                                         perc.neigh.gs, BATOT, min.obs,
Georges Kunstler's avatar
Georges Kunstler committed
128
129
130
                                         min.perc.neigh.sp = 0.90,
                                         min.perc.neigh.gs = 0.95,
                                         min.BA.G = -60,
131
                                         max.BA.G = 500){
Georges Kunstler's avatar
Georges Kunstler committed
132
## remove tree with NA
133
data.tree <- subset(data.tree, subset = (!is.na(data.tree[["BA.G"]])) &
Georges Kunstler's avatar
Georges Kunstler committed
134
135
                                     (!is.na(data.tree[["D"]])) )
## remove tree with zero
136
data.tree <- subset(data.tree, subset = data.tree[["BA.G"]]>min.BA.G &
137
                    data.tree[["BA.G"]]<max.BA.G &
Georges Kunstler's avatar
Georges Kunstler committed
138
                    data.tree[["D"]]>9.9)
Georges Kunstler's avatar
Georges Kunstler committed
139
140
141
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
                    !is.na(data.tree[[BATOT]]))
142
data.tree <- data.tree[select.temp, ]
143
# select  obs abov min perc.neigh species
Georges Kunstler's avatar
Georges Kunstler committed
144
145
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
                    min.perc.neigh.sp &
146
147
                    !is.na(data.tree[[perc.neigh.sp]]))
# select  obs abov min perc.neigh genus
Georges Kunstler's avatar
Georges Kunstler committed
148
149
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
                    min.perc.neigh.gs &
150
                    !is.na(data.tree[[perc.neigh.gs]]))
Georges Kunstler's avatar
Georges Kunstler committed
151
# select  species with minimum obs
152
data.tree <- subset(data.tree, subset = data.tree[["sp"]] %in%
153
154
                    names(table(factor(data.tree[["sp"]])))[
                                   table(factor(data.tree[["sp"]]))>min.obs])
Georges Kunstler's avatar
Georges Kunstler committed
155
156
157
158
159
160
161
162
return(data.tree)
}

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

### get variables for lmer
163
fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn,
Georges Kunstler's avatar
Georges Kunstler committed
164
                                                   abs.CWM.tntf, tf,
165
                                                   min.BA.G = 60){
Georges Kunstler's avatar
Georges Kunstler committed
166
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]] + min.BA.G))
167
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
168
169
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
170
171
172
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
173

174
175
176
177
178
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])

 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
179
180
181
                   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
182
sumTnTfBn.diff <- sumTnBn-sumTfBn
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
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
199
200
}

Georges Kunstler's avatar
Georges Kunstler committed
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
fun.get.the.variables.for.lmer.tree.id.cat <- 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']])

#get the three cwm per cat
vec.CWM.tn <- paste(CWM.tn, c('A_EV', 'A_D', 'C'), sep = '.')
vec.abs.CWM.tntf <- paste(abs.CWM.tntf, c('A_EV', 'A_D', 'C'), sep = '.')

sumTnTfBn.abs.A_EV<- data.tree[[vec.abs.CWM.tntf[1]]]
sumTnTfBn.abs.A_D <- data.tree[[vec.abs.CWM.tntf[2]]]
sumTnTfBn.abs.C   <- data.tree[[vec.abs.CWM.tntf[3]]]

sumTnBn.A_EV<- data.tree[[vec.CWM.tn[1]]]
sumTnBn.A_D <- data.tree[[vec.CWM.tn[2]]]
sumTnBn.C   <- data.tree[[vec.CWM.tn[3]]]

sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *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.abs.A_EV =
                   fun.center.and.standardized.var(sumTnTfBn.abs.A_EV),
            sumTnTfBn.abs.A_D =
                   fun.center.and.standardized.var(sumTnTfBn.abs.A_D),
            sumTnTfBn.abs.C =
                   fun.center.and.standardized.var(sumTnTfBn.abs.C),
            Tf.A_EV= fun.center.and.standardized.var(data.tree[[tf]] *
                  as.numeric(data.tree[['SLA.cat']] == 1)),
            Tf.A_D= fun.center.and.standardized.var(data.tree[[tf]] *
                as.numeric(data.tree[['SLA.cat']] == 2)),
            Tf.C = fun.center.and.standardized.var(data.tree[[tf]] *
                as.numeric(data.tree[['SLA.cat']] == 3)),
            sumTnBn.A_EV= fun.center.and.standardized.var(sumTnBn.A_EV),
            sumTnBn.A_D= fun.center.and.standardized.var(sumTnBn.A_D),
            sumTnBn.C = fun.center.and.standardized.var(sumTnBn.C),
            sumTfBn.A_EV= fun.center.and.standardized.var(sumTfBn.A_EV),
            sumTfBn.A_D= fun.center.and.standardized.var(sumTfBn.A_D),
            sumTfBn.C = fun.center.and.standardized.var(sumTfBn.C),
            sumBn = fun.center.and.standardized.var(sumBn)))
}

Georges Kunstler's avatar
Georges Kunstler committed
263

Georges Kunstler's avatar
Georges Kunstler committed
264
265
266
267
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
268
fun.data.for.lmer <-  function(data.tree, trait, min.obs = 10,
Georges Kunstler's avatar
Georges Kunstler committed
269
                               type.filling = 'species', cat.TF) {
270
271
if(! trait %in%  c("SLA", "Leaf.N", "Seed.mass",
                   "SLA", "Wood.density", "Max.height"))
Georges Kunstler's avatar
Georges Kunstler committed
272
273
    stop("need trait to be in SLA Leaf.N Seed.mass
         SLA Wood.density or  Max.height ")
Georges Kunstler's avatar
Georges Kunstler committed
274
# get vars names
275
276
277
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
278
BATOT <- "BATOT"
279
280
perc.neigh.sp <- paste(trait, "perc", 'species', sep = ".")
perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
Georges Kunstler's avatar
Georges Kunstler committed
281
282
data.tree <- fun.select.data.for.analysis(data.tree,
                                          abs.CWM.tntf,
Georges Kunstler's avatar
Georges Kunstler committed
283
                                          perc.neigh.sp,
284
                                          perc.neigh.gs,
Georges Kunstler's avatar
Georges Kunstler committed
285
286
287
                                          BATOT,
                                          min.obs)
#= DATA LIST FOR LMER
Georges Kunstler's avatar
Georges Kunstler committed
288
289
290
291
292
293
if(!cat.TF) lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,
                                                    BATOT,
                                                    CWM.tn,
                                                    abs.CWM.tntf,
                                                    tf)
if(cat.TF) lmer.data <- fun.get.the.variables.for.lmer.tree.id.cat(data.tree,
Georges Kunstler's avatar
Georges Kunstler committed
294
295
296
297
                                                    BATOT,
                                                    CWM.tn,
                                                    abs.CWM.tntf,
                                                    tf)
Georges Kunstler's avatar
Georges Kunstler committed
298
299
300
301
302
    return(lmer.data)
}