lmer.run.R 8.77 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
53
54


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

fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
   Start <- Sys.time()
55
   lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
56
57
58
   end <- Sys.time()
   print(end - Start)
   print(summary(lmer.output))
59
   saveRDS(lmer.output,file=file.path(path.out, "results.no.std.rds"))
Georges Kunstler's avatar
Georges Kunstler committed
60
61
62
63
64
65
66
67
68
69
70
}

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)
71
    print(path.out)
Georges Kunstler's avatar
Georges Kunstler committed
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
    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)
      }else{
        fun.call.lmer.and.save(formula=model$lmer.formula,df.lmer,path.out)
      }
}
#========================================================================

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


#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
                                  min.obs, sample.size, type.filling,  
                                  base.dir = "output/processed/"){
    ### load data
103
    data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"), stringsAsFactors = FALSE)
Georges Kunstler's avatar
Georges Kunstler committed
104
105
106
    fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}

107
108
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
109
110
111
112
113
114
115
116
117
118
119
120
## 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])
121
122
# 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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
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
150
            Tf=data.tree[[tf]],
Georges Kunstler's avatar
Georges Kunstler committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
            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
174
            Tf=data.tree[[tf]],
Georges Kunstler's avatar
Georges Kunstler committed
175
176
177
178
179
180
181
182
183
184
185
186
            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
187
188
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
189
190
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
191
192
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
193
194
195
196
197
198
199
200
201
202
203
204
#= 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)
}