lmer.output.figs.R 14.6 KB
Newer Older
1
#!/usr/bin/env Rscript
Georges Kunstler's avatar
Georges Kunstler committed
2
rm(list=ls())
3
source("R/analysis/lmer.output-fun.R")
Georges Kunstler's avatar
Georges Kunstler committed
4
source("R/utils/plot.R")
5

Georges Kunstler's avatar
Georges Kunstler committed
6
## load results
7
list.lmer  <- readRDS("output/lmer.list.out.rds")
8
9
## get names of all parameters
names.param <- unique(unlist(lapply(list.lmer,function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
10
11
12
names.model <- unique(unlist(lapply(list.lmer,function(list.res) (list.res$files.details[5]))))
names.set <- unique(unlist(lapply(list.lmer,function(list.res) (list.res$files.details[1]))))
# turn to a data frame
13
DF.results <- do.call("rbind",lapply(list.lmer,fun.format.in.data.frame,names.param=names.param))
14
# add id
15
16
17
18
19
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")

## load climatic data
site.clim.all <- read.csv( file.path("output/processed", "all.sites.clim.csv"), stringsAsFactors = FALSE)
site.clim.all$id <- paste(site.clim.all$set,site.clim.all$ecocode,sep=".")
Georges Kunstler's avatar
Georges Kunstler committed
20
21
22
## par(mfrow=c(1,2),mar=c(5, 9, 4, 1) +0.1)
## boxplot(site.clim.all$clim.all.MAT~site.clim.all$id,horizontal=TRUE,las=2,cex.axis=0.8)
## boxplot(site.clim.all$clim.all.MAP~site.clim.all$id,horizontal=TRUE,las=2,cex.axis=0.8)
23
24
25
26
27
28
mean.MAT <- tapply(site.clim.all$clim.all.MAT,INDEX=site.clim.all$id,FUN=mean,na.rm=TRUE)
mean.MAP <- tapply(site.clim.all$clim.all.MAP,INDEX=site.clim.all$id,FUN=mean,na.rm=TRUE)
data.clim.ecocode <- data.frame(id=names(mean.MAT),MAT=mean.MAT,MAP=mean.MAP)
### merge climate and lmer results
DF.results <- merge(DF.results,data.clim.ecocode,by="id")

29
DF.results$id2 <- paste(DF.results$id,DF.results$trait,DF.results$filling)
Georges Kunstler's avatar
Georges Kunstler committed
30

Georges Kunstler's avatar
Georges Kunstler committed
31
### save
32
saveRDS(DF.results,file='output/lmer.DF.results.merged.rds')
Georges Kunstler's avatar
Georges Kunstler committed
33

34
DF.results <- readRDS(file='output/lmer.DF.results.merged.rds')
Georges Kunstler's avatar
Georges Kunstler committed
35

36
### Analysis of the results
Georges Kunstler's avatar
Georges Kunstler committed
37
38
39
40
41

DF.R2m.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"R2m"))
DF.R2c.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"R2c"))
DF.AIC.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"AIC"))
DF.delta.AIC <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.delta.AIC,DF.results))
Georges Kunstler's avatar
Georges Kunstler committed
42

43
44
DF.var.fixed <- fun.ratio.var.fixed.effect(DF.results)
DF.results <- cbind(DF.results,DF.R2m.diff,DF.R2c.diff,DF.AIC.diff,DF.delta.AIC,DF.var.fixed)
Georges Kunstler's avatar
Georges Kunstler committed
45
46
47



48
### report best model based on AIC
Georges Kunstler's avatar
Georges Kunstler committed
49
50
51
DF.best.and.all.AIC <- do.call('rbind',lapply(unique(DF.results$id2),FUN=fun.AIC,DF.results))
DF.best.and.all.AICc <- do.call('rbind',lapply(unique(DF.results$id2),FUN=fun.AICc,DF.results))

52
53
table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model)

54
55
56
57
table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model,
      DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$trait,
      DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set)

Georges Kunstler's avatar
Georges Kunstler committed
58
59
60
61
62
63
t(table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model,
      DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set))/(apply(table(DF.best.and.all.AIC[
                              DF.best.and.all.AIC$filling=='species',]$best.model,
      DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set),
      MARGIN=2,sum))

64

Georges Kunstler's avatar
Georges Kunstler committed
65
## AIC weights
66
67
68
69
AIC.weights <- do.call('rbind',
                       lapply(1:nrow(DF.best.and.all.AICc),
                              FUN=function(i,DF) exp((min(DF[i,])-DF[i,])/2)/sum(exp((min(DF[i,])-DF[i,])/2)),
                              DF.best.and.all.AIC[,9:15]))
Georges Kunstler's avatar
Georges Kunstler committed
70
71
72
DF.AIC.weights <- data.frame(DF.best.and.all.AICc[,1],AIC.weights)
names(DF.AIC.weights) <- c('id2',paste('AIC.weight',names(DF.AIC.weights)[-1],sep='.'))
DF.best.and.all.AIC <- merge(DF.best.and.all.AIC,DF.AIC.weights,by='id2')
Georges Kunstler's avatar
Georges Kunstler committed
73

Georges Kunstler's avatar
Georges Kunstler committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87


## assume ER best model if either E or R best model
f <- function(i,DF){
d.AIC <- as.vector(DF$delta.AIC[DF$id2==i]    )
best <- DF$model[DF$id2==i][DF$delta.AIC[DF$id2==i]==0]
if(best %in% c('lmer.LOGLIN.E','lmer.LOGLIN.R'))  {
    d.AIC[ DF$model[DF$id2==i]=='lmer.LOGLIN.ER'] <- 0}
return(d.AIC)
}
d.AIC <- do.call('c',lapply(unique(DF.results$id2),FUN=f,DF.results))
DF.results$delta.AIC <- d.AIC


88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
## compute a global AIC
fun.global.aic <- function(DF.results){
DF.results <- DF.results[DF.results$nobs>1000,] # select set ecocode with more than 1000 obs
DF.results <- DF.results[DF.results$id2 %in%  names(table(DF.results$id2))[table(DF.results$id2)==7],] # select set ecocode with 7 model tested
## species
DF.results.sp <- DF.results[DF.results$filling=='species',] # select set ecocode with more than 1000 obs
test.same.n.model.ecoregion <- apply(table(DF.results.sp$trait,DF.results.sp$model),MARGIN=1,function(x) all(x==x[1]))
if(!all(test.same.n.model.ecoregion))
    stop(paste('error not the same number of ecoregion for traits',
               names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
list.sp <- list(AIC.tot=tapply(DF.results.sp$AIC,INDEX=list(DF.results.sp$trait,DF.results.sp$model),FUN=sum),
N.ecoregion=tapply(DF.results.sp$AIC,INDEX=list(DF.results.sp$trait,DF.results.sp$model),FUN=length))
## genus
DF.results.ge <- DF.results[DF.results$filling=='genus',] # select set ecocode with more than 1000 obs
test.same.n.model.ecoregion <- apply(table(DF.results.ge$trait,DF.results.ge$model),MARGIN=1,function(x) all(x==x[1]))
if(!all(test.same.n.model.ecoregion))
    stop(paste('error not the same number of ecoregion for traits',
               names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
list.ge <- list(AIC.tot=tapply(DF.results.ge$AIC,INDEX=list(DF.results.ge$trait,DF.results.ge$model),FUN=sum),
N.ecoregion=tapply(DF.results.ge$AIC,INDEX=list(DF.results.ge$trait,DF.results.ge$model),FUN=length))
return(list(sp=list.sp,ge=list.ge))
}    
global.AIC.list <- fun.global.aic(DF.results)

#### compute percentage of variance explained by var
113
114
115
116
117
118
119
DF.results$abs.perc.var <- DF.results$sumTnTfBn.abs.VAR/DF.results$sumBn.VAR
DF.results$R.perc.var <- DF.results$sumTfBn.VAR/DF.results$sumBn.VAR
DF.results$E.perc.var <- DF.results$sumTnBn.VAR/DF.results$sumBn.VAR
DF.results$ER.perc.var <- DF.results$effect.response.var/DF.results$sumBn.VAR


## print AIC table in markdown format
Georges Kunstler's avatar
Georges Kunstler committed
120
library('pander')
121
fun.AIC.print.pandoc.table.trait <- function(DF.best.and.all.AIC,trait.select){
122
    DF.best.and.all.AIC$best.model <- factor(DF.best.and.all.AIC$best.model)
123
124
125
126
aa <- as.data.frame.matrix( t(with(subset(DF.best.and.all.AIC,subset=trait %in% trait.select),
                                 table(best.model,set))))
AIC.pandoc.table <- data.frame(set=rownames(aa),aa[,c(4,6,1,5,2,3)],row.names=NULL)
names(AIC.pandoc.table) <- gsub("lmer.LOGLIN.","", names(AIC.pandoc.table))
127
128
print(pandoc.table(AIC.pandoc.table,caption=paste("Best models per data set for trait",trait.select)))}

129
130
131
132
133
134
135

fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="SLA")
fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="Leaf.N")
fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="Wood.density")
fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="Max.height")

fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select=c("SLA","Leaf.N","Wood.density","Max.height"))
Georges Kunstler's avatar
Georges Kunstler committed
136

Georges Kunstler's avatar
Georges Kunstler committed
137

138
139
140
#############################################
#############################################
### DO THE PLOT
Georges Kunstler's avatar
Georges Kunstler committed
141

142
143
144
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
var.y.l <- list('R2m.simplecomp','R2m.simplecomp')
145
146
147

pdf('figs/R2m.MAP.pdf',width=12,height=8)
fun.plot.panel.lmer.res.x.y(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
148
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
149
                            var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.06))
Georges Kunstler's avatar
Georges Kunstler committed
150
151
dev.off()

152
153
154
155
156
157
158
159
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
var.y.l <- list('effect.response.var','sumTnTfBn.abs.VAR')
pdf('figs/perc.var.MAP.pdf',width=12,height=8)
fun.plot.panel.lmer.res.x.y(models=models,
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
                            var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.17))
dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
160

Georges Kunstler's avatar
Georges Kunstler committed
161
162
163
164
165
166
167
168
## models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
## names(models) <- c('Effect/response','Absolute distance')
## var.y.l <- list('ER.perc.var','abs.perc.var')
## pdf('figs/perc.var.relative.BATOT.MAP.pdf',width=12,height=8)
## fun.plot.panel.lmer.res.x.y(models=models,
##                             traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
##                             var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,100))
## dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
169
170


171

172
173
par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(2,2))
boxplot(effect.response.var~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
174
par(mar=c(5.1,9.1,4.1,2.1))
175
boxplot(sumTfBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
176
par(mar=c(5.1,9.1,4.1,2.1))
177
boxplot(sumTnBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
178
par(mar=c(5.1,9.1,4.1,2.1))
179
180
181
182
boxplot(sumTnTfBn.abs.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))

par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(2,2))
boxplot(sumBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
183
par(mar=c(5.1,9.1,4.1,2.1))
184
185
186
boxplot(Tf.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
par(mar=c(5.1,9.1,4.1,2.1))
boxplot(logD.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
187
188
189
190
191


par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(1,2))
boxplot(ER.perc.var~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,100))
boxplot(abs.perc.var~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,100))
192

193
boxplot(sumBn.VAR~model+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
194
195
196
197
boxplot(R2m~trait+model+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
boxplot(R2m~set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)

par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(1,3))
198
boxplot(sumTnBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
199
200
201
202
203
abline(v=0,col='red')
boxplot(sumTfBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
abline(v=0,col='red')
boxplot(sumBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
abline(v=0,col='red')
204
205


206

Georges Kunstler's avatar
Georges Kunstler committed
207

Georges Kunstler's avatar
Georges Kunstler committed
208

209
210


Georges Kunstler's avatar
Georges Kunstler committed
211
### plot parameters
Georges Kunstler's avatar
Georges Kunstler committed
212

Georges Kunstler's avatar
Georges Kunstler committed
213

214

215
216
217
218
219
220
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response effect','Effect/response response')
list.params <- list(c(Response='sumTnBn'),
                    c(Effect='sumTfBn'))

pdf('figs/parameters.MAP.ER.all.pdf',width=9,height=7)
221
222
223
fun.plot.panel.lmer.parameters.c(models=models,
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
224
                               list.params=list.params,small.bar=0.02,
225
                                 threshold.delta.AIC=10000,ylim=c(-0.2,0.2),ylim.same=TRUE)
226
227
dev.off()

228
229
230
231
232
233
234
235
236
237
238
239
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response effect','Effect/response response')
list.params <- list(c(Response='sumTnBn'),
                    c(Effect='sumTfBn'))

pdf('figs/parameters.MAT.ER.all.pdf',width=9,height=7)
fun.plot.panel.lmer.parameters.c(models=models,
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAT',
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=10000)
dev.off()

240

241
242
243
244
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response effect','Effect/response response')
list.params <- list(c(Response='sumTnBn'),
                    c(Effect='sumTfBn'))
245

246
pdf('figs/parameters.MAP.ER.best.aic.pdf',width=9,height=7)
247
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
248
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
249
                               DF.results,var.x='MAP',
250
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=1)
251
252
dev.off()

253
254
255
models <- c('lmer.LOGLIN.E.Tf')
names(models) <- c('Effect effect')
list.params <- list(c(Response='sumTnBn'))
256

257
pdf('figs/parameters.MAP.E.pdf',width=9,height=5)
258
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
259
260
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
261
                               list.params=list.params,small.bar=0.02,threshold.delta.AIC=1000)
Georges Kunstler's avatar
Georges Kunstler committed
262
263
264
dev.off()


265
266
267
models <- c('lmer.LOGLIN.R.Tf')
names(models) <- c('Response response')
list.params <- list(c(Response='sumTfBn'))
Georges Kunstler's avatar
Georges Kunstler committed
268

269
270
pdf('figs/parameters.MAP.R.pdf',width=9,height=5)
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
271
272
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
273
                               list.params=list.params,small.bar=0.02,threshold.delta.AIC=100000)
Georges Kunstler's avatar
Georges Kunstler committed
274
275
276
dev.off()


277
278
279
280
281
282
283
284
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response effect','Effect/response response','Absolute distance')
list.params <- list(c(Response='sumTnBn'),
                    c(Effect='sumTfBn'),
                    c('sumTnTfBn.abs'))

pdf('figs/parameters.MAP.2models.pdf',width=12,height=10)
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
285
286
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
287
                               list.params=list.params,small.bar=0.02,threshold.delta.AIC=10000)
Georges Kunstler's avatar
Georges Kunstler committed
288
dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
289

Georges Kunstler's avatar
Georges Kunstler committed
290

291

292
293
294
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')
list.params <- list('sumBn')
295

296
297
298
299
300
301
pdf('figs/parameters.MAP.sumBn.ER.pdf',width=8,height=3)
fun.plot.panel.lmer.parameters.c(models=models,
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
                               list.params=list.params,small.bar=0.02,ylim=c(-0.4,0.13),threshold.delta.AIC=10000)
dev.off()
302
303
304
305
306
307
308
309
310
311
312
313
314
315

## models <- c('lmer.LOGLIN.AD')
## names(models) <- c('Absolute distance')
## list.params <- list(c('sumTnTfBn.abs'))

## pdf('figs/parameters.boxplot.AD.pdf',width=8,height=3)
## fun.plot.panel.lmer.parameters.boxplot(models=models,
##                                traits = c('Wood.density','SLA','Leaf.N','Max.height'),
##                                DF.results,var.x='MAP',
##                                list.params=list.params,small.bar=0.02,ylim=c(-0.1,0.13))
## dev.off()



316