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

Georges Kunstler's avatar
Georges Kunstler committed
6
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
kunstler's avatar
kunstler committed
7
                                                 data.type = 'simple',
8
                                                 ...){
Georges Kunstler's avatar
Georges Kunstler committed
9
  for (m in model.files)
Georges Kunstler's avatar
Georges Kunstler committed
10
    (run.model.for.set.one.trait (m, fun.model, trait,
kunstler's avatar
kunstler committed
11
12
                                  data.type = data.type,
                                   ...))
Georges Kunstler's avatar
Georges Kunstler committed
13
14
15
}


Georges Kunstler's avatar
Georges Kunstler committed
16
run.model.for.set.one.trait <- function(model.file, fun.model, trait,
kunstler's avatar
kunstler committed
17
                                        data.type, ...){
Georges Kunstler's avatar
Georges Kunstler committed
18
    fun.model <- match.fun(fun.model)
kunstler's avatar
kunstler committed
19
20
21
    try(fun.model(model.file, trait,
                  data.type = data.type,
                   ...))
Georges Kunstler's avatar
Georges Kunstler committed
22
23
24
25
}


#=====================================================================
kunstler's avatar
kunstler committed
26
# Run lmer() (in package lme4) for one trait, one model
Georges Kunstler's avatar
Georges Kunstler committed
27
#=====================================================================
28

kunstler's avatar
kunstler committed
29
path.model <- "R/analysis/model.lmer"
kunstler's avatar
kunstler committed
30
31
32
model.files.lmer.Tf.0 <- file.path(path.model,
      c("model.lmer.LOGLIN.r.set.species.R",
        "model.lmer.LOGLIN.r.set.fixed.biomes.species.R"))
Kunstler Georges's avatar
Kunstler Georges committed
33
34
35
model.files.lmer.Tf.0b <- file.path(path.model,
      c("model.lmer.LOGLIN.MAT.MAP.r.set.species.R",
        "model.lmer.LOGLIN.MAT.MAP.r.set.fixed.biomes.species.R"))
kunstler's avatar
kunstler committed
36
37
38
model.files.lmer.Tf.1 <- file.path(path.model,
      c("model.lmer.LOGLIN.ER.AD.Tf.r.set.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.r.set.fixed.biomes.species.R"))
39
model.files.lmer.Tf.2 <- file.path(path.model,
kunstler's avatar
kunstler committed
40
41
      c("model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.fixed.biomes.species.R"))
kunstler's avatar
kunstler committed
42
model.files.lmer.Tf.3 <- file.path(path.model,
kunstler's avatar
kunstler committed
43
44
      c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.fixed.biomes.species.R"))
Kunstler Georges's avatar
Kunstler Georges committed
45
46
47
model.files.lmer.Tf.4 <- file.path(path.model,
      c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.fixed.biomes.species.R"))
Kunstler Georges's avatar
Kunstler Georges committed
48
49
50
model.files.lmer.Tf.intra.0 <- file.path(path.model,
      c("model.lmer.LOGLIN.intra.r.set.species.R",
        "model.lmer.LOGLIN.intra.r.set.fixed.biomes.species.R"))
kunstler's avatar
kunstler committed
51
52
53
54
55
56
model.files.lmer.Tf.intra.1 <- file.path(path.model,
      c("model.lmer.LOGLIN.ER.AD.Tf.intra.r.set.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.intra.r.set.fixed.biomes.species.R"))
model.files.lmer.Tf.intra.2 <- file.path(path.model,
      c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species.R",
        "model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species.R"))
kunstler's avatar
kunstler committed
57

Georges Kunstler's avatar
Georges Kunstler committed
58
59


kunstler's avatar
kunstler committed
60
61
call.lmer.and.save <- function(formula, df.lmer, path.out,
                               var.sample, ecocode.var, select.biome, list.sd){
Georges Kunstler's avatar
Georges Kunstler committed
62
    lmer.output <- lmer(formula = formula, data = df.lmer, REML = FALSE,
kunstler's avatar
kunstler committed
63
                        control = lmerControl(optCtrl = list(maxfun = 40000) ) )
Georges Kunstler's avatar
Georges Kunstler committed
64
    print(warnings())
65
66
67
68
    print(summary(lmer.output))
    relgrad <- with(lmer.output@optinfo$derivs,solve(Hessian,gradient))
    print('test convergence of relative gradient of hessian')
    print(max(abs(relgrad)) )
69
    if(is.na(select.biome)){
kunstler's avatar
kunstler committed
70
71
       name.file <- paste(var.sample,ecocode.var,
                          "results.nolog.all.rds", sep ='.')
Georges Kunstler's avatar
Georges Kunstler committed
72
       }
73
    if(!is.na(select.biome)){
kunstler's avatar
kunstler committed
74
75
76
       name.file <- paste(var.sample, ecocode.var,
                          gsub(' ', '_',select.biome),
                          "results.nolog.all.rds", sep ='.')
Georges Kunstler's avatar
Georges Kunstler committed
77
       }
78
  saveRDS(list(output = lmer.output, relgrad = relgrad, list.sd = list.sd),
Georges Kunstler's avatar
Georges Kunstler committed
79
          file = file.path(path.out, name.file))
Georges Kunstler's avatar
Georges Kunstler committed
80
81
}

82
83


Georges Kunstler's avatar
Georges Kunstler committed
84
85
run.lmer <- function (model.file, trait,
                      min.obs = 10, sample.size = NA,
Kunstler Georges's avatar
Kunstler Georges committed
86
                      ecocode.var = 'koppen',
kunstler's avatar
kunstler committed
87
                      data.type ='simple',
kunstler's avatar
kunstler committed
88
                      var.sample = NA,
Georges Kunstler's avatar
Georges Kunstler committed
89
                      select.biome = NA,
kunstler's avatar
kunstler committed
90
91
                      merge.biomes.TF = FALSE,
                      sample.vec.TF = FALSE) {
Georges Kunstler's avatar
Georges Kunstler committed
92
93
94
95
    require(lme4)
    source(model.file, local = TRUE)
    model <- load.model()
    #= Path for output
kunstler's avatar
kunstler committed
96
    path.out <- output.dir('lmer', model$name, trait, data.type)
97
    print(path.out)
Georges Kunstler's avatar
Georges Kunstler committed
98
99
100
    dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
    cat("run lmer for model", model.file, " for trait",
         trait, "\n")
kunstler's avatar
kunstler committed
101

kunstler's avatar
kunstler committed
102
103
    list.df.lmer <- load.data.for.lmer(trait,
                                  data.type = data.type,
Georges Kunstler's avatar
Georges Kunstler committed
104
                                  sample.size. = sample.size,
kunstler's avatar
kunstler committed
105
                                  ecocode.var. = ecocode.var,
Georges Kunstler's avatar
Georges Kunstler committed
106
                                  var.sample. = var.sample,
Georges Kunstler's avatar
Georges Kunstler committed
107
                                  select.biome. = select.biome,
kunstler's avatar
kunstler committed
108
109
                                  merge.biomes.TF = merge.biomes.TF,
                                  sample.vec.TF. = sample.vec.TF)
kunstler's avatar
kunstler committed
110
111
                                        # return a list of DF
     cat("Ok data with Nobs", nrow(list.df.lmer[[1]]),
Georges Kunstler's avatar
Georges Kunstler committed
112
         "\n")
kunstler's avatar
kunstler committed
113
114
    print(names(list.df.lmer[[1]]))
    print(dim(list.df.lmer[[1]]))
Georges Kunstler's avatar
Georges Kunstler committed
115
      #= Run model
kunstler's avatar
kunstler committed
116
117
118
119
120
121
        call.lmer.and.save(formula = model$lmer.formula.tree.id,
                           list.df.lmer[[1]], path.out = path.out,
                           var.sample = var.sample,
                           ecocode.var = ecocode.var,
                           select.biome = select.biome,
                           list.sd = list.df.lmer[[2]])
Georges Kunstler's avatar
Georges Kunstler committed
122
123
 gc()
}
Georges Kunstler's avatar
Georges Kunstler committed
124
125
126



Georges Kunstler's avatar
Georges Kunstler committed
127
#========================================================================
128
129
output.dir <- function (type , model, trait, set) {
  file.path("output", type, set, trait,  model)
Georges Kunstler's avatar
Georges Kunstler committed
130
131
132
133
134
135
}


#============================================================
# Function to prepare data for lmer
#============================================================
136
137
138
139

## add sample equivalent per ecocode
add.sampling.prob.by.var.sample <-  function(df, var.sample){
    vec.prob <-  1/table(df[[var.sample]])
140
141
    df1 <- data.frame(ecocode = names(vec.prob),
                      prob.sample = as.numeric(vec.prob))
142
    names(df1) <-  c(var.sample, 'prob.sample')
kunstler's avatar
kunstler committed
143
    df2 <- left_join(df, df1, by = var.sample)
144
    return(df2)
145
146
}

147
148
resample.plot.by.var <- function(df, var.sample, sample.size){
    df2 <- add.sampling.prob.by.var.sample(df, var.sample)
kunstler's avatar
kunstler committed
149
    df1 <- df2 %>% group_by(plot) %>% summarise(prob.sample = mean(prob.sample))
150
151
152
153
154
    num.sample.plot <- sample(1:nrow(df1), size = sample.size,
                              prob = df1$prob.sample)
    vec.sample.plot <- df1[['plot']][num.sample.plot]
    return(vec.sample.plot)
}
155

156

157

kunstler's avatar
kunstler committed
158
load.data.for.lmer <- function(trait, data.type,
159
                               base.dir = "output/processed",
Georges Kunstler's avatar
Georges Kunstler committed
160
                               sample.size. ,
kunstler's avatar
kunstler committed
161
                               ecocode.var.,
Kunstler Georges's avatar
Kunstler Georges committed
162
                               var.sample. = 'koppen',
163
                               select.biome. = NA,
164
                               select.set. = NA,
165
                               sample.vec.TF. = FALSE,
kunstler's avatar
kunstler committed
166
                               merge.biomes.TF = FALSE){
kunstler's avatar
kunstler committed
167
  require(dplyr)
kunstler's avatar
kunstler committed
168
169
  if (! data.type %in% c('simple', 'intra',  'all.census'))
      stop('not good data.type')
170
  df <- readRDS( file.path(base.dir,paste('data', trait, data.type, 'rds',
kunstler's avatar
kunstler committed
171
                                        sep = '.')))
kunstler's avatar
kunstler committed
172
  df <- dplyr::mutate(df,sp.name = ifelse(is.na(sp.name), 'missing.sp',
kunstler's avatar
kunstler committed
173
                                   sp.name))
174

kunstler's avatar
kunstler committed
175
  if(merge.biomes.TF) df <- merge.biomes(df)
176
  if(!is.na(select.biome.)) {
kunstler's avatar
kunstler committed
177
    df <-  dplyr::filter(df, biomes == select.biome.)
178
179
    }
  if(!is.na(select.set.)) {
kunstler's avatar
kunstler committed
180
    df <-  dplyr::filter(df, set == select.set.)
181
    }
kunstler's avatar
kunstler committed
182
  if(!is.na(sample.size.)){
kunstler's avatar
kunstler committed
183
   if(sample.size. < length(unique(df$plot))){
184
    if(!sample.vec.TF.){
185
186
     if(sample.size. > length(unique(df$plot))){
         sample.size. <- length(unique(df$plot))}
187
188
189
190
         plot.selected <- resample.plot.by.var(df, var.sample., sample.size.)
         sample.vec <- (1:nrow(df))[df$plot %in% plot.selected]
         df <- df[sample.vec, ]
         saveRDS(sample.vec,file = file.path(base.dir,paste('sample.vec',
kunstler's avatar
kunstler committed
191
                                         trait, data.type,'rds',
192
                                         sep = '.')))
193
194
         print(paste('sub-sampled by',var.sample.))
     }else{
kunstler's avatar
kunstler committed
195
196
197
198
199
         sample.vec <- readRDS(file =
                     file.path(base.dir,paste('sample.vec',
                               trait,
                               data.type,'rds',
                               sep = '.')))
200
         df <- df[sample.vec, ]
Kunstler Georges's avatar
Kunstler Georges committed
201
202
         print(paste('sub-sampled by',var.sample.,
                     'by loading sample.vec', nrow(df)))
203
     }
204
   }
205
206
207
208
  }
  if(data.type != 'all.census') {
    df <-  select.one.census.per.plot(df)
  }
kunstler's avatar
kunstler committed
209
  list.sd <- get.sd.lmer(df, trait, min.obs = 10)
kunstler's avatar
kunstler committed
210
  print('sd ok')
kunstler's avatar
kunstler committed
211
212
  res <- format.data.for.lmer(df, trait,
                              data.type = data.type,
kunstler's avatar
kunstler committed
213
                              ecocode.var = ecocode.var.)
214
return( list(res, list.sd))
Georges Kunstler's avatar
Georges Kunstler committed
215
}
216

kunstler's avatar
kunstler committed
217
merge.biomes <- function(data){
Georges Kunstler's avatar
Georges Kunstler committed
218
219
220
221
data$biomes <-  factor(data$biomes)
levels.name <- levels(data$biomes)
levels.name[levels.name == 'Tundra'] <- 'Boreal forest'
levels.name[levels.name == 'Subtropical desert'] <- 'Temperate grassland desert'
kunstler's avatar
kunstler committed
222
levels.name[levels.name == "Tropical rain forest"] <- "Tropical forest savanna"
Georges Kunstler's avatar
Georges Kunstler committed
223
levels(data$biomes) <- levels.name
kunstler's avatar
kunstler committed
224
print('biomes merged')
Georges Kunstler's avatar
Georges Kunstler committed
225
226
227
return(data)
}

kunstler's avatar
kunstler committed
228
229
select.one.census.per.plot <- function(df){
    require(dplyr)
230
    cat(dim(df))
kunstler's avatar
kunstler committed
231
232
233
234
235
236
237
    df <- df %>% mutate(plot.c = paste(plot, census))
    df.select <- df %>%
                 sample_n(., nrow(df)) %>% group_by( plot) %>%
                 summarise( plot.cs = first(plot.c)) %>%
                 dplyr::select(plot.cs)
    plots.select <- as.vector(df.select$plot.cs)
    df2 <-  filter(df, plot.c %in% plots.select)
238
239
    print(dim(df2))
    if(!(nrow(df2) < nrow(df))) stop('selection census error')
kunstler's avatar
kunstler committed
240
241
242
return(df2)
}

243

Georges Kunstler's avatar
Georges Kunstler committed
244

kunstler's avatar
kunstler committed
245
select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp,
Georges Kunstler's avatar
Georges Kunstler committed
246
247
248
                                         perc.neigh.gs, BATOT, min.obs,
                                         min.perc.neigh.sp = 0.90,
                                         min.perc.neigh.gs = 0.95,
kunstler's avatar
kunstler committed
249
250
                                         min.BA.G = -40,
                                         max.BA.G = 180){
Georges Kunstler's avatar
Georges Kunstler committed
251
252
253

select.temp <-(1:nrow(data.tree))[!is.na(data.tree[["BA.G"]])      &
                                  !is.na(data.tree[["D"]])         &
kunstler's avatar
kunstler committed
254
                                  !is.na(data.tree[["biomes"]])         &
kunstler's avatar
kunstler committed
255
256
                                  !is.na(data.tree[["MAT"]])         &
                                  !is.na(data.tree[["MAP"]])         &
Georges Kunstler's avatar
Georges Kunstler committed
257
258
259
260
261
262
263
264
265
266
267
268
                                 data.tree[["BA.G"]]>min.BA.G       &
                                 data.tree[["BA.G"]]<max.BA.G       &
                                 data.tree[["D"]]>9.9               &
                                 !is.na(data.tree[[abs.CWM.tntf]])  &
                                 !is.na(data.tree[[BATOT]])         &
                                 data.tree[[perc.neigh.sp]] >
                                      min.perc.neigh.sp             &
                                 !is.na(data.tree[[perc.neigh.sp]]) &
                                 data.tree[[perc.neigh.gs]] >
                                      min.perc.neigh.gs             &
                                 !is.na(data.tree[[perc.neigh.gs]])]

Georges Kunstler's avatar
Georges Kunstler committed
269
270

## remove tree with NA
Georges Kunstler's avatar
Georges Kunstler committed
271
272
273
274
275
276
277
278
279
280
data.tree <- data.tree[select.temp ,]
# select  species with minimum obs
select.temp <- (1:nrow(data.tree))[data.tree[["sp"]] %in%
                                   names(table(factor(data.tree[["sp"]])))[
                                     table(factor(data.tree[["sp"]]))>min.obs]]
data.tree <- data.tree[select.temp, ]
DT <-  drop(data.tree)
rm(data.tree)
gc()
return(DT)
Georges Kunstler's avatar
Georges Kunstler committed
281
282
}

kunstler's avatar
kunstler committed
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
select.data.trait <-  function(data.tree, trait,  min.obs = 10) {
if(! trait %in%  c("SLA", "Leaf.N", "Seed.mass",
                   "Wood.density", "Max.height"))
    stop("trait need to be in SLA, Leaf.N, Seed.mass,
         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.sp <- paste(trait, "perc", 'species', sep = ".")
perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
data.tree <- select.data.for.analysis(data.tree,
                                      abs.CWM.tntf,
                                      perc.neigh.sp,
                                      perc.neigh.gs,
                                      BATOT,
                                      min.obs)
    return(data.tree)
}
kunstler's avatar
kunstler committed
307

Georges Kunstler's avatar
Georges Kunstler committed
308

kunstler's avatar
kunstler committed
309
310
311
312
load.and.save.data.for.lmer <- function(trait,
                                        min.obs= 10,
                                        data.type = 'simple',
                                        base.dir = "output/processed"){
kunstler's avatar
kunstler committed
313
fname <- 'data.all.no.log.all.census.rds'
kunstler's avatar
kunstler committed
314
315
316
317
data.tree.tot <- readRDS(file.path(base.dir, fname))
df <- select.data.trait(data.tree.tot, trait)
saveRDS(df,file = file.path(base.dir,paste('data', trait, data.type, 'rds',
                                           sep = '.')))
kunstler's avatar
kunstler committed
318
}
Georges Kunstler's avatar
Georges Kunstler committed
319

320

kunstler's avatar
kunstler committed
321
322
###################
## FORMAT DATA
323

kunstler's avatar
kunstler committed
324
325
326
### get variables for analysis three type simple cat multi
scale.d <- function(x, ...){
return(as.vector(scale(x, ...)))
327
328
}

kunstler's avatar
kunstler committed
329
330
331
scale.nc <- function(x, center = FALSE){
return(as.vector(scale(x,
                       center = center,
kunstler's avatar
kunstler committed
332
                       scale = sd(x, na.rm = TRUE))))
kunstler's avatar
kunstler committed
333
334
335
336
}



337

kunstler's avatar
kunstler committed
338
get.variables <- function(data.tree, BATOT, CWM.tn,
Kunstler Georges's avatar
Kunstler Georges committed
339
                          abs.CWM.tntf, tf, ecocode.var = 'koppen',
kunstler's avatar
kunstler committed
340
                         min.BA.G = 40, min.MAT = 10){
kunstler's avatar
kunstler committed
341
342
logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
logD <- scale.d(log(data.tree[["D"]]))
kunstler's avatar
kunstler committed
343
344
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
MAP <- scale.d(log(data.tree[["MAP"]]))
345
species.id <- as.character(factor(data.tree[["sp.name"]]))
Georges Kunstler's avatar
Georges Kunstler committed
346
plot.id <- as.character(factor(data.tree[["plot"]]))
kunstler's avatar
kunstler committed
347
tree.id <- as.character(factor(data.tree[["tree.id"]]))
Georges Kunstler's avatar
Georges Kunstler committed
348
349

set.id <- as.character(factor(data.tree[["set"]]))
Kunstler Georges's avatar
Kunstler Georges committed
350
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
351
biomes.id <- factor(data.tree[['biomes']])
Georges Kunstler's avatar
Georges Kunstler committed
352
353
354
355

 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
                   data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
356
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
357
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
358
sumBn <- data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
359
360
361
362
363
364
365
return(data.frame(logG =  logG,
                  logD = logD,
                  MAT = MAT,
                  MAP = MAP,
                  species.id = species.id,
                  set.id = set.id,
                  ecocode.id = ecocode.id,
366
                  biomes.id = biomes.id,
Georges Kunstler's avatar
Georges Kunstler committed
367
                  plot.id = plot.id,
kunstler's avatar
kunstler committed
368
                  tree.id = tree.id,
kunstler's avatar
kunstler committed
369
                  sumTnTfBn.abs = scale.nc(sumTnTfBn.abs, center = FALSE),
kunstler's avatar
kunstler committed
370
                  Tf = scale.d(data.tree[[tf]]),
kunstler's avatar
kunstler committed
371
372
373
                  sumTnBn = scale.nc(sumTnBn, center = FALSE),
                  sumTfBn = scale.nc(sumTfBn, center = FALSE),
                  sumBn = scale.nc(sumBn, center = FALSE),
Georges Kunstler's avatar
Georges Kunstler committed
374
375
376
                  stringsAsFactors = FALSE))
}

kunstler's avatar
kunstler committed
377

kunstler's avatar
kunstler committed
378
get.variables.intra <- function(data.tree, BATOT, CWM.tn,
Kunstler Georges's avatar
Kunstler Georges committed
379
                          abs.CWM.tntf, tf, ecocode.var = 'koppen',
kunstler's avatar
kunstler committed
380
                         min.BA.G = 40, min.MAT = 10){
kunstler's avatar
kunstler committed
381
382
logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
logD <- scale.d(log(data.tree[["D"]]))
kunstler's avatar
kunstler committed
383
384
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
MAP <- scale.d(log(data.tree[["MAP"]]))
kunstler's avatar
kunstler committed
385
386
species.id <- as.character(factor(data.tree[["sp.name"]]))
plot.id <- as.character(factor(data.tree[["plot"]]))
kunstler's avatar
kunstler committed
387
tree.id <- as.character(factor(data.tree[["tree.id"]]))
kunstler's avatar
kunstler committed
388

kunstler's avatar
kunstler committed
389
set.id <- as.character(factor(data.tree[["set"]]))
Kunstler Georges's avatar
Kunstler Georges committed
390
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
kunstler's avatar
kunstler committed
391
392
biomes.id <- factor(data.tree[['biomes']])

kunstler's avatar
kunstler committed
393
394
395
396
397
398
 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
                   (data.tree[[BATOT]] - data.tree[['BAintra']])
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumBn.inter <- (data.tree[[BATOT]] - data.tree[['BAintra']])
kunstler's avatar
kunstler committed
399
400
401
402
403
404
405
406
407
return(data.frame(logG =  logG,
                  logD = logD,
                  MAT = MAT,
                  MAP = MAP,
                  species.id = species.id,
                  set.id = set.id,
                  ecocode.id = ecocode.id,
                  biomes.id = biomes.id,
                  plot.id = plot.id,
kunstler's avatar
kunstler committed
408
                  tree.id = tree.id,
kunstler's avatar
kunstler committed
409
410
411
412
413
414
415
                  sumTnTfBn.abs = scale.nc(sumTnTfBn.abs, center = FALSE),
                  Tf = scale.d(data.tree[[tf]]),
                  sumTnBn = scale.nc(sumTnBn, center = FALSE),
                  sumTfBn = scale.nc(sumTfBn, center = FALSE),
                  sumBn.inter = scale.nc(sumBn.inter, center = FALSE),
                  sumBn.intra = scale.nc(data.tree[['BAintra']],
                                         center = FALSE),
kunstler's avatar
kunstler committed
416
417
418
419
                  stringsAsFactors = FALSE))
}


420
421

#============================================================
kunstler's avatar
kunstler committed
422
# Function to prepare data for lmer with CWM data
423
#============================================================
kunstler's avatar
kunstler committed
424
format.data.for.lmer <-  function(data.tree, trait, min.obs = 10,
Kunstler Georges's avatar
Kunstler Georges committed
425
                                  data.type = 'simple', ecocode.var = 'koppen') {
426
if(! trait %in%  c("SLA", "Leaf.N", "Seed.mass",
kunstler's avatar
kunstler committed
427
                   "SLA", "Wood.density", "Max.height"))
428
    stop("need trait to be in SLA Leaf.N Seed.mass
kunstler's avatar
kunstler committed
429
         SLA Wood.density or  Max.height")
430
431
432
433
434
# 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"
Georges Kunstler's avatar
Georges Kunstler committed
435
#= DATA LIST FOR LMER
436
if(data.type == 'simple' | data.type == 'all.census') {
kunstler's avatar
kunstler committed
437
    print('simple')
kunstler's avatar
kunstler committed
438
439
440
441
442
443
    lmer.data <- get.variables(data.tree,
                               BATOT,
                               CWM.tn,
                               abs.CWM.tntf,
                               tf,
                               ecocode.var)
Georges Kunstler's avatar
Georges Kunstler committed
444
}
kunstler's avatar
kunstler committed
445
446
447
if(data.type =='intra') {
    print('intra')
    lmer.data <- get.variables.intra(data.tree,
kunstler's avatar
kunstler committed
448
449
450
451
452
                               BATOT,
                               CWM.tn,
                               abs.CWM.tntf,
                               tf,
                               ecocode.var)
kunstler's avatar
kunstler committed
453
    }
Georges Kunstler's avatar
Georges Kunstler committed
454
455
456
    return(lmer.data)
}

kunstler's avatar
kunstler committed
457
458
## SD for plot
get.sd.lmer <-  function(data.tree, trait, min.obs = 10) {
459
if(! trait %in%  c("SLA", "Leaf.N", "Seed.mass",
kunstler's avatar
kunstler committed
460
                   "SLA", "Wood.density", "Max.height"))
461
    stop("need trait to be in SLA Leaf.N Seed.mass
kunstler's avatar
kunstler committed
462
         SLA Wood.density or  Max.height")
463
464
465
466
467
468
# 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"
#= DATA SD LIST FOR LMER
kunstler's avatar
kunstler committed
469
list.sd <- compute.sd.mean.var.lmer(data.tree, tf, CWM.tn,
kunstler's avatar
kunstler committed
470
                                    abs.CWM.tntf, min.BA.G = 40)
471
472
473
    return(list.sd)
}

kunstler's avatar
kunstler committed
474
compute.sd.mean.var.lmer <- function(data.tree, tf, CWM.tn,
kunstler's avatar
kunstler committed
475
                                         abs.CWM.tntf, min.BA.G = 40){
kunstler's avatar
kunstler committed
476
477
478
479
480
481
482
483
list.sd <- list(
'mean.logG' = mean(log(data.tree[["BA.G"]] + min.BA.G+1), na.rm = TRUE),
'sd.logG' = sd(log(data.tree[["BA.G"]] + min.BA.G+1), na.rm = TRUE),
'mean.logD' = mean(log(data.tree[["D"]] ), na.rm = TRUE),
'sd.logD' = sd(log(data.tree[["D"]] ), na.rm = TRUE),
'mean.tf' = mean(data.tree[[tf]] , na.rm = TRUE),
'sd.tf' = sd(data.tree[[tf]] , na.rm = TRUE),
'sd.sumBn' = sd( data.tree[['BATOT']], na.rm = TRUE),
kunstler's avatar
kunstler committed
484
485
486
'sd.sumBn.inter' = sd( data.tree[['BATOT']] - data.tree[['BAintra']],
    na.rm = TRUE),
'sd.sumBn.intra' = sd( data.tree[['BAintra']], na.rm = TRUE),
kunstler's avatar
kunstler committed
487
488
489
490
491
492
493
494
'sd.sumTfBn' = sd(data.tree[[tf]]* data.tree[['BATOT']], na.rm = TRUE),
'sd.sumTnTfBn.abs' = sd(data.tree[[abs.CWM.tntf]]*
                   data.tree[['BATOT']], na.rm = TRUE),
'sd.sumTnBn' = sd(data.tree[[CWM.tn]]*data.tree[['BATOT']], na.rm = TRUE)
)
return(list.sd)
}

495

Georges Kunstler's avatar
Georges Kunstler committed
496