Commit d9f7a50a by Kunstler Georges

### rerun deid

parent 3bd7fd04
 ... ... @@ -3,80 +3,45 @@ library(lme4) library(MuMIn) read.lmer.output <- function(file.name){ tryCatch(readRDS(file.name), error = function(cond)return(NULL)) # Choose a return value in case of error read.lmer.output <- function(file.name){ #OK tryCatch(readRDS(file.name), error = function(cond)return(NULL)) # Choose a return value in case of error } summarise.lmer.output <- function(x){ list( nobs = nobs(x), R2m =r.squaredGLMM(x)[['R2m']], R2c =r.squaredGLMM(x)[['R2c']], 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)) summarise.lmer.all.output <- function(x, random.name ){#OK list(nobs = nobs(x), R2m =r.squaredGLMM(x)[['R2m']], R2c =r.squaredGLMM(x)[['R2c']], AIC = AIC(x), deviance = deviance(x), conv = x@optinfo\$conv, fixed.coeff.E = fixef(x), fixed.coeff.Std.Error = sqrt(diag(vcov(x))), set.BLUP = coef(x)[[random.name]]) } summarise.lmer.all.output <- function(x, random.name ){ list( nobs = nobs(x), R2m =r.squaredGLMM(x)[['R2m']], R2c =r.squaredGLMM(x)[['R2c']], 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)[[random.name]]) } summarise.lmer.output.list <- function(f ){ list.output.lmer <- read.lmer.output(f) output.lmer <- list.output.lmer[['output']] relgrad <- list.output.lmer[['relgrad']] list.sd <- list.output.lmer[['lits.sd']] if(!is.null(output.lmer)){ res <- list(files.details = files.details(f) , lmer.summary = summarise.lmer.output( output.lmer), terms = terms(output.lmer), vcov = vcov(output.lmer), relgrad = relgrad, list.sd = list.sd) }else{ res <- NULL } return(res) } summarise.lmer.output.all.list <- function(f ){ summarise.lmer.output.all.list <- function(f){#OK list.output.lmer <- read.lmer.output(f) output.lmer <- list.output.lmer[['output']] relgrad <- list.output.lmer[['relgrad']] list.sd <- list.output.lmer[['list.sd']] if(!is.null(output.lmer)){ details <- files.details.all(f) source(file.path('R','analysis', 'model.lmer', paste('model',details\$model,'R', sep = '.'))) model <- load.model() res <- list(files.details = details, lmer.summary = summarise.lmer.all.output( output.lmer, random.name = model\$var.BLUP), lmer.summary = summarise.lmer.all.output(output.lmer, random.name = model\$var.BLUP), terms = terms(output.lmer), vcov = vcov(output.lmer), relgrad = relgrad, list.sd = list.sd) }else{ res <- NULL res <- NULL } rm(list.output.lmer) gc() ... ... @@ -85,16 +50,8 @@ summarise.lmer.output.all.list <- function(f ){ files.details <- function(x){ s <- data.frame(t(strsplit(x, "/", fixed = TRUE)[[1]]), stringsAsFactors = FALSE) names(s) <- c("d1", "d2", "data.type", "ecocode", "trait", "model", "file" ) s[-(1:2)] } files.details.all <- function(x){ files.details.all <- function(x){#OK s <- data.frame(t(strsplit(x, "/", fixed = TRUE)[[1]]), stringsAsFactors = FALSE) names(s) <- c("d1", "d2", "data.type", "trait", ... ... @@ -102,89 +59,38 @@ files.details.all <- function(x){ s[-(1:2)] } BLUP.CI <- function(fm, var.BLUP = 'species.id'){ ## NOT WORKING WHEN MULTIPLE term per factor cV <- ranef(fm, condVar = TRUE, whichel = var.BLUP) return(cV) } #### R squared functions # 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/ Rsquared.glmm.lmer <- function(i){ # 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 VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1]))) # Get residual variance VarResid <- attr(VarCorr(i), "sc")^2 # 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 Rsquared.mat <- data.frame(Class = class(i), Family = "Gaussian", Marginal = Rm, Conditional = Rc, AIC = AIC(i)) return(Rsquared.mat) } variance.fixed.glmm.lmer <- function(i){ # Get variance of for each fixed effects by multiplying coefficients by design matrix var.vec <- apply(fixef(i) * t(i@pp\$X), MARGIN = 1, var) # 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 VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1]))) # Get residual variance VarResid <- attr(VarCorr(i), "sc")^2 var.vec <- var.vec/(VarF+VarRand+VarResid) names(var.vec) <- paste(names(var.vec), "VAR", sep = ".") return(var.vec) } variance.fixed.glmm.lmer.effect.and.response <- function(i){ if(sum(c("sumTfBn", "sumTnBn") %in% names(fixef(i)))==2){ # Get variance of for each fixed effects by multiplying coefficients by design matrix var.vec <- var(as.vector(fixef(i)[c("sumTfBn", "sumTnBn")] %*% t(i@pp\$X[, c("sumTfBn", "sumTnBn")]))) # 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 VarRand <- colSums(do.call(rbind, lapply(VarCorr(i), function(j) j[1]))) # Get residual variance VarResid <- attr(VarCorr(i), "sc")^2 var.vec <- var.vec/(VarF+VarRand+VarResid) }else{ var.vec <- NA } names(var.vec) <- paste("effect.response", "VAR", sep = ".") return(var.vec) } #### Function to proceed all output and save in a list format.all.output.lmer <- function(file.name, list.file.name, models, traits = c("SLA", "Wood.density", "Max.height", "Leaf.N", "Seed.mass"), data.type = "simple"){ files <- c() for (trait in traits){ for(model in models){ source(model, local = TRUE) model.obj <- load.model() pathout <- output.dir('lmer', model.obj\$name, trait, data.type) files <- c(files,file.path(pathout,file.name)) } } ## function to turn lmer output from list to DF 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))] <- list.res\$lmer.summary\$fixed.coeff.E dat.t[, length(names.param)+match(names(list.res\$lmer.summary\$fixed.coeff.E), names(dat.t))] <- list.res\$lmer.summary\$fixed.coeff.Std.Error dat.t[, match(names(list.res\$lmer.summary\$fixed.var), names(dat.t))] <- list.res\$lmer.summary\$fixed.var res <- data.frame(list.res\$files.details, list.res\$lmer.summary[1:7], dat.t) return(res) out <- lapply(files, summarise.lmer.output.all.list) names(out) <- lapply(lapply(files,files.details.all), function(x) paste(as.vector(x[names(x) != 'file']), collapse="_")) ### remove missing out <- out[!unlist(lapply(out,FUN=function(x) is.null(x\$lmer.summary)))] saveRDS(out,file=file.path('output' , list.file.name)) rm(out) gc() } ... ... @@ -192,10 +98,13 @@ return(res) ########################### ### FUNCTIONS TO PLOTS RESULTS fun.plot.error.bar.horiz <- function(x, y, sd, small.bar = (max(x) - min(x))/40, ...){ segments( unlist(x - 1.96*sd), y, unlist(x +1.96*sd), y, ...) segments( unlist(x - 1.96*sd), y-small.bar, unlist(x -1.96*sd), y+small.bar, ...) segments( unlist(x + 1.96*sd), y-small.bar, unlist(x +1.96*sd), y+small.bar, ...) fun.plot.error.bar.horiz <- function(x, y, sd, #OK small.bar = (max(x) - min(x))/40, ...){ segments(unlist(x - 1.96*sd), y, unlist(x +1.96*sd), y, ...) segments(unlist(x - 1.96*sd), y-small.bar, unlist(x -1.96*sd), y+small.bar, ...) segments(unlist(x + 1.96*sd), y-small.bar, unlist(x +1.96*sd), y+small.bar, ...) } fun.col.param <- function(){ ... ... @@ -952,88 +861,13 @@ ggplot(DF, aes_string(x = var.x, y = var.y, colour = var.quant)) + fun.axis.one.by.one <- function(x, side, labels, cols.vec , cex.axis = 2){ fun.axis.one.by.one <- function(x, side, labels, cols.vec , cex.axis = 2){ # OK axis(side, x, labels = bquote(.(labels[x])), las = 1, cex.axis = cex.axis, mgp = c(1.5,0.55,0),col.axis = cols.vec[x]) } plot.param <- function(list.res, model = 'lmer.LOGLIN.ER.AD.Tf', traits = c('Wood.density', 'SLA', 'Max.height'), traits.names = c(Wood.density= 'Wood density', SLA = 'Specific leaf area', Max.height = 'Maximum height'), param.vec = c("Tf","sumBn", "sumTnBn", "sumTfBn", "sumTnTfBn.abs"), param.names = c('Direct trait', 'Compet', 'Compet effect x trait', 'Compet response x trait', 'Compet x trait dissimilarity'), param.print = 1:5, traits.print = 1:3, col.names = fun.col.param(), data.type = "simple", add.param.descrip.TF = TRUE, intra.TF = FALSE, ...){ if(!intra.TF) x.line <- -0.63 if(intra.TF) x.line = -0.7 m <- matrix(c(1:3), 1, 3) big.m <- 2.8 small.m <- 0.42 wid <- c(big.m ,0, small.m) + rep((14-big.m-small.m)/3, each= 3) layout(m, widths=wid) for (i in traits[traits.print]){ list.temp <- list.res[[paste(data.type, "_", i , "_", model, sep = '')]]\$lmer.summary param.mean <- list.temp\$fixed.coeff.E[param.vec] param.mean[!names(param.mean) %in% c("Tf", "sumTfBn")] <- -param.mean[!names(param.mean) %in% c("Tf", "sumTfBn")] param.std <- list.temp\$fixed.coeff.Std.Error names(param.std) <- names(list.temp\$fixed.coeff.E) param.std <- param.std[param.vec] param.BLUP <- list.temp\$set.BLUP if(i == traits[1]) { par(mai=c(1.2, big.m,0.6,0), xpd = TRUE) }else{ par(mai=c(1.2,0,0.6,0), xpd = TRUE) } if(i == traits[length(traits)]) { par(mai=c(1.2,0,0.6,small.m), xpd = TRUE) } plot(param.mean[param.print], (1:length(param.vec))[param.print], yaxt = 'n', xlab = NA, ylab = NA, pch = 16, cex = 2, cex.lab = 1.7, cex.axis = 1.5, ylim = range(1-0.21, length(param.vec)+0.21), ...) mtext(traits.names[i], side=3, cex =1.7, line = 1) box(lwd= 2) lines(c(0, 0), c(par()\$usr[3], par()\$usr[4])) if(i == traits[1]) {lapply(1:length(param.vec), fun.axis.one.by.one, side = 2, labels = param.names, cols.vec = col.names[param.vec]) if(add.param.descrip.TF){ fun.param.descrip(c(-0.1,0.1), length(param.vec), x.line) } } if(i == traits[2]){ mtext('Standardized coefficients', side=1, cex =1.7, line = 4) print(i) } fun.plot.error.bar.horiz(param.mean[param.print], (1:length(param.vec))[param.print], param.std[param.print]) } } ... ... @@ -1048,30 +882,6 @@ extract.param <- function(trait, list.res, return(param.mean) } extract.param.b <- function(trait, list.res, model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species', param.vec = c("logD", "Tf","sumBn", "sumTnBn", "sumTfBn", "sumTnTfBn.abs"), data.type = 'simple', biomes = biomes.factor.selected() ){ require(reshape2) list.temp <- list.res[[paste0(data.type, "_", trait , "_", model)]] res <- lapply(param.vec, fun.get.fixed.biomes, list.temp, biomes) param <- do.call('rbind', lapply(res, function(x) x[[1]])) rownames(param) <- param.vec vec.Neg <- !row.names(param) %in% c("(Intercept)", "logD", "Tf", 'MAT', 'MAP', "sumTfBn") param[vec.Neg,] <- - param[vec.Neg,] df.param <- data.frame(param=rownames(param), param) df.melt <- melt(df.param, id.vars = 'param') param.sd <- do.call('rbind', lapply(res, function(x) x[[2]])) df.param.sd <- data.frame(param=rownames(param), param.sd) names(df.param.sd) <- names(df.param) df.melt.sd <- melt(df.param.sd, id.vars = 'param') return(list(df.melt, df.melt.sd)) } extract.param.sd <- function(trait, list.res, ... ... @@ -1081,71 +891,67 @@ extract.param.sd <- function(trait, list.res, data.type = 'simple'){ list.temp <- list.res[[paste0(data.type, "_", trait , "_", model)]]\$lmer.summary param.sd <- list.temp\$fixed.coeff.Std.Error names(param.sd) <- names(list.temp\$fixed.coeff.E) return(param.sd[param.vec]) param.sd <- list.temp\$fixed.coeff.Std.Error names(param.sd) <- names(list.temp\$fixed.coeff.E) return(param.sd[param.vec]) } extract.R2c <- function(trait, list.res, extract.R2c <- function(trait, list.res, #OK model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species', data.type = 'simple'){ list.temp <- list.res[[paste0(data.type,"_", trait , "_", model)]]\$lmer.summary return(list.temp\$R2c) return(list.temp\$R2c) } extract.R2m <- function(trait, list.res, extract.R2m <- function(trait, list.res, #OK model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species', data.type = 'simple'){ list.temp <- list.res[[paste0(data.type,"_", trait , "_", model)]]\$lmer.summary return(list.temp\$R2m) } extract.AIC <- function(trait, list.res, model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species', data.type = 'simple'){ list.temp <- list.res[[paste0(data.type,"_", trait , "_", model)]]\$lmer.summary return(list.temp\$AIC) return(list.temp\$R2m) } ## get fixed biomes fun.get.fixed.biomes <- function(var, list, fun.get.fixed.biomes <- function(var, list, #OK biomes.vec ){ param <- list\$lmer.summary\$fixed.coeff.E remaining <- grep(paste0(':',var), names(param)) boreal <- seq_len(length(param))[names(param) == var] select <- c(boreal, remaining) ## create design matrix ff <- ~ biomes.id mm <- model.matrix(ff, data.frame(biomes.id = biomes.vec)) # compute mean and std based on Bolker post param.biomes <- drop(mm %*% param[select]) names(param.biomes) <- biomes.vec std.biomes <- sqrt(diag(mm %*% tcrossprod(list\$vcov[select, select],mm))) return(list(fixed.biomes = param.biomes, fixed.biomes.std = std.biomes)) } biomes.factor.selected <- function(idx.select = c(4:7, 9)){ all.fact <- factor(c("Subtropical desert", "Temperate grassland desert", "Tundra", "Boreal forest", "Temperate forest", "Temperate rain forest", "Tropical forest savanna", "Tropical rain forest", "Woodland shrubland")) levels.name <- levels(all.fact) levels.name[levels.name == 'Tundra'] <- 'Boreal forest' levels.name[levels.name == 'Subtropical desert'] <- 'Temperate grassland desert' levels.name[levels.name == "Tropical rain forest"] <- "Tropical forest savanna" levels(all.fact) <- levels.name param <- list\$lmer.summary\$fixed.coeff.E remaining <- grep(paste0(':',var), names(param)) boreal <- seq_len(length(param))[names(param) == var] select <- c(boreal, remaining) ## create design matrix ff <- ~ biomes.id mm <- model.matrix(ff, data.frame(biomes.id = biomes.vec)) # compute mean and std based on Bolker post param.biomes <- drop(mm %*% param[select]) names(param.biomes) <- biomes.vec std.biomes <- sqrt(diag(mm %*% tcrossprod(list\$vcov[select, select],mm))) return(list(fixed.biomes = param.biomes, fixed.biomes.std = std.biomes)) } biomes.factor.selected <- function(idx.select = c(4:7, 9)){#ok all.fact <- factor(c("Subtropical desert", "Temperate grassland desert", "Tundra", "Boreal forest", "Temperate forest", "Temperate rain forest", "Tropical forest savanna", "Tropical rain forest", "Woodland shrubland")) levels.name <- levels(all.fact) levels.name[levels.name == 'Tundra'] <- 'Boreal forest' levels.name[levels.name == 'Subtropical desert'] <- 'Temperate grassland desert' levels.name[levels.name == "Tropical rain forest"] <- "Tropical forest savanna" levels(all.fact) <- levels.name return(all.fact[idx.select]) } fun.biomes.names <- function(){ fun.biomes.names <- function(){#OK c("Subtropical desert" = "desert", "Temperate grassland desert" = "desert", "Tundra" = "tundra", ... ... @@ -1158,231 +964,99 @@ c("Subtropical desert" = "desert", } plot.param.biomes.fixed <- function(list.res, model, biomes.names = fun.biomes.names(), biomes = biomes.factor.selected(), traits = c('Wood.density', 'SLA', 'Max.height'), traits.names = c(Wood.density= 'Wood density', SLA = 'Specific leaf area', Max.height = 'Maximum height'), param.vec = c("Tf", "sumTnBn", "sumTfBn", "sumTnTfBn.abs"), param.names = c('Direct trait', 'Compet effect x trait', 'Compet response x trait', 'Compet x trait dissimilarity'), param.print = 1:5, col.names = fun.col.param() , data.type = "all.no.log", col.vec, pch.vec, names.bio, legend.pos = 1, add.param.descrip.TF = 1, ...){ if(! add.param.descrip.TF %in% 0:2) stop("add.param.descrip.TF need to be in 0 1 2") col.vec[2] <- col.vec[1] biomes.c <- as.character(biomes) big.m <- 3.0 small.m <- 0 legend.m <- 1.9 if (legend.pos==1){ wid <- c(big.m ,0, small.m) + rep((14-big.m-small.m-legend.m)/3, each= 3) m <- matrix(c(1:4), 1, 4) layout(m, widths=c(wid, legend.m)) } if (legend.pos==2){ wid <- c(big.m ,0, small.m) + rep((14-big.m-small.m)/3, each= 3) m <- matrix(c(1:3, 4, 4, 4), 2, 3, byrow = TRUE) layout(m, widths=c(wid), height = c(5,1)) } if(legend.pos==0) { wid <- c(big.m ,0, small.m) + rep((14-big.m-small.m)/3, each= 3) m <- matrix(c(1:3), 1, 3) layout(m, widths=c(wid)) } for (i in traits){ list.temp <- list.res[[paste0(data.type, "_", i , "_", model)]] ## NEED VAR for (n.vars in seq_len(length(param.vec))){ list.fixed <- fun.get.fixed.biomes(param.vec[n.vars], list.temp, biomes.vec = biomes) param.mean <- list.fixed\$fixed.biomes param.std <- list.fixed\$fixed.biomes.std if(i == traits[1]) { par(mai=c(0.8, big.m,0.7,0), xpd = TRUE) }else{ par(mai=c(0.8,0,0.7,0), xpd = TRUE) } if(i == traits[length(traits)]) { par(mai=c(0.8,0,0.7,small.m), xpd = TRUE) } seq.jitter <- seq(25, -25, length.out = length(biomes))/120 if(n.vars == 1){ plot(param.mean, seq.jitter+n.vars, yaxt = 'n', xlab = NA, ylab = NA, , pch = pch.vec[biomes.c] , cex = 2, cex.axis = 1.5, cex.lab = 1.5, ylim = range(1-0.21, length(param.vec)+0.21), col = col.vec[biomes.c], ...) if(i == traits[2]) mtext('Standardized coefficients', side=1, cex =1.5, line = 4) } if(n.vars != 1){ points(param.mean, seq.jitter+n.vars, col = col.vec[biomes.c], pch = pch.vec[biomes.c], cex = 2)} mtext(traits.names[i], side=3, cex =1.5, line = 1) box(lwd= 2) lines(c(0, 0), c(1-0.35, length(param.vec)+0.35), lwd= 1) if(i == traits[1]) { lapply(1:length(param.vec), fun.axis.one.by.one, side = 2, labels = param.names, cols.vec = col.names[param.vec], cex.axis = 3) if(add.param.descrip.TF == 2){