lmer.output-fun.R 28.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
#### function to analyse lmer output


library(lme4)


read.lmer.output <- function(file.name){
  tryCatch(readRDS(file.name), error=function(cond)return(NULL))   # Choose a return value in case of error
}



summarise.lmer.output <- function(x){
14
15
16
17
18
 list( nobs = nobs(x),
       R2m =Rsquared.glmm.lmer(x)$Marginal,
       R2c =Rsquared.glmm.lmer(x)$Conditional,
       AIC = AIC(x),
       deviance = deviance(x),
19
       conv=x@optinfo$conv,
20
21
22
23
24
       effect.response.var=variance.fixed.glmm.lmer.effect.and.response(x),
       fixed.coeff.E=fixef(x),
       fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
       fixed.var=variance.fixed.glmm.lmer(x))
}
25

Georges Kunstler's avatar
Georges Kunstler committed
26
27
28
29
30
31
32
33
34
35
36
37
38
summarise.lmer.all.output <- function(x){
 list( nobs = nobs(x),
       R2m =Rsquared.glmm.lmer(x)$Marginal,
       R2c =Rsquared.glmm.lmer(x)$Conditional,
       AIC = AIC(x),
       deviance = deviance(x),
       conv=x@optinfo$conv,
       effect.response.var=variance.fixed.glmm.lmer.effect.and.response(x),
       fixed.coeff.E=fixef(x),
       fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
       fixed.var=variance.fixed.glmm.lmer(x),
       set.BLUP=coef(x)$set.id)
}
39
40
41
42
43
44
45
46
47
48
49
50
51

summarise.lmer.output.list <- function(f ){
    output.lmer <- read.lmer.output(f)
    if(!is.null(output.lmer)){
	res <- list(files.details=files.details(f),
		   lmer.summary=summarise.lmer.output( output.lmer))
    }else{
        res <- NULL
    }
    return(res)
}


52
53
54
55
summarise.lmer.output.all.list <- function(f ){
    output.lmer <- read.lmer.output(f)
    if(!is.null(output.lmer)){
	res <- list(files.details=files.details.all(f),
Georges Kunstler's avatar
Georges Kunstler committed
56
		   lmer.summary=summarise.lmer.all.output( output.lmer))
57
58
59
60
61
62
63
    }else{
        res <- NULL
    }
    return(res)
}


64
65

files.details <- function(x){
Georges Kunstler's avatar
Georges Kunstler committed
66
67
68
69
	s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]),
                        stringsAsFactors= FALSE)
	names(s)  <- c("d1", "d2", "set", "ecocode", "trait",
                       "filling", "model", "file" )
70
71
72
73
	s[-(1:2)]
}


74
files.details.all <- function(x){
Georges Kunstler's avatar
Georges Kunstler committed
75
76
77
78
	s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]),
                        stringsAsFactors= FALSE)
	names(s)  <- c("d1", "d2", "set", "trait",
                       "filling", "model", "file" )
79
80
81
	s[-(1:2)]
}

82

83
#### R squared functions
84

Georges Kunstler's avatar
Georges Kunstler committed
85
86
87
88
89
90
# Function rsquared.glmm requires models to be input as a list
# (can include fixed-
# effects only models,but not a good idea to mix models of
# class "mer" with models
# of class "lme") FROM
#  http://jslefche.wordpress.com/2013/03/13/r2-for-linear-mixed-effects-models/
91
92

Rsquared.glmm.lmer <- function(i){
93
94
95
# Get variance of fixed effects by multiplying coefficients by design matrix
      VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
      # Get variance of random effects by extracting variance components
Georges Kunstler's avatar
Georges Kunstler committed
96
      VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1])))
97
      # Get residual variance
Georges Kunstler's avatar
Georges Kunstler committed
98
      VarResid <- attr(VarCorr(i), "sc")^2
99
100
101
102
103
      # Calculate marginal R-squared (fixed effects/total variance)
      Rm <- VarF/(VarF+VarRand+VarResid)
      # Calculate conditional R-squared (fixed effects+random effects/total variance)
      Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid)
      # Bind R^2s into a matrix and return with AIC values
Georges Kunstler's avatar
Georges Kunstler committed
104
      Rsquared.mat <- data.frame(Class=class(i), Family="Gaussian", Marginal=Rm,
Georges Kunstler's avatar
Georges Kunstler committed
105
                              Conditional=Rc, AIC=AIC(i))
106
107
108
109
110
      return(Rsquared.mat)
}


variance.fixed.glmm.lmer <- function(i){
111

112
# Get variance of for each fixed effects by multiplying coefficients by design matrix
Georges Kunstler's avatar
Georges Kunstler committed
113
var.vec <- apply(fixef(i) * t(i@pp$X), MARGIN=1, var)
114
115
116
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
Georges Kunstler's avatar
Georges Kunstler committed
117
VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1])))
118
# Get residual variance
Georges Kunstler's avatar
Georges Kunstler committed
119
VarResid <- attr(VarCorr(i), "sc")^2
120
var.vec <- var.vec/(VarF+VarRand+VarResid)
Georges Kunstler's avatar
Georges Kunstler committed
121
names(var.vec) <- paste(names(var.vec), "VAR", sep=".")
122
123
124
125
return(var.vec)
}

variance.fixed.glmm.lmer.effect.and.response <- function(i){
Georges Kunstler's avatar
Georges Kunstler committed
126
if(sum(c("sumTfBn", "sumTnBn") %in% names(fixef(i)))==2){
127
# Get variance of for each fixed effects by multiplying coefficients by design matrix
Georges Kunstler's avatar
Georges Kunstler committed
128
129
var.vec <- var(as.vector(fixef(i)[c("sumTfBn", "sumTnBn")] %*%
                         t(i@pp$X[, c("sumTfBn", "sumTnBn")])))
130
131
132
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
Georges Kunstler's avatar
Georges Kunstler committed
133
VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1])))
134
# Get residual variance
Georges Kunstler's avatar
Georges Kunstler committed
135
VarResid <- attr(VarCorr(i), "sc")^2
136
137
138
var.vec <- var.vec/(VarF+VarRand+VarResid)
}else{
var.vec <- NA
Georges Kunstler's avatar
Georges Kunstler committed
139
}
Georges Kunstler's avatar
Georges Kunstler committed
140
names(var.vec) <- paste("effect.response", "VAR", sep=".")
141
142
return(var.vec)
}
143
144


Georges Kunstler's avatar
Georges Kunstler committed
145
146
147


## function to turn lmer output from list to DF
Georges Kunstler's avatar
Georges Kunstler committed
148
149
150
151
152
fun.format.in.data.frame <- function(list.res, names.param){
dat.t <- data.frame(t(rep(NA, 3*length(names.param))))
names(dat.t) <-  c(names.param, paste(names.param, "Std.Error", sep=".")
                   , paste(names.param, "VAR", sep="."))
dat.t[, match(names(list.res$lmer.summary$fixed.coeff.E), names(dat.t))] <-
Georges Kunstler's avatar
Georges Kunstler committed
153
    list.res$lmer.summary$fixed.coeff.E
Georges Kunstler's avatar
Georges Kunstler committed
154
dat.t[, length(names.param)+match(names(list.res$lmer.summary$fixed.coeff.E),
Georges Kunstler's avatar
Georges Kunstler committed
155
                                 names(dat.t))] <-
Georges Kunstler's avatar
Georges Kunstler committed
156
    list.res$lmer.summary$fixed.coeff.Std.Error
Georges Kunstler's avatar
Georges Kunstler committed
157
dat.t[, match(names(list.res$lmer.summary$fixed.var), names(dat.t))] <-
158
    list.res$lmer.summary$fixed.var
Georges Kunstler's avatar
Georges Kunstler committed
159
res <- data.frame(list.res$files.details, list.res$lmer.summary[1:7], dat.t)
Georges Kunstler's avatar
Georges Kunstler committed
160
161
return(res)
}
162
163
164
165
166
167
168


#################################################################
#################################################################
##### FUNCTION TO analyse the results AIC effect size

## function to compute delat R2
Georges Kunstler's avatar
Georges Kunstler committed
169
fun.compute.criteria.diff <- function(i, DF.results, criteria.selected){
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
select.simple.compet <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i] &
                    DF.results$model=='lmer.LOGLIN.simplecomp.Tf'
select.no.compet <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i] &
                    DF.results$model=='lmer.LOGLIN.nocomp.Tf'
if(sum(select.simple.compet)==1){
diff.criteria.simple.compet <- DF.results[[criteria.selected]][i] - DF.results[[criteria.selected]][
                                                                                select.simple.compet]
}else{
diff.criteria.simple.compet <- NA
}
if(sum(select.no.compet)==1){
diff.criteria.no.compet <- DF.results[[criteria.selected]][i] - DF.results[[criteria.selected]][
                                                                              select.no.compet]
}else{
diff.criteria.no.compet <- NA
}

Georges Kunstler's avatar
Georges Kunstler committed
191
192
df.res <- data.frame(diff.criteria.simple.compet, diff.criteria.no.compet)
names(df.res) <- paste(criteria.selected, c('simplecomp', 'nocomp'), sep=".")
193
194
195
196
197
198

return(df.res)
}



Georges Kunstler's avatar
Georges Kunstler committed
199
fun.compute.delta.AIC <- function(i, DF.results){
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
select.model.trait.fill <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i]
if(sum(select.model.trait.fill)>0){
delta.AIC <- DF.results[['AIC']][i] - min(DF.results[['AIC']][select.model.trait.fill])
if (sum( DF.results$nobs[select.model.trait.fill]!=DF.results$nobs[i])>0)
    stop('no same number of observation')
}else{
delta.AIC <- NA
}
df.res <- data.frame(delta.AIC=delta.AIC)
return(df.res)
}

## function to compute ratio of variance explained by a trait variable over the variance explained by the BATOT
fun.ratio.var.fixed.effect <- function(DF.results){
Georges Kunstler's avatar
Georges Kunstler committed
216
mat.ratio <- DF.results[, c('sumTnTfBn.abs.VAR', 'sumTfBn.VAR', 'sumTnBn.VAR',
217
                           'effect.response.var')]/
Georges Kunstler's avatar
Georges Kunstler committed
218
219
    DF.results[, 'sumBn.VAR']
names(mat.ratio) <- c('abs.dist', 'Response', 'Effect', 'Effect.Response')
220
221
222
223
return(mat.ratio)
}

### FUNCTION TO REPORT BEST MODEL PER ECOREGION AND TRAITS
Georges Kunstler's avatar
Georges Kunstler committed
224
fun.AIC <- function(id2.one, DF.results){
Georges Kunstler's avatar
Georges Kunstler committed
225
226
227
 models <- c('lmer.LOGLIN.nocomp.Tf',  'lmer.LOGLIN.simplecomp.Tf',
             'lmer.LOGLIN.HD.Tf',
             'lmer.LOGLIN.E.Tf', 'lmer.LOGLIN.R.Tf', 'lmer.LOGLIN.ER.Tf',
228
             'lmer.LOGLIN.AD.Tf')
Georges Kunstler's avatar
Georges Kunstler committed
229
 best <- as.vector(DF.results[DF.results$id2==id2.one, c('id2', 'trait', 'set',
Georges Kunstler's avatar
Georges Kunstler committed
230
231
232
233
                                  'ecocode', 'filling', 'MAT', 'MAP', 'model')])[
                           which.min(DF.results$AIC[DF.results$id2==id2.one]), ]
 AIC.all <- as.vector(DF.results[DF.results$id2==id2.one, c('AIC')])
 names(AIC.all) <- as.vector(DF.results[DF.results$id2==id2.one, c('model')])
234
 AIC.all <- AIC.all[models]-min(AIC.all)
Georges Kunstler's avatar
Georges Kunstler committed
235
 res <- data.frame((best), t(AIC.all))
Georges Kunstler's avatar
Georges Kunstler committed
236
 names(res) <- c('id2', 'trait', 'set', 'ecocode', 'filling', 'MAT', 'MAP',
Georges Kunstler's avatar
Georges Kunstler committed
237
                 'best.model', paste('AIC', models, sep='.'))
238
239
240
 return(res)
}

Georges Kunstler's avatar
Georges Kunstler committed
241
fun.AICc <- function(id2.one, DF.results){
Georges Kunstler's avatar
Georges Kunstler committed
242
243
 models <- c('lmer.LOGLIN.nocomp.Tf',  'lmer.LOGLIN.simplecomp.Tf',
             'lmer.LOGLIN.HD.Tf', 'lmer.LOGLIN.E.Tf', 'lmer.LOGLIN.R.Tf',
Georges Kunstler's avatar
Georges Kunstler committed
244
245
246
             'lmer.LOGLIN.ER.Tf', 'lmer.LOGLIN.AD.Tf')
 Deviance.all <- DF.results[DF.results$id2==id2.one, 'deviance']
 names(Deviance.all) <- DF.results[DF.results$id2==id2.one, 'model']
247
 Deviance.all <- Deviance.all[models]
Georges Kunstler's avatar
Georges Kunstler committed
248
249
 nobs.all <- DF.results[DF.results$id2==id2.one, 'nobs']
 names(nobs.all) <- DF.results[DF.results$id2==id2.one, 'model']
250
 nobs.all <- nobs.all[models]
Georges Kunstler's avatar
Georges Kunstler committed
251
 n.param <- c(2, 3, 4, 4, 4, 5, 4)
252
 AICc <-  Deviance.all+2*n.param*(nobs.all)/(nobs.all-n.param-1)
Georges Kunstler's avatar
Georges Kunstler committed
253
254
255
 id2.n <- unique(DF.results[DF.results$id2==id2.one, c('id2')])
 res <- data.frame(id2.n, models[which.min(AICc)], t(AICc), row.names=NULL)
 names(res) <- c('id2', 'best.model', models)
256
257
258
259
260
261
 return(res)
}

#################################333

### function to get the data for a given model with criteria to select
Georges Kunstler's avatar
Georges Kunstler committed
262
263
fun.select.ecoregions.trait <- function(DF.results, trait.name, model.selected,
                                      nobs.min=1000, filling.selected="species",
264
                                      threshold.delta.AIC = 100000){
265
266
267
268
DF.results[DF.results$nobs>nobs.min &
           DF.results$filling==filling.selected &
           DF.results$trait==trait.name &
           DF.results$model %in% model.selected &
Georges Kunstler's avatar
Georges Kunstler committed
269
           DF.results$delta.AIC<threshold.delta.AIC, ]
270
271
272
}

### function to get the data for a given model with criteria to select only site with competition
Georges Kunstler's avatar
Georges Kunstler committed
273
274
fun.select.ecoregions.trait.compet <- function(DF.results, trait.name, model.selected,
                                      nobs.min=1000, filling.selected="species",
275
                                      threshold.delta.AIC = 100000){
276
277
278
279
280
DF.results[DF.results$nobs>nobs.min &
           DF.results$filling==filling.selected &
           DF.results$trait==trait.name &
           DF.results$model %in% model.selected &
           DF.results$sumBn < 0.0 &
Georges Kunstler's avatar
Georges Kunstler committed
281
282
           ## DF.results$delta.AIC==0, ]
           DF.results$delta.AIC<threshold.delta.AIC, ]
283
284
285
286
287
288
289
}




#########################
##### FUNCTIONS FOR PLOTS
Georges Kunstler's avatar
Georges Kunstler committed
290
291
292
fun.plot.lmer.res.x.y <- function(model.selected, trait.name, DF.results, var.x, var.y, threshold.delta.AIC, ...){
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name, model.selected=model.selected, threshold.delta.AIC=10000000)
plot(df.selected[[var.x]], df.selected[[var.y]], ...)
293
294
}

Georges Kunstler's avatar
Georges Kunstler committed
295
296
297
298
fun.plot.lmer.res.x.y.2 <- function(model.selected, trait.name, DF.results, var.x, var.y, col.vec, pch.AIC=TRUE, cex.AIC=TRUE, col.set=TRUE, ...){
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name, model.selected=model.selected, threshold.delta.AIC=10000000)
if(pch.AIC) {pch.vec <- c(1, 16)[as.numeric(df.selected[['delta.AIC']]==0)+1]}else{pch.vec <- 1}
if(cex.AIC) {cex.vec <- c(1, 1.5)[as.numeric(df.selected[['delta.AIC']]==0)+1]}else{cex.vec <- 1}
299
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
Georges Kunstler's avatar
Georges Kunstler committed
300
301
302
plot(df.selected[[var.x]], df.selected[[var.y]],
     pch=pch.vec,
     cex=cex.vec,
Georges Kunstler's avatar
Georges Kunstler committed
303
     col=col.vec2, ...)
304
305
}

Georges Kunstler's avatar
Georges Kunstler committed
306
307
308
fun.plot.lmer.res.boxplot <- function(model.selected, trait.name, DF.results, var.y, ...){
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name, model.selected=model.selected)
boxplot(df.selected[[var.y]], ...)
309
310
}

Georges Kunstler's avatar
Georges Kunstler committed
311
312
313
fun.plot.lmer.res.boxplot.compare.model <- function(model.selected, trait.name, DF.results, var.y, ...){
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name, model.selected=model.selected)
boxplot(df.selected[[var.y]]~df.selected[['model']], ...)
314
315
316
317
318
}




Georges Kunstler's avatar
Georges Kunstler committed
319
320
321
322
fun.plot.error.bar <- function(x,  y,  sd,  small.bar = (max(x) - min(x))/40){
segments( x,  unlist(y - 1.96*sd),  x ,  unlist(y +1.96*sd))
segments( x - small.bar,  unlist(y - 1.96*sd),  x + small.bar ,  unlist(y -1.96*sd))
segments( x - small.bar,  unlist(y + 1.96*sd),  x + small.bar,  unlist(y +1.96*sd))
323
324
}

Georges Kunstler's avatar
Georges Kunstler committed
325
326
327

fun.plot.param.error.bar <- function(df.selected, var.x, param, small.bar, col.vec){
segments( df.selected[[var.x]], df.selected[[param]] - 1.96*df.selected[[paste(
Georges Kunstler's avatar
Georges Kunstler committed
328
    param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
329
330
331
          df.selected[[var.x]], df.selected[[param]] + 1.96*df.selected[[paste(
              param, "Std.Error", sep=".")]], col=col.vec)
segments( df.selected[[var.x]]-small.bar, df.selected[[param]] - 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
332
         df.selected[[paste(param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
333
334
335
          df.selected[[var.x]]+small.bar, df.selected[[param]] - 1.96*
         df.selected[[paste(param, "Std.Error", sep=".")]], col=col.vec)
segments( df.selected[[var.x]]-small.bar, df.selected[[param]] + 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
336
         df.selected[[paste(param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
337
338
339
340
          df.selected[[var.x]]+small.bar, df.selected[[param]] + 1.96*
         df.selected[[paste(param, "Std.Error", sep=".")]], col=col.vec)
}

Georges Kunstler's avatar
Georges Kunstler committed
341
342
343
fun.plot.all.param.x.y.c <- function(model.selected, trait.name, DF.results, var.x,
                                     params, small.bar, threshold.delta.AIC,
                                     col.vec, col.set=TRUE, ylim=NA,
Georges Kunstler's avatar
Georges Kunstler committed
344
                                     ylim.same=FALSE, add.name = FALSE,  ...){
Georges Kunstler's avatar
Georges Kunstler committed
345
346
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name,
                                           model.selected=model.selected,
347
                                           threshold.delta.AIC=threshold.delta.AIC)
348
if(col.set) {col.vec2 <- col.vec[as.character(df.selected[['set']])]}else{
349
350
    col.vec2 <- 1}
if(!ylim.same) {ylim <- range(c(df.selected[[params[1]]] - 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
351
                           df.selected[[paste(params[1], "Std.Error", sep=".")]],
352
                df.selected[[params[1]]] + 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
353
354
               df.selected[[paste(params[1], "Std.Error", sep=".")]]), na.rm=TRUE)}
plot(df.selected[[var.x]], df.selected[[params[1]]], col=col.vec2, ylim=ylim, ...)
Georges Kunstler's avatar
Georges Kunstler committed
355
356
if (add.name) text(df.selected[[var.x]],  df.selected[[params[1]]],
                   labels = paste(df.selected[['set']],
Georges Kunstler's avatar
Georges Kunstler committed
357
                       df.selected[['ecocode']]),  cex = 0.5,  pos = 2)
Georges Kunstler's avatar
Georges Kunstler committed
358
fun.plot.param.error.bar(df.selected, var.x, param=params[1],
Georges Kunstler's avatar
Georges Kunstler committed
359
360
361
362
363
364
                         small.bar, col=col.vec2)
}


fun.plot.param.error.bar.std <- function(df.selected, var.x, param, small.bar, col.vec){
segments( df.selected[[var.x]], df.selected[[param]] + df.selected[['sumBn']]  - 1.96*df.selected[[paste(
Georges Kunstler's avatar
Georges Kunstler committed
365
    param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
366
367
368
          df.selected[[var.x]], df.selected[[param]]  + df.selected[['sumBn']] + 1.96*df.selected[[paste(
              param, "Std.Error", sep=".")]], col=col.vec)
segments( df.selected[[var.x]]-small.bar, df.selected[[param]] + df.selected[['sumBn']] - 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
369
         df.selected[[paste(param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
370
371
372
          df.selected[[var.x]]+small.bar, df.selected[[param]]  + df.selected[['sumBn']] - 1.96*
         df.selected[[paste(param, "Std.Error", sep=".")]], col=col.vec)
segments( df.selected[[var.x]]-small.bar, df.selected[[param]]  + df.selected[['sumBn']] + 1.96*
Georges Kunstler's avatar
Georges Kunstler committed
373
         df.selected[[paste(param, "Std.Error", sep=".")]],
Georges Kunstler's avatar
Georges Kunstler committed
374
375
376
377
          df.selected[[var.x]]+small.bar, df.selected[[param]]  + df.selected[['sumBn']] + 1.96*
         df.selected[[paste(param, "Std.Error", sep=".")]], col=col.vec)
}

Georges Kunstler's avatar
Georges Kunstler committed
378
379
380
381
fun.plot.all.param.x.y.std <- function(model.selected, trait.name,
                                       DF.results, var.x,
                                     params, small.bar, threshold.delta.AIC,
                                     col.vec, col.set=TRUE, ylim=NA,
Georges Kunstler's avatar
Georges Kunstler committed
382
                                     ylim.same=FALSE, add.name = FALSE,  ...){
Georges Kunstler's avatar
Georges Kunstler committed
383
384
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name,
                                           model.selected=model.selected,
385
                                           threshold.delta.AIC=threshold.delta.AIC)
386
if(col.set) {col.vec2 <- col.vec[df.selected[['set']]]}else{
387
    col.vec2 <- 1}
Georges Kunstler's avatar
Georges Kunstler committed
388
plot(df.selected[[var.x]], df.selected[[params[1]]] + df.selected[['sumBn']],
Georges Kunstler's avatar
Georges Kunstler committed
389
     col = col.vec2,  ...)
Georges Kunstler's avatar
Georges Kunstler committed
390
391
392
if (add.name) text(df.selected[[var.x]],
                   df.selected[[params[1]]] + df.selected[['sumBn']],
                   labels = paste(df.selected[['set']],
Georges Kunstler's avatar
Georges Kunstler committed
393
                       df.selected[['ecocode']]),  cex = 0.5,  pos = 2)
Georges Kunstler's avatar
Georges Kunstler committed
394
fun.plot.param.error.bar.std(df.selected, var.x, param=params[1],
Georges Kunstler's avatar
Georges Kunstler committed
395
                         small.bar, col=col.vec2)
396
397
398
399
400

}



Georges Kunstler's avatar
Georges Kunstler committed
401
fun.plot.all.param.boxplot <- function(model.selected, trait.name, DF.results,
Georges Kunstler's avatar
Georges Kunstler committed
402
                                       params, small.bar, ...){
Georges Kunstler's avatar
Georges Kunstler committed
403
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name,
404
                                           model.selected=model.selected)
405
 if(length(params)>1){
Georges Kunstler's avatar
Georges Kunstler committed
406
407
408
DF.t <- data.frame(param=rep(names(params), each=nrow(df.selected)), value=c(
                            df.selected[[params[1]]], df.selected[[params[2]]]))
boxplot(DF.t[['value']]~DF.t[['param']], ...)
409
}else{
Georges Kunstler's avatar
Georges Kunstler committed
410
boxplot(df.selected[[params[1]]], ...)
411
412
413
}
}

Georges Kunstler's avatar
Georges Kunstler committed
414
fun.plot.all.param.er.diff.MAP <- function(model.selected, trait.name,
Georges Kunstler's avatar
Georges Kunstler committed
415
                                           DF.results, ...){
Georges Kunstler's avatar
Georges Kunstler committed
416
df.selected <- fun.select.ecoregions.trait(DF.results, trait.name=trait.name,
417
                                           model.selected=model.selected)
Georges Kunstler's avatar
Georges Kunstler committed
418
plot(df.selected[['MAP']], df.selected[['sumTnBn']]-df.selected[['sumTfBn']], ...)
419
420
421
422
}



Georges Kunstler's avatar
Georges Kunstler committed
423
fun.plot.panel.lmer.res.x.y <- function(models, traits, DF.results, var.x, var.y.l,
Georges Kunstler's avatar
Georges Kunstler committed
424
                                        express, ...){
425
426
427
    ncols = length(traits)
    nrows = length(models)
    list.models <- as.list(names(models))
Georges Kunstler's avatar
Georges Kunstler committed
428
    names(list.models) <- rep('model', length(list.models))
429
430
    DF.results$set <- factor(DF.results$set)
    col.vec <- niceColors(n=nlevels(DF.results$set))
Georges Kunstler's avatar
Georges Kunstler committed
431
    par(mfrow = c(nrows,  ncols),  mar = c(2, 2, 1, 1),  oma = c(4, 4, 4, 1) )
432
433
    for(i in 1:nrows)
        for(j in 1:ncols){
Georges Kunstler's avatar
Georges Kunstler committed
434
              fun.plot.lmer.res.x.y.2(models[i], traits[j],
Georges Kunstler's avatar
Georges Kunstler committed
435
436
                          DF.results, var.x, var.y=var.y.l[[i]], col.vec, ...)
              abline(h=0, lty=3)
437
            if(i==1 )
Georges Kunstler's avatar
Georges Kunstler committed
438
                mtext(traits[j],  side=3,  line =1)
439
            if(i==nrows)
Georges Kunstler's avatar
Georges Kunstler committed
440
                mtext(var.x,  side=1,  line =4)
441
            if(j==1)
Georges Kunstler's avatar
Georges Kunstler committed
442
                mtext(paste('Effect size', list.models[i]),  side=2,  line =4,
443
                      cex=0.9)
444
            if(i==nrows & j==ncols)
Georges Kunstler's avatar
Georges Kunstler committed
445
                legend('topright', legend=levels(DF.results$set), pch=16,
Georges Kunstler's avatar
Georges Kunstler committed
446
                       col=col.vec, bty='n', ncol=2)
447
448
449
450
        }
}


Georges Kunstler's avatar
Georges Kunstler committed
451
fun.plot.panel.lmer.res.boxplot <- function(models, traits, DF.results, var.y,
Georges Kunstler's avatar
Georges Kunstler committed
452
                                            express, ...){
453
454
455
    ncols = length(traits)
    nrows = length(models)
    list.models <- as.list(names(models))
Georges Kunstler's avatar
Georges Kunstler committed
456
457
    names(list.models) <- rep('model', length(list.models))
    par(mfrow = c(nrows,  ncols),  mar = c(2, 2, 1, 1),  oma = c(4, 4, 4, 1) )
458
459
    for(i in 1:nrows)
        for(j in 1:ncols){
Georges Kunstler's avatar
Georges Kunstler committed
460
              fun.plot.lmer.res.boxplot(models[i], traits[j],
Georges Kunstler's avatar
Georges Kunstler committed
461
462
                          DF.results, var.y[[j]], ...)
              abline(h=0, lty=3)
463
            if(i==1 )
Georges Kunstler's avatar
Georges Kunstler committed
464
                mtext(traits[j],  side=3,  line =1)
465
            if(j==1)
Georges Kunstler's avatar
Georges Kunstler committed
466
                mtext(paste("Effect size",  list.models[i], sep=" "),  side=2,
Georges Kunstler's avatar
Georges Kunstler committed
467
                      line =4, cex=0.9)
468
469
470
471
        }
}


Georges Kunstler's avatar
Georges Kunstler committed
472
fun.plot.panel.lmer.res.boxplot.compare <- function(models, traits, DF.results,
Georges Kunstler's avatar
Georges Kunstler committed
473
                                                    var.y, express, ...){
474
475
    ncols = length(traits)
    list.models <- as.list(names(models))
Georges Kunstler's avatar
Georges Kunstler committed
476
477
    names(list.models) <- rep('model', length(list.models))
    par(mfrow = c(1,  ncols),  mar = c(2, 2, 1, 1),  oma = c(4, 4, 4, 1) )
478
    for(j in 1:ncols){
Georges Kunstler's avatar
Georges Kunstler committed
479
              fun.plot.lmer.res.boxplot.compare.model(models, traits[j],
Georges Kunstler's avatar
Georges Kunstler committed
480
481
482
                          DF.results, var.y, ...)
              abline(h=0, lty=3)
                mtext(traits[j],  side=3,  line =1, cex=2)
483
            if(j==1 )
Georges Kunstler's avatar
Georges Kunstler committed
484
                mtext("Effect size",  side=2,  line =4, cex=1.5)
485
486
487
        }
}

Georges Kunstler's avatar
Georges Kunstler committed
488
489
490
fun.plot.panel.lmer.parameters.c <- function(models, traits, DF.results, var.x,
                                             list.params, threshold.delta.AIC,
                                             small.bar=10, ylim=NA,
Georges Kunstler's avatar
Georges Kunstler committed
491
492
                                             ylim.same=FALSE,
                                             col.vec2=col.vec, ...){
493
494
    ncols = length(traits)
    nrows = length(models)
Georges Kunstler's avatar
Georges Kunstler committed
495
    par(mfrow = c(nrows,  ncols),  mar = c(2, 2, 1, 1),  oma = c(4, 4, 4, 1) )
496
    DF.results$set <- factor(DF.results$set)
497
498

## ### TO COMPARE THE PARAMTERS WE NEED TO DIVIDE THE ABS>DIST per two
499
    ## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
500
501
    for(i in 1:nrows)
        for(j in 1:ncols){
Georges Kunstler's avatar
Georges Kunstler committed
502
              try(fun.plot.all.param.x.y.c(models[i], traits[j], DF.results,
Georges Kunstler's avatar
Georges Kunstler committed
503
504
505
506
507
                                           var.x,
                                       params=list.params[[i]],
                                       small.bar=small.bar,
                                       threshold.delta.AIC=threshold.delta.AIC,
                                       col.vec=col.vec2, col.set=TRUE,
Georges Kunstler's avatar
Georges Kunstler committed
508
509
                                       ylim=ylim, ylim.same=ylim.same, ...))
              try(abline(h=0, lty=3))
510
              if(length(list.params[[i]])>1)
Georges Kunstler's avatar
Georges Kunstler committed
511
512
                  legend("topright", names(list.params[[i]]),
                         pch=rep(1, length(list.params[[i]])),
Georges Kunstler's avatar
Georges Kunstler committed
513
                         col=1:length(list.params[[i]]), bty='n', cex=1)
514
            if(i==1 )
Georges Kunstler's avatar
Georges Kunstler committed
515
                mtext(traits[j],  side=3,  line =1)
516
            if(i==nrows)
Georges Kunstler's avatar
Georges Kunstler committed
517
                mtext(var.x,  side=1,  line =4)
518
            if(j==1)
Georges Kunstler's avatar
Georges Kunstler committed
519
520
                mtext(paste('param', names(models)[i]),  side=2,
                      line =4, cex=1)
521
            if(i==nrows & j==ncols)
Georges Kunstler's avatar
Georges Kunstler committed
522
                legend('topleft', legend=levels(DF.results$set), pch=16,
Georges Kunstler's avatar
Georges Kunstler committed
523
                       col=col.vec2[levels(DF.results$set)], bty='n', ncol=2)
524
525
526
527
528
        }
}



Georges Kunstler's avatar
Georges Kunstler committed
529
530
531
fun.plot.panel.lmer.parameters.boxplot <- function(models, traits, DF.results,
                                                   list.params,
                                                   small.bar=10, ...){
532
533
    ncols = length(traits)
    nrows = length(models)
Georges Kunstler's avatar
Georges Kunstler committed
534
    par(mfrow = c(nrows,  ncols),  mar = c(2, 2, 1, 1),  oma = c(4, 4, 4, 1) )
535
536
537
538
### TO COMPARE THE PARAMTERS WE NEED TO DIVIDE THE ABS>DIST per two
    ## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
    for(i in 1:nrows)
        for(j in 1:ncols){
Georges Kunstler's avatar
Georges Kunstler committed
539
540
541
542
              fun.plot.all.param.boxplot(models[i], traits[j], DF.results,
                                         params=list.params[[i]],
                                         small.bar=small.bar, ...)
              abline(h=0, lty=3)
543
            if(i==1 )
Georges Kunstler's avatar
Georges Kunstler committed
544
                mtext(traits[j],  side=3,  line =1)
545
            if(j==1)
Georges Kunstler's avatar
Georges Kunstler committed
546
                mtext(paste('param', names(models)[i]),  side=2,  line =4, cex=1)
547
548
549
        }
}

550
551
552


## create a data base from the global data
Georges Kunstler's avatar
Georges Kunstler committed
553
fun.merge.results.global <-  function(list.res,
554
555
556
                                      names.param ){
  names.param.vec <- unique(names.param)
  df.details  <- data.frame(list.res$files.details[1:4],
557
                  list.res$lmer.summary[1:6])
558
559
560
561
562
563
564
565
566
567
568
569
  dat.t <- data.frame(matrix(rep(NA, 3 * length(names.param.vec)), nrow = 1,
                             ncol = 3 * length(names.param.vec)))
  names(dat.t) <-  c(names.param.vec,
                     paste(names.param.vec, "Std.Error", sep="."),
                     paste(names.param.vec, "VAR", sep="."))

dat.t[, names.param2[ names(list.res$lmer.summary$fixed.coeff.E)]] <-
                                             list.res$lmer.summary$fixed.coeff.E
dat.t[, paste(names.param2[  names(list.res$lmer.summary$fixed.coeff.E)], "Std.Error", sep=".")] <-
                                             list.res$lmer.summary$fixed.coeff.Std.Error
dat.t[, paste(names.param2[  names(list.res$lmer.summary$fixed.coeff.E)], "VAR", sep=".")] <-
                                             list.res$lmer.summary$fixed.var
Georges Kunstler's avatar
Georges Kunstler committed
570
571
res <- data.frame(df.details,  dat.t)
return(res)
572
}
573
574
575
576
577
578



#################################
## compute a global AIC
fun.global.aic <- function(DF.results){
Georges Kunstler's avatar
Georges Kunstler committed
579
DF.results <- DF.results[DF.results$nobs>1000,  ]
580
581
# select set ecocode with more than 1000 obs
DF.results <- DF.results[DF.results$id2 %in%  names(table(DF.results$id2))[
Georges Kunstler's avatar
Georges Kunstler committed
582
                         table(DF.results$id2) == 7],   ]
583
584
# select set ecocode with 7 model tested
## species
Georges Kunstler's avatar
Georges Kunstler committed
585
DF.results.sp <- DF.results[DF.results$filling == 'species',  ]
586
                               # select set ecocode with more than 1000 obs
Georges Kunstler's avatar
Georges Kunstler committed
587
588
589
test.same.n.model.ecoregion <- apply(table(DF.results.sp$trait,
                                           DF.results.sp$model),
                                     MARGIN = 1,
590
591
                                     function(x) all(x == x[1]))
if(!all(test.same.n.model.ecoregion))
Georges Kunstler's avatar
Georges Kunstler committed
592
    stop(paste('error not the same number of ecoregion for traits',
593
               names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
Georges Kunstler's avatar
Georges Kunstler committed
594
595
596
597
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,
Georges Kunstler's avatar
Georges Kunstler committed
598
                                           DF.results.sp$model),  FUN = length))
599
## genus
Georges Kunstler's avatar
Georges Kunstler committed
600
DF.results.ge <- DF.results[DF.results$filling == 'genus',  ]
601
   # select set ecocode with more than 1000 obs
Georges Kunstler's avatar
Georges Kunstler committed
602
603
604
test.same.n.model.ecoregion <- apply(table(DF.results.ge$trait,
                                           DF.results.ge$model),
                                     MARGIN = 1,
605
606
                                     function(x) all(x == x[1]))
if(!all(test.same.n.model.ecoregion))
Georges Kunstler's avatar
Georges Kunstler committed
607
    stop(paste('error not the same number of ecoregion for traits',
608
             names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
Georges Kunstler's avatar
Georges Kunstler committed
609
610
611
612
613
614
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),
615
616
    FUN = length))
return(list(sp = list.sp, ge = list.ge))
Georges Kunstler's avatar
Georges Kunstler committed
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
}


## start code to predict regression line with CI

library(lme4)
library(ggplot2) # Plotting
data("Orthodont",package="MEMSS")
fm1 <- lmer(
    formula = distance ~ age*Sex + (age|Subject)
    , data = Orthodont)


fun.prepare.data.pred <-  function(data, var, n.length = 100){
if( !is.numeric(data[[var]]) ) stop('var need to be a numeric variables')
q.var <- quantile(data[[var]], probs = c(0.05, 0.95))
seq.var <- seq(q.var[1], q.var[2], length = n.length)
vec.class <- sapply(data, class)
vec.num <- names(vec.class)[vec.class == 'numeric']
means.num <- lapply(vec.num, function(x, data) mean(data[[x]], na.rm =TRUE),
                    data)
names(means.num) <- vec.num
vec.no.num <- names(vec.class)[vec.class != 'numeric']
values.no.num <- lapply(vec.no.num, function(x, data) data[[x]][1],
                        data)
names(values.no.num) <- vec.no.num
list.res <- c(means.num,values.no.num)
list.res[[var]] <- seq.var
new.data <- expand.grid(list.res)
return(new.data)
}




fun.predict.ci.lmer <- function(fm, newdata){
mm <- model.matrix(terms(fm1),newdat)
newdat$distance <- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
newdat <- data.frame(newdat,
                     plo = newdat$distance-1.96*sqrt(pvar1),
                     phi = newdat$distance+1.96*sqrt(pvar1))
return(newdat)
}


newdat <- fun.prepare.data.pred(data = newdat, 'age')
newdat <-  fun.predict.ci.lmer(fm1, newdat)

## TODO TES THAT !

#plot confidence

g0 <- ggplot(newdat, aes(x=age, y=distance, colour=Sex))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
    opts(title="CI based on fixed-effects uncertainty ONLY")