lmer.output.figs.R 14.5 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
54
55
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
56
57
58
59
60
61
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))

62

Georges Kunstler's avatar
Georges Kunstler committed
63
## AIC weights
64
65
66
67
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
68
69
70
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
71

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


## 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


86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
## 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
111
112
113
114
115
116
117
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
118
library('pander')
119
fun.AIC.print.pandoc.table.trait <- function(DF.best.and.all.AIC,trait.select){
120
    DF.best.and.all.AIC$best.model <- factor(DF.best.and.all.AIC$best.model)
121
122
123
124
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))
125
126
print(pandoc.table(AIC.pandoc.table,caption=paste("Best models per data set for trait",trait.select)))}

127
128
129
130
131
132
133

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
134

Georges Kunstler's avatar
Georges Kunstler committed
135

136
137
138
#############################################
#############################################
### DO THE PLOT
Georges Kunstler's avatar
Georges Kunstler committed
139

140
141
142
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')
143
144
145

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
146
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
147
                            var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.06))
Georges Kunstler's avatar
Georges Kunstler committed
148
149
dev.off()

150
151
152
153
154
155
156
157
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
158

Georges Kunstler's avatar
Georges Kunstler committed
159
160
161
162
163
164
165
166
## 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
167
168


169

170
171
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))
172
par(mar=c(5.1,9.1,4.1,2.1))
173
boxplot(sumTfBn.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(sumTnBn.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
178
179
180
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))
181
par(mar=c(5.1,9.1,4.1,2.1))
182
183
184
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))
185
186
187
188
189


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))
190

191
boxplot(sumBn.VAR~model+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
192
193
194
195
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))
196
boxplot(sumTnBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
197
198
199
200
201
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')
202
203


204

Georges Kunstler's avatar
Georges Kunstler committed
205

Georges Kunstler's avatar
Georges Kunstler committed
206

207
208


Georges Kunstler's avatar
Georges Kunstler committed
209
### plot parameters
Georges Kunstler's avatar
Georges Kunstler committed
210

Georges Kunstler's avatar
Georges Kunstler committed
211

212

213
214
215
216
217
218
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)
219
220
221
fun.plot.panel.lmer.parameters.c(models=models,
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
222
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=10000)
223
224
dev.off()

225
226
227
228
229
230
231
232
233
234
235
236
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()

237

238
239
240
241
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'))
242

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

250
251
252
models <- c('lmer.LOGLIN.E.Tf')
names(models) <- c('Effect effect')
list.params <- list(c(Response='sumTnBn'))
253

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


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

266
267
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
268
269
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
270
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=100000)
Georges Kunstler's avatar
Georges Kunstler committed
271
272
273
dev.off()


274
275
276
277
278
279
280
281
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
282
283
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
284
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=10000)
Georges Kunstler's avatar
Georges Kunstler committed
285
dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
286

Georges Kunstler's avatar
Georges Kunstler committed
287

288

289
290
291
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')
list.params <- list('sumBn')
292

293
294
295
296
297
298
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()
299
300
301
302
303
304
305
306
307
308
309
310
311
312

## 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()



313