Skip to content

GitLab

Commit f9153d45 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

start global analysis and look at shannon

parent 1bffcab2
Expand all Hide whitespace changes
Inline Side-by-side
R/analysis/glmer.global.std.output.R 0 → 100644
View file @ f9153d45
#!/usr/bin/env Rscript
source("R/analysis/glmer.output-fun.R")
files <- list.files("output/glmer", recursive=TRUE, full.names =TRUE, pattern = "glmer.results.global.rds")
out <- lapply(files, summarise.glmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/glmer.global.std.list.out.rds')
names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$glmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
saveRDS(DF.results,file='output/glmer.global.std.DF.rds')
R/analysis/glmer.output.R
View file @ f9153d45
#!/usr/bin/env Rscript
source("R/analysis/glmer.output-fun.R")
files <- list.files("output/glmer", recursive=TRUE, full.names =TRUE, pattern = "rds")
files <- list.files("output/glmer", recursive=TRUE, full.names =TRUE, pattern = "glmer.results.rds")
out <- lapply(files, summarise.glmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/glmer.list.out.rds')
......
R/analysis/glmer.run.R
View file @ f9153d45
......@@ -56,7 +56,7 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.no.std.rds"))
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.global.rds"))
}
run.glmer <- function (model.file, trait, set, ecoregion,
......@@ -92,10 +92,10 @@ output.dir.glmer <- function (model, trait, set, ecoregion,type.filling) {
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait, set, ecoregion,
min.obs, sample.size, type.filling,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"), stringsAsFactors = FALSE)
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.global.csv"), stringsAsFactors = FALSE)
fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
}
......@@ -105,14 +105,14 @@ fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) )
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
# select species with minimum obs
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs])
# select obs abov min perc.neigh
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
return(data.tree)
}
......@@ -159,6 +159,7 @@ sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(dead=dead,
logD=logD,
species.id=species.id,
......@@ -176,12 +177,14 @@ return(data.frame(dead=dead,
# that will be used in a simple linear model
#============================================================
fun.data.for.glmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height")) stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height"))
stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
# get vars names
CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
CWM.tn <- paste(trait,"CWM",'fill',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".")
data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs)
#= DATA LIST FOR JAGS
......
R/analysis/lmer.global.std.nolog.output.R 0 → 100644
View file @ f9153d45
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
sets <- c('BCI','Canada','France','Fushan','NVS','Paracou',
'Spain','US','Swiss','Sweden','NSW','Mbaiki','Luquillo','Japan')
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (set in sets){
ecoregions <- get.ecoregions.for.set(set)
for (ecoregion in ecoregions){
for (trait in traits){
for (model in c(model.files.lmer.Tf.1,model.files.lmer.Tf.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model=model.obj$name, trait,
set, ecoregion,type.filling)
files <- c(files,file.path(pathout,"results.global.nolog.rds"))
}
}
}
}
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),
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='output/lmer.list.out.global.std.nolog.rds')
names.param <- unique(unlist(lapply(out,
function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out, fun.format.in.data.frame,
names.param = names.param))
DF.results$id <- paste(DF.results$set, DF.results$ecocode, sep=".")
saveRDS(DF.results,file='output/lmer.global.std.nolog.DF.rds')
R/analysis/lmer.nolog.all.output.R
View file @ f9153d45
......@@ -7,7 +7,7 @@ source("R/analysis/lmer.run.nolog.all.R")
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (trait in trait){
for (trait in traits){
for(model in c(model.files.lmer.Tf.1,model.files.lmer.Tf.2)){
source(model, local = TRUE)
model.obj <- load.model()
......
R/analysis/lmer.nolog.output.figs.R
View file @ f9153d45
This diff is collapsed.
R/analysis/lmer.output-fun.R
View file @ f9153d45
......@@ -191,7 +191,8 @@ 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){
mat.ratio <- DF.results[,c('sumTnTfBn.abs.VAR','sumTfBn.VAR','sumTnBn.VAR','effect.response.var')]/
mat.ratio <- DF.results[,c('sumTnTfBn.abs.VAR','sumTfBn.VAR','sumTnBn.VAR',
'effect.response.var')]/
DF.results[,'sumBn.VAR']
names(mat.ratio) <- c('abs.dist','Response','Effect','Effect.Response')
return(mat.ratio)
......@@ -199,19 +200,26 @@ return(mat.ratio)
### FUNCTION TO REPORT BEST MODEL PER ECOREGION AND TRAITS
fun.AIC <- function(id2.one,DF.results){
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','lmer.LOGLIN.AD.Tf')
best <- as.vector(DF.results[DF.results$id2==id2.one,c('id2','trait','set','ecocode','filling','MAT','MAP','model')])[which.min(DF.results$AIC[DF.results$id2==id2.one]),]
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',
'lmer.LOGLIN.AD.Tf')
best <- as.vector(DF.results[DF.results$id2==id2.one,c('id2','trait','set',
'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')])
AIC.all <- AIC.all[models]-min(AIC.all)
res <- data.frame((best),t(AIC.all))
names(res) <- c('id2','trait','set','ecocode','filling','MAT','MAP','best.model',paste('AIC',models,sep='.'))
names(res) <- c('id2','trait','set','ecocode','filling','MAT','MAP',
'best.model',paste('AIC',models,sep='.'))
return(res)
}
fun.AICc <- function(id2.one,DF.results){
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','lmer.LOGLIN.AD.Tf')
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','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']
Deviance.all <- Deviance.all[models]
......@@ -293,43 +301,104 @@ print(tapply(df.selected[['AIC.simplecomp']]<0,
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(param,"Std.Error",sep=".")]],
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*df.selected[[paste(param,"Std.Error",sep=".")]],
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*df.selected[[paste(param,"Std.Error",sep=".")]],
df.selected[[var.x]]+small.bar,df.selected[[param]] + 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
}
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,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected,threshold.delta.AIC=threshold.delta.AIC)
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
ylim <- range(c(df.selected[[params[1]]] - 1.96*df.selected[[paste(params[1],"Std.Error",sep=".")]],
df.selected[[params[1]]] + 1.96*df.selected[[paste(params[1],"Std.Error",sep=".")]]),na.rm=TRUE)
segments( df.selected[[var.x]],df.selected[[param]] - 1.96*df.selected[[paste(
param,"Std.Error",sep=".")]],
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*
df.selected[[paste(param,"Std.Error",sep=".")]],
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*
df.selected[[paste(param,"Std.Error",sep=".")]],
df.selected[[var.x]]+small.bar,df.selected[[param]] + 1.96*
df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
}
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,
ylim.same=FALSE,add.name = FALSE, ...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,
model.selected=model.selected,
threshold.delta.AIC=threshold.delta.AIC)
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{
col.vec2 <- 1}
if(!ylim.same) {ylim <- range(c(df.selected[[params[1]]] - 1.96*
df.selected[[paste(params[1],"Std.Error",sep=".")]],
df.selected[[params[1]]] + 1.96*
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,...)
fun.plot.param.error.bar(df.selected,var.x,param=params[1],small.bar,col=col.vec2)
}
fun.plot.all.param.boxplot <- function(model.selected,trait.name,DF.results,params,small.bar,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
if (add.name) text(df.selected[[var.x]], df.selected[[params[1]]],
labels = paste(df.selected[['set']],
df.selected[['ecocode']]), cex = 0.5, pos = 2)
fun.plot.param.error.bar(df.selected,var.x,param=params[1],
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(
param,"Std.Error",sep=".")]],
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*
df.selected[[paste(param,"Std.Error",sep=".")]],
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*
df.selected[[paste(param,"Std.Error",sep=".")]],
df.selected[[var.x]]+small.bar,df.selected[[param]] + df.selected[['sumBn']] + 1.96*
df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
}
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,
ylim.same=FALSE,add.name = FALSE, ...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,
model.selected=model.selected,
threshold.delta.AIC=threshold.delta.AIC)
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{
col.vec2 <- 1}
plot(df.selected[[var.x]],df.selected[[params[1]]] + df.selected[['sumBn']],
col = col.vec2, ...)
if (add.name) text(df.selected[[var.x]],
df.selected[[params[1]]] + df.selected[['sumBn']],
labels = paste(df.selected[['set']],
df.selected[['ecocode']]), cex = 0.5, pos = 2)
fun.plot.param.error.bar.std(df.selected,var.x,param=params[1],
small.bar,col=col.vec2)
}
fun.plot.all.param.boxplot <- function(model.selected,trait.name,DF.results,
params,small.bar,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,
model.selected=model.selected)
if(length(params)>1){
DF.t <- data.frame(param=rep(names(params),each=nrow(df.selected)),value=c(df.selected[[params[1]]],df.selected[[params[2]]]))
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']],...)
}else{
boxplot(df.selected[[params[1]]],...)
}
}
fun.plot.all.param.er.diff.MAP <- function(model.selected,trait.name,DF.results,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
fun.plot.all.param.er.diff.MAP <- function(model.selected,trait.name,
DF.results,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,
model.selected=model.selected)
plot(df.selected[['MAP']],df.selected[['sumTnBn']]-df.selected[['sumTfBn']],...)
}
fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y.l,express,...){
fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y.l,
express,...){
ncols = length(traits)
nrows = length(models)
list.models <- as.list(names(models))
......@@ -347,7 +416,8 @@ fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y.l,e
if(i==nrows)
mtext(var.x, side=1, line =4)
if(j==1)
mtext(paste('Effect size',list.models[i]), side=2, line =4,cex=0.9)
mtext(paste('Effect size',list.models[i]), side=2, line =4,
cex=0.9)
if(i==nrows & j==ncols)
legend('topright',legend=levels(DF.results$set),pch=16,
col=col.vec,bty='n',ncol=2)
......@@ -355,7 +425,8 @@ fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y.l,e
}
fun.plot.panel.lmer.res.boxplot <- function(models,traits,DF.results,var.y,express,...){
fun.plot.panel.lmer.res.boxplot <- function(models,traits,DF.results,var.y,
express,...){
ncols = length(traits)
nrows = length(models)
list.models <- as.list(names(models))
......@@ -369,12 +440,14 @@ fun.plot.panel.lmer.res.boxplot <- function(models,traits,DF.results,var.y,expre
if(i==1 )
mtext(traits[j], side=3, line =1)
if(j==1)
mtext(paste("Effect size", list.models[i],sep=" "), side=2, line =4,cex=0.9)
mtext(paste("Effect size", list.models[i],sep=" "), side=2,
line =4,cex=0.9)
}
}
fun.plot.panel.lmer.res.boxplot.compare <- function(models,traits,DF.results,var.y,express,...){
fun.plot.panel.lmer.res.boxplot.compare <- function(models,traits,DF.results,
var.y,express,...){
ncols = length(traits)
list.models <- as.list(names(models))
names(list.models) <- rep('model',length(list.models))
......@@ -389,7 +462,10 @@ fun.plot.panel.lmer.res.boxplot.compare <- function(models,traits,DF.results,var
}
}
fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,list.params,threshold.delta.AIC,small.bar=10,...){
fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,
list.params,threshold.delta.AIC,
small.bar=10,ylim=NA,
ylim.same=FALSE,...){
ncols = length(traits)
nrows = length(models)
par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
......@@ -400,10 +476,13 @@ fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,list
## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
for(i in 1:nrows)
for(j in 1:ncols){
fun.plot.all.param.x.y.c(models[i],traits[j],DF.results,var.x,params=list.params[[i]],
try(fun.plot.all.param.x.y.c(models[i],traits[j],DF.results,var.x,
params=list.params[[i]],
small.bar=small.bar,
threshold.delta.AIC=threshold.delta.AIC,col.vec=col.vec,col.set=TRUE,...)
abline(h=0,lty=3)
threshold.delta.AIC=threshold.delta.AIC,
col.vec=col.vec,col.set=TRUE,
ylim=ylim,ylim.same=ylim.same,...))
try(abline(h=0,lty=3))
if(length(list.params[[i]])>1)
legend("topright",names(list.params[[i]]),
pch=rep(1,length(list.params[[i]])),
......@@ -414,9 +493,9 @@ fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,list
mtext(var.x, side=1, line =4)
if(j==1)
mtext(paste('param',names(models)[i]), side=2, line =4,cex=1)
## if(i==nrows & j==ncols)
## legend('topright',legend=levels(DF.results$set),pch=16,
## col=col.vec,bty='n',ncol=2)
if(i==nrows & j==ncols)
legend('bottomright',legend=levels(DF.results$set),pch=16,
col=col.vec,bty='n',ncol=2)
}
}
......@@ -439,3 +518,25 @@ fun.plot.panel.lmer.parameters.boxplot <- function(models,traits,DF.results,list
}
}
## create a data base from the global data
fun.merge.results.global <- function(list.res,
names.param = c( "(Intercept)", "Tf",
"logD", "sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs",
"sumTnTfBn.diff")){
df.details <- data.frame(list.res$files.details[1:4],
list.res$lmer.summary[1:6])
fun.rep.df <- function(x, df = df.details) df
n_rep <- nrow(list.res$lmer.summary$ecocode.BLUP)
dat.t <- data.frame(matrix(rep(NA,n_rep * length(names.param)), nrow = n_rep,
ncol = length(names.param)))
names(dat.t) <- c(names.param)
dat.t[, match(names(list.res$lmer.summary$ecocode.BLUP), names(dat.t))] <-
list.res$lmer.summary$ecocode.BLUP
df.res <- data.frame(ecocode = rownames(list.res$lmer.summary$ecocode.BLUP),
do.call( "rbind", lapply(1:n_rep, fun.rep.df)),
dat.t)
return(df.res)
}
R/analysis/lmer.output.all-fun.R 0 → 100644
View file @ f9153d45
#### function to analyse lmer output for GLOBAL ANALYSIS
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){
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),
ecocode.BLUP=coef(x)$ecocode.i)
}
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),
lmer.summary=summarise.lmer.output( output.lmer))
}else{
res <- NULL
}
return(res)
}
files.details.all <- function(x){
s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
names(s) <- c("d1", "d2", "set", "trait", "filling", "model", "file" )
s[-(1:2)]
}
#### 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)
}
coef(.)$ecocode.i
ranef(.,whichel='ecocode.id',condVar=TRUE)
R/analysis/lmer.output.figs.R
View file @ f9153d45
......@@ -222,7 +222,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.02,
threshold.delta.AIC=10000)
threshold.delta.AIC=10000,ylim=c(-0.2,0.2),ylim.same=TRUE)
dev.off()
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf')
......@@ -258,7 +258,7 @@ pdf('figs/parameters.MAP.E.pdf',width=9,height=5)
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(-.15,.3),threshold.delta.AIC=1000)
list.params=list.params,small.bar=0.02,threshold.delta.AIC=1000)
dev.off()
......@@ -270,7 +270,7 @@ pdf('figs/parameters.MAP.R.pdf',width=9,height=5)
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(-.15,.3),threshold.delta.AIC=100000)
list.params=list.params,small.bar=0.02,threshold.delta.AIC=100000)
dev.off()
......@@ -284,7 +284,7 @@ pdf('figs/parameters.MAP.2models.pdf',width=12,height=10)
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(-.15,.3),threshold.delta.AIC=10000)
list.params=list.params,small.bar=0.02,threshold.delta.AIC=10000)
dev.off()