Skip to content

GitLab

Commit 025bdc23 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

changed to estimate with un std data plus advance on processing data

parent 01d2e58c
Hide whitespace changes
Inline Side-by-side
R/analysis.tar.gz
View file @ 025bdc23
No preview for this file type
R/analysis/glmer.output-fun.R
View file @ 025bdc23
......@@ -16,6 +16,7 @@ summarise.glmer.output <- function(x){
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))),
......@@ -109,7 +110,7 @@ dat.t[,match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
list.res$glmer.summary$fixed.coeff.Std.Error
dat.t[,match(names(list.res$glmer.summary$fixed.var),names(dat.t))] <-
list.res$glmer.summary$fixed.var
res <- data.frame(list.res$files.details,list.res$glmer.summary[1:6],dat.t)
res <- data.frame(list.res$files.details,list.res$glmer.summary[1:7],dat.t)
return(res)
}
......
R/analysis/glmer.run.R
View file @ 025bdc23
......@@ -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.rds"))
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.no.std.rds"))
}
run.glmer <- function (model.file, trait, set, ecoregion,
......@@ -95,7 +95,7 @@ load.and.prepare.data.for.glmer <- function(trait, set, ecoregion,
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.csv"), stringsAsFactors = FALSE)
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"), stringsAsFactors = FALSE)
fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
}
......
R/analysis/lmer.output-fun.R
View file @ 025bdc23
......@@ -16,6 +16,7 @@ summarise.lmer.output <- function(x){
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))),
......@@ -115,7 +116,7 @@ dat.t[,length(names.param)+match(names(list.res$lmer.summary$fixed.coeff.E),name
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:6],dat.t)
res <- data.frame(list.res$files.details,list.res$lmer.summary[1:7],dat.t)
return(res)
}
......@@ -211,7 +212,7 @@ fun.AICc <- function(id2.one,DF.results){
### function to get the data for a given model with criteria to select
fun.select.ecoregions.trait <- function(DF.results,trait.name,model.selected,
nobs.min=1000,filling.selected="species",
threshold.delta.AIC){ ## NEED TO TEST WITH DIFFERENT TRESHOLD HERE
threshold.delta.AIC){
DF.results[DF.results$nobs>nobs.min &
DF.results$filling==filling.selected &
DF.results$trait==trait.name &
......@@ -222,7 +223,7 @@ DF.results[DF.results$nobs>nobs.min &
### function to get the data for a given model with criteria to select only site with competition
fun.select.ecoregions.trait.compet <- function(DF.results,trait.name,model.selected,
nobs.min=1000,filling.selected="species",
threshold.delta.AIC){ ## NEED TO TEST WITH DIFFERENT TRESHOLD HERE
threshold.delta.AIC){
DF.results[DF.results$nobs>nobs.min &
DF.results$filling==filling.selected &
DF.results$trait==trait.name &
......@@ -232,35 +233,21 @@ DF.results[DF.results$nobs>nobs.min &
DF.results$delta.AIC<threshold.delta.AIC,]
}
### function to get the data for a given model with criteria to select only site with facilitation
fun.select.ecoregions.trait.facilit <- function(DF.results,trait.name,model.selected,
nobs.min=1000,filling.selected="fill",
threshold.delta.AIC){ ## NEED TO TEST WITH DIFFERENT TRESHOLD HERE
perc <- "perc.genus"
sp <- "CWM.species"
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 &
## DF.results$delta.AIC==0,]
DF.results$delta.AIC<threshold.delta.AIC,]
}
#########################
##### FUNCTIONS FOR PLOTS
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)
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]],...)
}
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=100000)
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}
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec <- 1}
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
plot(df.selected[[var.x]],df.selected[[var.y]],
pch=pch.vec,
cex=cex.vec,
......@@ -286,43 +273,23 @@ print(tapply(df.selected[['AIC.simplecomp']]<0,
fun.plot.param.error.bar <- function(df.selected,var.x,param,small.bar,...){
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=".")]],...)
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=".")]])
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=".")]])
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,...){
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)
plot(df.selected[[var.x]],df.selected[[params[1]]],...)
fun.plot.param.error.bar(df.selected,var.x,param=params[1],small.bar)
if(length(params)>1){
col <- 2
for (i in params[-1]){
points(df.selected[[var.x]],df.selected[[i]],col=col)
fun.plot.param.error.bar(df.selected,var.x,param=i,small.bar,col=col)
col <- col+1
}
}
}
fun.plot.all.param.x.y.f <- function(model.selected,trait.name,DF.results,var.x,params,small.bar,...){
df.selected <- fun.select.ecoregions.trait.facilit(DF.results,trait.name=trait.name,model.selected=model.selected)
plot(df.selected[[var.x]],df.selected[[params[1]]],...)
fun.plot.param.error.bar(df.selected,var.x,param=params[1],small.bar)
if(length(params)>1){
col <- 2
for (i in params[-1]){
points(df.selected[[var.x]],df.selected[[i]],col=col)
fun.plot.param.error.bar(df.selected,var.x,param=i,small.bar,col=col)
col <- col+1
}
}
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
plot(df.selected[[var.x]],df.selected[[params[1]]],col=col.vec2,...)
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(length(params)>1){
......@@ -340,7 +307,7 @@ 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,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))
......@@ -351,7 +318,7 @@ fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y,exp
for(i in 1:nrows)
for(j in 1:ncols){
fun.plot.lmer.res.x.y.2(models[i],traits[j],
DF.results,var.x,var.y,col.vec,...)
DF.results,var.x,var.y=var.y.l[[i]],col.vec,...)
abline(h=0,lty=3)
if(i==1 )
mtext(traits[j], side=3, line =1)
......@@ -404,35 +371,16 @@ fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,list
ncols = length(traits)
nrows = length(models)
par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
## ### 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){
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,...)
abline(h=0,lty=3)
if(length(list.params[[i]])>1)
legend("topright",names(list.params[[i]]),
pch=rep(1,length(list.params[[i]])),
col=1:length(list.params[[i]]),bty='n',cex=1)
if(i==1 )
mtext(traits[j], side=3, line =1)
if(i==nrows)
mtext(var.x, side=1, line =4)
if(j==1)
mtext(paste('param',names(models)[i]), side=2, line =4,cex=1)
}
}
DF.results$set <- factor(DF.results$set)
col.vec <- niceColors(n=nlevels(DF.results$set))
fun.plot.panel.lmer.parameters.f <- function(models,traits,DF.results,var.x,list.params,small.bar=10,...){
ncols = length(traits)
nrows = length(models)
par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
## ### TO COMPARE THE PARAMTERS WE NEED TO DIVIDE THE ABS>DIST per two
## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
## 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.f(models[i],traits[j],DF.results,var.x,params=list.params[[i]],small.bar=small.bar,...)
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)
if(length(list.params[[i]])>1)
legend("topright",names(list.params[[i]]),
......@@ -444,6 +392,9 @@ fun.plot.panel.lmer.parameters.f <- 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)
}
}
......
R/analysis/lmer.output.figs.R
View file @ 025bdc23
......@@ -42,8 +42,6 @@ DF.AIC.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.di
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)
......@@ -59,7 +57,18 @@ names(DF.AIC.weights) <- c('id2',paste('AIC.weight',names(DF.AIC.weights)[-1],se
DF.best.and.all.AIC <- merge(DF.best.and.all.AIC,DF.AIC.weights,by='id2')
## TODO NOT WORKING need to print AIC table in markdown format
#### 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
library('pander')
fun.AIC.print.pandoc.table.trait <- function(DF.best.and.all.AIC,trait.select){
DF.best.and.all.AIC$best.model <- factor(DF.best.and.all.AIC$best.model)
......@@ -81,124 +90,179 @@ fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select=c("SLA","Leaf.
#############################################
#############################################
### DO THE PLOT
models <- c('lmer.LOGLIN.AD.Tf','lmer.LOGLIN.ER.Tf')
names(models) <- c('Absolute distance','Effect/response')
pdf('figs/R2.boxplot.two.pdf',width=16,height=5)
fun.plot.panel.lmer.res.boxplot.compare(models=models,
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,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
var.y='R2m.ratio',ylim=c(-0.015,1.5),
col=c(rgb(87, 157, 25,maxColorValue=255),rgb(255,102,51,maxColorValue=255)),
cex.axis=1.7)
var.x='MAP',var.y='R2m.simplecomp',ylim=c(-0.015,0.06))
dev.off()
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()
models <- c('lmer.LOGLIN.AD')
names(models) <- c('Absolute distance')
pdf('figs/R2.boxplot.one.pdf',width=8,height=3)
fun.plot.panel.lmer.res.boxplot.compare(models=models,
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,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
var.y='R2c.simplecomp',ylim=c(-0.015,0.06),
col=c(rgb(87, 157, 25,maxColorValue=255)), cex.axis=1.7)
var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,100))
dev.off()
models <- c('lmer.LOGLIN.ER')
names(models) <- c('Effect/response')
pdf('figs/R2.MAP.pdf',width=8,height=4)
models <- c('lmer.LOGLIN.AD.Tf')
names(models) <- c()
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='R2c.simplecomp',ylim=c(-0.015,0.06))
dev.off()
var.x='MAP',var.y='sumTnTfBn.abs.VAR',ylim=c(-0.015,0.2))
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
pdf('figs/R2.MAP.two.pdf',width=10,height=7)
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')
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='R2m.simplecomp',ylim=c(-0.015,0.08))
dev.off()
var.x='MAP',var.y='ER.perc.var',ylim=c(-0.015,100))
models <- c('lmer.LOGLIN.AD.Tf')
names(models) <- c('AB')
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='R2m.nocomp',ylim=c(-0.015,0.1))
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))
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')
pdf('figs/R2.boxplot.pdf',width=8,height=6)
fun.plot.panel.lmer.res.boxplot(models=models,
pdf('figs/R2.MAP.two.pdf',width=10,height=7)
fun.plot.panel.lmer.res.x.y(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
var.y='R2c.simplecomp',ylim=c(-0.015,0.06))
var.x='MAP',var.y='R2m.simplecomp',ylim=c(-0.015,0.08))
dev.off()
### plot parameters
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.E.Tf','lmer.LOGLIN.R.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Effect','Response','Absolute distance')
list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'),
c('sumTnBn'),
c('sumTfBn'),
c('sumTnTfBn.abs'))
pdf('figs/parameters.MAP.4models.pdf',width=9,height=9)
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)
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(-.1,.2))
list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=10000)
dev.off()
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'),
c('sumTnTfBn.abs'))
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.2models.pdf',width=9,height=5)
pdf('figs/parameters.MAP.ER.best.aic.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.02,ylim=c(-.15,.2),threshold.delta.AIC=10)
list.params=list.params,small.bar=0.02,ylim=c(-.15,.25),threshold.delta.AIC=1)
dev.off()
models <- c('lmer.LOGLIN.E.Tf')
names(models) <- c('Effect effect')
list.params <- list(c(Response='sumTnBn'))
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('BATOT')
list.params <- list(c(Bn='sumBn'))
pdf('figs/parameters.BATOT.MAP.pdf',width=8,height=4)
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(-0.25,0.04))
list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=1000)
dev.off()
models <- c('lmer.LOGLIN.ER')
names(models) <- c('Effect/response')
list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'))
models <- c('lmer.LOGLIN.R.Tf')
names(models) <- c('Response response')
list.params <- list(c(Response='sumTfBn'))
pdf('figs/parameters.boxplot.ER.pdf',width=8,height=3)
fun.plot.panel.lmer.parameters.boxplot(models=models,
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(-0.1,0.13),col=c('red','black'))
list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=100000)
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,
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,
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))
list.params=list.params,small.bar=0.02,ylim=c(-.15,.3),threshold.delta.AIC=10000)
dev.off()
## 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()
R/analysis/lmer.run.R
View file @ 025bdc23
......@@ -56,7 +56,7 @@ fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
end <- Sys.time()
print(end - Start)
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.rds"))
saveRDS(lmer.output,file=file.path(path.out, "results.no.std.rds"))
}
run.lmer <- function (model.file, trait, set, ecoregion,
......@@ -100,7 +100,7 @@ load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
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.csv"), stringsAsFactors = FALSE)
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"), stringsAsFactors = FALSE)
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
......
R/analysis/plots.resid.regression.lines.R
View file @ 025bdc23
......@@ -14,15 +14,26 @@ fun.compute.resid <- function(i){
return(fitted(i) +resid(i) -i@pp$X %*%fixef(i))
}
fun.plot.colors.density <- function(x,y,...){
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
plot(x,y, col=colors.dens, pch=20,...)
}
fun.points.colors.density <- function(x,y,...){
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
points(x,y, col=colors.dens, pch=20,...)
}
fun.boxplot.breaks <- function(x,y,Nclass=15,...){
breakss <- seq(min(x),max(x),length=Nclass+1)
x.cut <- cut(x,breakss)
mid.points <- breakss[-(Nclass+1)]+(breakss[2]-breakss[1])/2
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
plot(x,y, col=colors.dens, pch=20,...)
boxplot(y~x.cut,at=mid.points,add=TRUE,names=NA)
boxplot(y~x.cut,at=mid.points,add=TRUE,names=NA,outline=FALSE)
}
seq.from.to.quantile <- function(x,length.out=20,probs=c(0.001,0.999)){
......@@ -30,71 +41,160 @@ qq <- quantile(x,probs)
return(seq(from=qq[1],to=qq[2],length.out=length.out))
}
fun.plot.residual.plus.regression.lines <- function(df.lmer,res.fix.no,res.fix.simple,
fun.plot.residual.plus.regression.lines.all <- function(df.lmer,res.fix.no,res.fix.simple,
ERcomp,ADcomp,trait,set,ecoregion,
type.filling){
par(mfrow=c(2,3),oma=c(0,0,2,0))
## Effect /reponse
fun.boxplot.breaks(df.lmer$sumBn,res.fix.no,cex=0.1,
fun.plot.colors.density(df.lmer$sumBn,res.fix.no,cex=0.1,
xlab="sum of basal area",
ylab="growth residual",
main="Effect/response model")
fun.boxplot.breaks(df.lmer$sumBn,res.fix.no)
mtext(paste(trait,set,ecoregion,type.filling), side=3,line=0.1,outer=TRU