glmer.nolog.output.figs.R 5.86 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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#!/usr/bin/env Rscript
rm(list = ls())
source("R/analysis/lmer.output-fun.R")
source("R/utils/plot.R")

## load results
DF.results <- readRDS('output/glmer.global.std.DF.rds')
table(DF.results$conv)
DF.results$id <- paste(DF.results$set, DF.results$ecocode, sep = ".")

## climate
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 = ".")
mean.MAT <- tapply(site.clim.all$MAT, INDEX = site.clim.all$id,
                   FUN = mean, na.rm = TRUE)
mean.MAP <- tapply(site.clim.all$MAP, INDEX = site.clim.all$id,
                   FUN = mean, na.rm = TRUE)
unique.set <- tapply(site.clim.all$set, INDEX = site.clim.all$id,
                   FUN = unique, na.rm = TRUE)

data.clim.ecocode <- data.frame(id = names(mean.MAT),
                                MAT = mean.MAT,
                                MAP = mean.MAP,
                                set = unique.set, stringsAsFactors = FALSE)

biomes.vec <- fun.overly.plot.on.biomes(MAP = data.clim.ecocode$MAP/10,
                                        MAT = data.clim.ecocode$MAT,
                          names.vec = 1:nrow(data.clim.ecocode))
data.clim.ecocode$biomes <- as.character(unlist(biomes.vec))


### merge climate and lmer results
DF.results <- merge(DF.results,
                    data.clim.ecocode[, c('id','MAT','MAP','biomes')] ,
                    by = "id")
DF.results$id2 <- paste(DF.results$id, DF.results$trait, DF.results$filling)

## add value for NVS
DF.results$biomes[DF.results$ecocode == 'MixedCool'] <-  'temperate_rainforest'


## ADD percentage traits
site.perc.all <- read.csv('output/processed/all.sites.perc.traits.csv', 
                          stringsAsFactors = FALSE)
site.perc.all$id <- paste(site.perc.all$set, site.perc.all$ecocode, sep = ".")
##
DF.results <- merge(DF.results, 
                    subset(site.perc.all, 
                           select = c('id', 'perc.gymno',
                               'perc.ev', 'BATOT.m')), 
                    by = "id")

### Analysis of the results

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

## best model
DF.best.and.all.AIC <- do.call('rbind', lapply(unique(DF.results$id2), 
                                              FUN = fun.AIC, DF.results))
table(DF.best.and.all.AIC$best.model, DF.best.and.all.AIC$trait)

# global AIC
global.AIC.list <- fun.global.aic(DF.results)


### plots

models <- c('glmer.LOGLIN.ER.Tf', 'glmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response effect', 'Effect/response response')
list.params <- list(c(Response = 'sumTnBn'), 
                    c(Effect = 'sumTfBn'))

## pdf('figs/parameters.glmer.MAP.ER.all.nolog.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 = 'MAP', 
                               list.params = list.params, 
                                 small.bar = 0.2, threshold.delta.AIC = 10000000, 
                                 ylim = c(-0.2, 0.2), ylim.same = TRUE, 
                                 log = 'x', lwd = 1.5)
## dev.off()

models <- c('glmer.LOGLIN.AD.Tf')
names(models) <- c('Absolut distance effect')
list.params <- list(c(Response = 'sumTnTfBn.abs'))

## pdf('figs/parameters.glmer.MAP.ER.all.nolog.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 = 'MAP', 
                               list.params = list.params, 
                                 small.bar = 0.2, threshold.delta.AIC = 10000000, 
                                 ylim = c(-0.005, 0.2), ylim.same = FALSE, 
                                 log = 'x', lwd = 1.5)
## dev.off()
 

models <- c('glmer.LOGLIN.ER.Tf')
names(models) <- c('ER Tf')
list.params <- list(c(Response = 'Tf'))
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.2, threshold.delta.AIC = 10000000, 
                                 ylim = c(-1, 1), ylim.same = FALSE, log = 'x')


models <- c('glmer.LOGLIN.ER.Tf')
names(models) <- c('ER Tf')
list.params <- list(c(Response = 'logD'))
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.2, threshold.delta.AIC = 10000000, 
                                 ylim = c(-1, 1), ylim.same = FALSE, log = 'x')