glmer.run.all.R 9.34 KB
Newer Older
1
2
3
4
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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
###########################
###########################
### FUNCTION TO RUN GLMER ESTIMATION
library(lme4)



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

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


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


#=====================================================================
# Run glmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================

model.files.glmer.Tf.1 <-
    c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <-
    c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.AD.Tf.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.R")


model.files.glmer.Tf.3 <-
    c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.MAT.MAP.R")
model.files.glmer.Tf.4 <-
    c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.R",
      "R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.MAT.MAP.R")



fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
   Start <- Sys.time()
Georges Kunstler's avatar
Georges Kunstler committed
61
62
   glmer.output <- glmer(formula=formula,data=df.lmer,
                         family=binomial(link = 'cloglog'))
63
64
65
66
67
68
69
70
   end <- Sys.time()
   print(end - Start)
   print(summary(glmer.output))
   saveRDS(glmer.output,file=file.path(path.out, "glmer.results.global.rds"))
}

run.glmer <- function (model.file, trait,
                      min.obs=5, sample.size=NA,
Georges Kunstler's avatar
Georges Kunstler committed
71
                      type.filling, fname = 'data.all.no.std.csv') {
72
73
74
75
    require(lme4)
    source(model.file, local = TRUE)
    model <- load.model()
    #= Path for output
Georges Kunstler's avatar
Georges Kunstler committed
76
77
78
    if(fname == 'data.all.no.std.csv') dir <- 'all.no.log'
    if(fname == 'data.all.csv') dir <- 'all.log'
   path.out <- output.dir.glmer(dir, model$name, trait,
79
80
81
82
83
84
85
86
                                 type.filling=type.filling)
    print(path.out)
    dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
    cat("run glmer for model",model.file," for set",
         'all',"trait",
         trait,"\n")
      df.glmer <- load.and.prepare.data.for.glmer(trait,
                                           min.obs, sample.size,
Georges Kunstler's avatar
Georges Kunstler committed
87
88
                                         type.filling=type.filling,
                                                  fname = fname) # return a DF
89
90
91
92
93
94
95
     cat("Ok data with Nobs",nrow(df.glmer),"\n")
        #= Run model
        fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
}
#========================================================================


Georges Kunstler's avatar
Georges Kunstler committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
output.dir.glmer <- function (dir, model, trait, type.filling) {
  file.path("output/glmer", dir,trait,type.filling,model)
}



fun.load.data.all <- function(base.dir,fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
                       stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
                           intern = TRUE))

data.tree.tot <- read.csv(file = file.path(base.dir, fname),
                     stringsAsFactors = FALSE,
                     nrows = nrows,
                     colClasses = classes)
return(data.tree.tot)
115
116
117
118
119
120
121
122
}


#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait,
                                  min.obs, sample.size, type.filling,
Georges Kunstler's avatar
Georges Kunstler committed
123
                                  fname = 'data.all.no.std.csv',
124
125
                                  base.dir = "output/processed/"){
    ### load data
Georges Kunstler's avatar
Georges Kunstler committed
126
    data.tree.tot <- fun.load.data.all(base.dir,fname)
127
128
129
130

    fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
}

Georges Kunstler's avatar
Georges Kunstler committed
131
132
133
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,
                                         perc.neigh.sp,
                                         perc.neigh.gs,
134
                                         BATOT,min.obs,
Georges Kunstler's avatar
Georges Kunstler committed
135
136
                                         min.perc.neigh.sp = 0.90,
                                         min.perc.neigh.gs = 0.95){
137
138
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
Georges Kunstler's avatar
Georges Kunstler committed
139
140
                                     (!is.na(data.tree[["D"]])) &
                                     (!is.na(data.tree[['year']])))
141
142
143
144
145
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
                    !is.na(data.tree[[BATOT]])),]
Georges Kunstler's avatar
Georges Kunstler committed
146
147
148
149
150
151
152
153
# select  obs abov min perc.neigh species
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
                    min.perc.neigh.sp &
                    !is.na(data.tree[[perc.neigh.sp]]))
# select  obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
                    min.perc.neigh.gs & !is.na(data.tree[[perc.neigh.gs]]))

154
155
156
157
158
159
160
161
162
163
164
165
# 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])
return(data.tree)
}

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

### get variables for lmer
Georges Kunstler's avatar
Georges Kunstler committed
166
fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,
167
168
169
                                                    abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
Georges Kunstler's avatar
Georges Kunstler committed
170
logyear <- log(data.tree[["year"]])
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
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']])
#= 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(dead=dead,
            logD=logD,
Georges Kunstler's avatar
Georges Kunstler committed
186
            logyear=logyear,
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
            MAT = MAT,
            MAP = MAP,
            species.id=species.id,
            tree.id=tree.id,
            plot.id=plot.id,
            set.id = set.id,
            ecocode.id = ecocode.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)))
}


#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
fun.data.for.glmer <-  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
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
218
219
220
221
222
223
224
perc.neigh.sp <- paste(trait, "perc", 'species', sep = ".")
perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
data.tree <- fun.select.data.for.analysis(data.tree,
                                          abs.CWM.tntf,
                                          perc.neigh.sp,
                                          perc.neigh.gs,
                                          BATOT,min.obs)
225
226
227
228
229
230
231
#= DATA LIST FOR JAGS
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
    return(glmer.data)
}