lmer.run.R 9.98 KB
Newer Older
Georges Kunstler's avatar
Georges Kunstler committed
1
2
3
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION
Georges Kunstler's avatar
Georges Kunstler committed
4
library(lme4)
Georges Kunstler's avatar
Georges Kunstler committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32


get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
  sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
}

run.models.for.set.all.traits  <- function(set,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, set, type.filling=type.filling, ...)
}

run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set),  ...){
  for (m in model.files)
    try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling, ...))
}


run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set),  ...){
    fun.model <- match.fun(fun.model)
  for (e in ecoregions)
    try(fun.model(model.file, trait, set, e, type.filling=type.filling,...))
}


#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
Georges Kunstler's avatar
Georges Kunstler committed
33
model.files.lmer.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.R",
Georges Kunstler's avatar
Georges Kunstler committed
34
                 "R/analysis/model.lmer/model.lmer.LOGLIN.R.R",
Georges Kunstler's avatar
Georges Kunstler committed
35
36
                 "R/analysis/model.lmer/model.lmer.LOGLIN.ER.R")
model.files.lmer.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
Georges Kunstler's avatar
Georges Kunstler committed
37
                 "R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",
Georges Kunstler's avatar
Georges Kunstler committed
38
39
40
                 "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R")

model.files.lmer.Tf.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.Tf.R",
41
                 "R/analysis/model.lmer/model.lmer.LOGLIN.R.Tf.R",
Georges Kunstler's avatar
Georges Kunstler committed
42
43
44
                 "R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R")
model.files.lmer.Tf.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",
                 "R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",
45
                 "R/analysis/model.lmer/model.lmer.LOGLIN.HD.Tf.R",
Georges Kunstler's avatar
Georges Kunstler committed
46
                 "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R")
Georges Kunstler's avatar
Georges Kunstler committed
47
48
49
50
51
52


fun.test.if.multi.census <- function(data){
return("tree.id" %in% names(data))
}

Georges Kunstler's avatar
Georges Kunstler committed
53
54
55
56
57
58
fun.call.lmer <- function(formula,df.lmer){
   lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
   return(lmer.output)
}


Georges Kunstler's avatar
Georges Kunstler committed
59
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
60
   lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
61
   print(summary(lmer.output))
Georges Kunstler's avatar
Georges Kunstler committed
62
   saveRDS(lmer.output,file=file.path(path.out, "results.rds"))
Georges Kunstler's avatar
Georges Kunstler committed
63
64
65
66
67
68
69
70
71
72
73
}

run.lmer <- function (model.file, trait, set, ecoregion,
                      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, set,
                                ecoregion,type.filling=type.filling)
74
    print(path.out)
Georges Kunstler's avatar
Georges Kunstler committed
75
76
77
78
79
80
81
82
83
84
85
86
87
    dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
    cat("run lmer for model",model.file," for set",
         set,"ecoregion",ecoregion,"trait",
         trait,"\n")
      df.lmer <- load.and.prepare.data.for.lmer(trait, set, ecoregion,
                                           min.obs, sample.size,
                                         type.filling=type.filling) # return a DF
     test.if.multi.census <- fun.test.if.multi.census(df.lmer)
     cat("Ok data with Nobs",nrow(df.lmer),
         "multiple census", test.if.multi.census ,"\n")
        #= Run model
      if(test.if.multi.census){
        fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out)
Georges Kunstler's avatar
Georges Kunstler committed
88
        fun.ci.boot(df.lmer,formula=model$lmer.formula.tree.id,path.out,level=0.95,nsim=500)
Georges Kunstler's avatar
Georges Kunstler committed
89
90
      }else{
        fun.call.lmer.and.save(formula=model$lmer.formula,df.lmer,path.out)
Georges Kunstler's avatar
Georges Kunstler committed
91
        fun.ci.boot(df.lmer,formula=model$lmer.formula,path.out,level=0.95,nsim=500)
Georges Kunstler's avatar
Georges Kunstler committed
92
93
      }
}
Georges Kunstler's avatar
Georges Kunstler committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118


## new function to compute boot ci
 fun.ci.boot <- function(df.lmer,formula,path.out,level=0.95,nsim=500){
     require(boot)
     require(multicore)
        bb <- boot(data=df.lmer, statistic=boot.fun,R= nsim,formula=formula)
        bci <- lapply(seq_along(bb$t0), boot.out = bb, boot::boot.ci, 
            type = "perc", conf = level)
        citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
        a <- (1 - level)/2
        a <- c(a, 1 - a)
        pct <- paste('CI',round(a, 3),sep='.')
        dimnames(citab) <- list(names(bb[["t0"]]), pct)
        saveRDS(citab,file=file.path(path.out, "results.ci.rds"))
}

boot.fun <- function(data, indices, formula){
 df.lmer <- data[indices,] # select obs. in bootstrap sample
 res <- fun.call.lmer(formula=formula,df.lmer)
 fixef(res)
 }



Georges Kunstler's avatar
Georges Kunstler committed
119
120
121
#========================================================================

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


#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
                                  min.obs, sample.size, type.filling,  
Georges Kunstler's avatar
Georges Kunstler committed
131
                                  base.dir = "output/processed/",std=FALSE){
Georges Kunstler's avatar
Georges Kunstler committed
132
    ### load data
133
134
135
136
137
if(std) {   data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"),
                              stringsAsFactors = FALSE)}else{
    data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.csv"),
                              stringsAsFactors = FALSE)}

Georges Kunstler's avatar
Georges Kunstler committed
138
139
140
    fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}

141
142
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
                                         min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
Georges Kunstler's avatar
Georges Kunstler committed
143
144
145
146
147
148
149
150
151
152
153
154
## 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 
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
                    !is.na(data.tree[[BATOT]])),]
# 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])
155
156
# 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]]))
Georges Kunstler's avatar
Georges Kunstler committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
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=50){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= 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,
            species.id=species.id,
            tree.id=tree.id,
            plot.species.id=plot.species.id,
            sumTnTfBn.diff=sumTnTfBn.diff,
            sumTnTfBn.abs=sumTnTfBn.abs,
Georges Kunstler's avatar
Georges Kunstler committed
184
            Tf=data.tree[[tf]],
Georges Kunstler's avatar
Georges Kunstler committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
            sumTnBn=sumTnBn,
            sumTfBn=sumTfBn,
            sumBn=sumBn))
}

fun.get.the.variables.for.lmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= 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,
            species.id=species.id,
            plot.species.id=plot.species.id,
            sumTnTfBn.diff=sumTnTfBn.diff,
            sumTnTfBn.abs=sumTnTfBn.abs,
Georges Kunstler's avatar
Georges Kunstler committed
208
            Tf=data.tree[[tf]],
Georges Kunstler's avatar
Georges Kunstler committed
209
210
211
212
213
214
215
216
217
218
219
220
            sumTnBn=sumTnBn,
            sumTfBn=sumTfBn,
            sumBn=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
221
222
CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
Georges Kunstler's avatar
Georges Kunstler committed
223
224
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
225
226
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)
Georges Kunstler's avatar
Georges Kunstler committed
227
228
229
230
231
232
233
234
235
236
237
238
#= DATA LIST FOR JAGS
if (length(table(table(data.tree[["tree.id"]])))>1){
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
if (length(table(table(data.tree[["tree.id"]])))<2){
lmer.data <- fun.get.the.variables.for.lmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
    return(lmer.data)
}