lmer.run.R 19.3 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's avatar
kunstler committed
33
34
35
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"))
36
model.files.lmer.Tf.2 <- file.path(path.model,
kunstler's avatar
kunstler committed
37
38
      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
39
model.files.lmer.Tf.3 <- file.path(path.model,
kunstler's avatar
kunstler committed
40
41
      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
42
43
44
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's avatar
kunstler committed
45
46
47
48
49
50
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
51

Georges Kunstler's avatar
Georges Kunstler committed
52
53


kunstler's avatar
kunstler committed
54
55
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
56
    lmer.output <- lmer(formula = formula, data = df.lmer, REML = FALSE,
kunstler's avatar
kunstler committed
57
                        control = lmerControl(optCtrl = list(maxfun = 40000) ) )
Georges Kunstler's avatar
Georges Kunstler committed
58
    print(warnings())
59
60
61
62
    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)) )
63
    if(is.na(select.biome)){
kunstler's avatar
kunstler committed
64
65
       name.file <- paste(var.sample,ecocode.var,
                          "results.nolog.all.rds", sep ='.')
Georges Kunstler's avatar
Georges Kunstler committed
66
       }
67
    if(!is.na(select.biome)){
kunstler's avatar
kunstler committed
68
69
70
       name.file <- paste(var.sample, ecocode.var,
                          gsub(' ', '_',select.biome),
                          "results.nolog.all.rds", sep ='.')
Georges Kunstler's avatar
Georges Kunstler committed
71
       }
72
  saveRDS(list(output = lmer.output, relgrad = relgrad, list.sd = list.sd),
Georges Kunstler's avatar
Georges Kunstler committed
73
          file = file.path(path.out, name.file))
Georges Kunstler's avatar
Georges Kunstler committed
74
75
}

76
77


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

kunstler's avatar
kunstler committed
96
97
    list.df.lmer <- load.data.for.lmer(trait,
                                  data.type = data.type,
Georges Kunstler's avatar
Georges Kunstler committed
98
                                  sample.size. = sample.size,
kunstler's avatar
kunstler committed
99
                                  ecocode.var. = ecocode.var,
Georges Kunstler's avatar
Georges Kunstler committed
100
                                  var.sample. = var.sample,
Georges Kunstler's avatar
Georges Kunstler committed
101
                                  select.biome. = select.biome,
kunstler's avatar
kunstler committed
102
103
                                  merge.biomes.TF = merge.biomes.TF,
                                  sample.vec.TF. = sample.vec.TF)
kunstler's avatar
kunstler committed
104
105
                                        # return a list of DF
     cat("Ok data with Nobs", nrow(list.df.lmer[[1]]),
Georges Kunstler's avatar
Georges Kunstler committed
106
         "\n")
kunstler's avatar
kunstler committed
107
108
    print(names(list.df.lmer[[1]]))
    print(dim(list.df.lmer[[1]]))
Georges Kunstler's avatar
Georges Kunstler committed
109
      #= Run model
kunstler's avatar
kunstler committed
110
111
112
113
114
115
        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
116
117
 gc()
}
Georges Kunstler's avatar
Georges Kunstler committed
118
119
120



Georges Kunstler's avatar
Georges Kunstler committed
121
#========================================================================
122
123
output.dir <- function (type , model, trait, set) {
  file.path("output", type, set, trait,  model)
Georges Kunstler's avatar
Georges Kunstler committed
124
125
126
127
128
129
}


#============================================================
# Function to prepare data for lmer
#============================================================
130
131
132
133

## add sample equivalent per ecocode
add.sampling.prob.by.var.sample <-  function(df, var.sample){
    vec.prob <-  1/table(df[[var.sample]])
134
135
    df1 <- data.frame(ecocode = names(vec.prob),
                      prob.sample = as.numeric(vec.prob))
136
    names(df1) <-  c(var.sample, 'prob.sample')
kunstler's avatar
kunstler committed
137
    df2 <- left_join(df, df1, by = var.sample)
138
    return(df2)
139
140
}

141
142
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
143
    df1 <- df2 %>% group_by(plot) %>% summarise(prob.sample = mean(prob.sample))
144
145
146
147
148
    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)
}
149

150

151

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

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

kunstler's avatar
kunstler committed
210
merge.biomes <- function(data){
Georges Kunstler's avatar
Georges Kunstler committed
211
212
213
214
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
215
levels.name[levels.name == "Tropical rain forest"] <- "Tropical forest savanna"
Georges Kunstler's avatar
Georges Kunstler committed
216
levels(data$biomes) <- levels.name
kunstler's avatar
kunstler committed
217
print('biomes merged')
Georges Kunstler's avatar
Georges Kunstler committed
218
219
220
return(data)
}

kunstler's avatar
kunstler committed
221
222
select.one.census.per.plot <- function(df){
    require(dplyr)
223
    cat(dim(df))
kunstler's avatar
kunstler committed
224
225
226
227
228
229
230
    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)
231
232
    print(dim(df2))
    if(!(nrow(df2) < nrow(df))) stop('selection census error')
kunstler's avatar
kunstler committed
233
234
235
return(df2)
}

236

Georges Kunstler's avatar
Georges Kunstler committed
237

kunstler's avatar
kunstler committed
238
select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp,
Georges Kunstler's avatar
Georges Kunstler committed
239
240
241
                                         perc.neigh.gs, BATOT, min.obs,
                                         min.perc.neigh.sp = 0.90,
                                         min.perc.neigh.gs = 0.95,
kunstler's avatar
kunstler committed
242
243
                                         min.BA.G = -40,
                                         max.BA.G = 180){
Georges Kunstler's avatar
Georges Kunstler committed
244
245
246

select.temp <-(1:nrow(data.tree))[!is.na(data.tree[["BA.G"]])      &
                                  !is.na(data.tree[["D"]])         &
kunstler's avatar
kunstler committed
247
                                  !is.na(data.tree[["biomes"]])         &
kunstler's avatar
kunstler committed
248
249
                                  !is.na(data.tree[["MAT"]])         &
                                  !is.na(data.tree[["MAP"]])         &
Georges Kunstler's avatar
Georges Kunstler committed
250
251
252
253
254
255
256
257
258
259
260
261
                                 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
262
263

## remove tree with NA
Georges Kunstler's avatar
Georges Kunstler committed
264
265
266
267
268
269
270
271
272
273
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
274
275
}

kunstler's avatar
kunstler committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
#============================================================
# 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
300

Georges Kunstler's avatar
Georges Kunstler committed
301

kunstler's avatar
kunstler committed
302
303
304
305
load.and.save.data.for.lmer <- function(trait,
                                        min.obs= 10,
                                        data.type = 'simple',
                                        base.dir = "output/processed"){
kunstler's avatar
kunstler committed
306
fname <- 'data.all.no.log.all.census.rds'
kunstler's avatar
kunstler committed
307
308
309
310
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
311
}
Georges Kunstler's avatar
Georges Kunstler committed
312

313

kunstler's avatar
kunstler committed
314
315
###################
## FORMAT DATA
316

kunstler's avatar
kunstler committed
317
318
319
### get variables for analysis three type simple cat multi
scale.d <- function(x, ...){
return(as.vector(scale(x, ...)))
320
321
}

kunstler's avatar
kunstler committed
322
323
324
scale.nc <- function(x, center = FALSE){
return(as.vector(scale(x,
                       center = center,
kunstler's avatar
kunstler committed
325
                       scale = sd(x, na.rm = TRUE))))
kunstler's avatar
kunstler committed
326
327
328
329
}



330

kunstler's avatar
kunstler committed
331
332
get.variables <- function(data.tree, BATOT, CWM.tn,
                          abs.CWM.tntf, tf, ecocode.var = 'wwf',
kunstler's avatar
kunstler committed
333
                         min.BA.G = 40, min.MAT = 10){
kunstler's avatar
kunstler committed
334
335
logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
logD <- scale.d(log(data.tree[["D"]]))
kunstler's avatar
kunstler committed
336
337
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
MAP <- scale.d(log(data.tree[["MAP"]]))
338
species.id <- as.character(factor(data.tree[["sp.name"]]))
Georges Kunstler's avatar
Georges Kunstler committed
339
plot.id <- as.character(factor(data.tree[["plot"]]))
kunstler's avatar
kunstler committed
340
tree.id <- as.character(factor(data.tree[["tree.id"]]))
Georges Kunstler's avatar
Georges Kunstler committed
341
342

set.id <- as.character(factor(data.tree[["set"]]))
Kunstler Georges's avatar
Kunstler Georges committed
343
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
344
biomes.id <- factor(data.tree[['biomes']])
Georges Kunstler's avatar
Georges Kunstler committed
345
346
347
348

 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
                   data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
349
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
350
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
351
sumBn <- data.tree[[BATOT]]
Georges Kunstler's avatar
Georges Kunstler committed
352
353
354
355
356
357
358
return(data.frame(logG =  logG,
                  logD = logD,
                  MAT = MAT,
                  MAP = MAP,
                  species.id = species.id,
                  set.id = set.id,
                  ecocode.id = ecocode.id,
359
                  biomes.id = biomes.id,
Georges Kunstler's avatar
Georges Kunstler committed
360
                  plot.id = plot.id,
kunstler's avatar
kunstler committed
361
                  tree.id = tree.id,
kunstler's avatar
kunstler committed
362
                  sumTnTfBn.abs = scale.nc(sumTnTfBn.abs, center = FALSE),
kunstler's avatar
kunstler committed
363
                  Tf = scale.d(data.tree[[tf]]),
kunstler's avatar
kunstler committed
364
365
366
                  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
367
368
369
                  stringsAsFactors = FALSE))
}

kunstler's avatar
kunstler committed
370

kunstler's avatar
kunstler committed
371
372
373
get.variables.intra <- function(data.tree, BATOT, CWM.tn,
                          abs.CWM.tntf, tf, ecocode.var = 'wwf',
                         min.BA.G = 40, min.MAT = 10){
kunstler's avatar
kunstler committed
374
375
logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
logD <- scale.d(log(data.tree[["D"]]))
kunstler's avatar
kunstler committed
376
377
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
MAP <- scale.d(log(data.tree[["MAP"]]))
kunstler's avatar
kunstler committed
378
379
species.id <- as.character(factor(data.tree[["sp.name"]]))
plot.id <- as.character(factor(data.tree[["plot"]]))
kunstler's avatar
kunstler committed
380
tree.id <- as.character(factor(data.tree[["tree.id"]]))
kunstler's avatar
kunstler committed
381

kunstler's avatar
kunstler committed
382
set.id <- as.character(factor(data.tree[["set"]]))
Kunstler Georges's avatar
Kunstler Georges committed
383
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
kunstler's avatar
kunstler committed
384
385
biomes.id <- factor(data.tree[['biomes']])

kunstler's avatar
kunstler committed
386
387
388
389
390
391
 #= 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
392
393
394
395
396
397
398
399
400
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
401
                  tree.id = tree.id,
kunstler's avatar
kunstler committed
402
403
404
405
406
407
408
                  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
409
410
411
412
                  stringsAsFactors = FALSE))
}


413
414

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

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

kunstler's avatar
kunstler committed
467
compute.sd.mean.var.lmer <- function(data.tree, tf, CWM.tn,
kunstler's avatar
kunstler committed
468
                                         abs.CWM.tntf, min.BA.G = 40){
kunstler's avatar
kunstler committed
469
470
471
472
473
474
475
476
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
477
478
479
'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
480
481
482
483
484
485
486
487
'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)
}

488

Georges Kunstler's avatar
Georges Kunstler committed
489