lmer.output.figs.R 11.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
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")

17

18
19
20
## 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
21
22
23
## 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)
24
25
26
27
28
29
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")

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

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

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

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

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
43

44
45
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
46
47
48



49
### report best model based on AIC
Georges Kunstler's avatar
Georges Kunstler committed
50
51
52
53
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))

## AIC weights
54
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:14]))
Georges Kunstler's avatar
Georges Kunstler committed
55
56
57
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
58
59


60
61
62
63
64
65
66
67
68
69
70
71
#### compute percentage of vqariance explained by var
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
72
library('pander')
73
fun.AIC.print.pandoc.table.trait <- function(DF.best.and.all.AIC,trait.select){
74
    DF.best.and.all.AIC$best.model <- factor(DF.best.and.all.AIC$best.model)
75
76
77
78
79
80
81
82
83
84
85
86
87
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))
print(pandoc.table(AIC.pandoc.table,caption=paste("Best models per data set for trait",trait.select)))
}

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
88

Georges Kunstler's avatar
Georges Kunstler committed
89

90
91
92
#############################################
#############################################
### DO THE PLOT
Georges Kunstler's avatar
Georges Kunstler committed
93

94
95
96
97
98
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')

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
99
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
100
                            var.x='MAP',var.y='R2m.simplecomp',ylim=c(-0.015,0.06))
Georges Kunstler's avatar
Georges Kunstler committed
101
102
dev.off()

103
104
105
106
107
108
109
110
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
111

112
113
114
115
116
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.MAP.pdf',width=12,height=8)
fun.plot.panel.lmer.res.x.y(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
117
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
118
                            var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,100))
Georges Kunstler's avatar
Georges Kunstler committed
119
120
121
dev.off()


122
123
124
models <- c('lmer.LOGLIN.AD.Tf')
names(models) <- c()

Georges Kunstler's avatar
Georges Kunstler committed
125
126
fun.plot.panel.lmer.res.x.y(models=models,
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
127
                            var.x='MAP',var.y='sumTnTfBn.abs.VAR',ylim=c(-0.015,0.2))
Georges Kunstler's avatar
Georges Kunstler committed
128
129


130
131
132
133

models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')

Georges Kunstler's avatar
Georges Kunstler committed
134
135
fun.plot.panel.lmer.res.x.y(models=models,
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
136
                            var.x='MAP',var.y='ER.perc.var',ylim=c(-0.015,100))
Georges Kunstler's avatar
Georges Kunstler committed
137

138
139
models <- c('lmer.LOGLIN.AD.Tf')
names(models) <- c('AB')
140
141
142

fun.plot.panel.lmer.res.x.y(models=models,
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
                            var.x='MAP',var.y='abs.perc.var',ylim=c(-0.015,100))

par(mar=c(5.1,9.1,4.1,2.1))
boxplot(effect.response.var~trait+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))
boxplot(sumTfBn.VAR~trait,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))
boxplot(sumTnBn.VAR~trait,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))
boxplot(sumTnTfBn.abs.VAR~trait+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,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))
158

159
160
161
162
163
164
165
boxplot(sumBn.VAR~model+trait,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))
boxplot(sumTnBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)


models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
166

167
168
pdf('figs/R2.MAP.two.pdf',width=10,height=7)
fun.plot.panel.lmer.res.x.y(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
169
                            traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
170
                            var.x='MAP',var.y='R2m.simplecomp',ylim=c(-0.015,0.08))
Georges Kunstler's avatar
Georges Kunstler committed
171
172
dev.off()

Georges Kunstler's avatar
Georges Kunstler committed
173

174
175


Georges Kunstler's avatar
Georges Kunstler committed
176
### plot parameters
Georges Kunstler's avatar
Georges Kunstler committed
177

Georges Kunstler's avatar
Georges Kunstler committed
178

179

180
181
182
183
184
185
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)
186
187
188
fun.plot.panel.lmer.parameters.c(models=models,
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
189
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=10000)
190
191
192
dev.off()


193
194
195
196
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'))
197

198
pdf('figs/parameters.MAP.ER.best.aic.pdf',width=9,height=7)
199
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
200
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
201
                               DF.results,var.x='MAP',
202
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=1)
203
204
dev.off()

205
206
207
models <- c('lmer.LOGLIN.E.Tf')
names(models) <- c('Effect effect')
list.params <- list(c(Response='sumTnBn'))
208

209
pdf('figs/parameters.MAP.E.pdf',width=9,height=5)
210
fun.plot.panel.lmer.parameters.c(models=models,
Georges Kunstler's avatar
Georges Kunstler committed
211
212
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
213
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=1000)
Georges Kunstler's avatar
Georges Kunstler committed
214
215
216
dev.off()


217
218
219
models <- c('lmer.LOGLIN.R.Tf')
names(models) <- c('Response response')
list.params <- list(c(Response='sumTfBn'))
Georges Kunstler's avatar
Georges Kunstler committed
220

221
222
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
223
224
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
225
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=100000)
Georges Kunstler's avatar
Georges Kunstler committed
226
227
228
dev.off()


229
230
231
232
233
234
235
236
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
237
238
                               traits = c('Wood.density','SLA','Leaf.N','Max.height'),
                               DF.results,var.x='MAP',
239
                               list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=10000)
Georges Kunstler's avatar
Georges Kunstler committed
240
dev.off()
Georges Kunstler's avatar
Georges Kunstler committed
241

Georges Kunstler's avatar
Georges Kunstler committed
242

243

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
## models <- c('lmer.LOGLIN.ER')
## names(models) <- c('Effect/response')
## list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'))

## pdf('figs/parameters.boxplot.ER.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),col=c('red','black'))
## dev.off()

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



268