Skip to content

GitLab

Commit 67ee0df5 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

minor change in function for output

parent 025bdc23
Hide whitespace changes
Inline Side-by-side
R/analysis/lmer.output-fun.R
View file @ 67ee0df5
......@@ -181,7 +181,8 @@ 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.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')
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')])
......@@ -192,14 +193,14 @@ fun.AIC <- function(id2.one,DF.results){
}
fun.AICc <- function(id2.one,DF.results){
models <- c('lmer.LOGLIN.nocomp.Tf', 'lmer.LOGLIN.simplecomp.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]
nobs.all <- DF.results[DF.results$id2==id2.one,'nobs']
names(nobs.all) <- DF.results[DF.results$id2==id2.one,'model']
nobs.all <- nobs.all[models]
n.param <- c(2,3,4,4,5,4)
n.param <- c(2,3,4,4,4,5,4)
AICc <- Deviance.all+2*n.param*(nobs.all)/(nobs.all-n.param-1)
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)
......
R/analysis/lmer.output.figs.R
View file @ 67ee0df5
......@@ -14,7 +14,6 @@ DF.results <- do.call("rbind",lapply(list.lmer,fun.format.in.data.frame,names.pa
# add id
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
## load climatic data
site.clim.all <- read.csv( file.path("output/processed", "all.sites.clim.csv"), stringsAsFactors = FALSE)
site.clim.all$id <- paste(site.clim.all$set,site.clim.all$ecocode,sep=".")
......@@ -50,24 +49,51 @@ DF.results <- cbind(DF.results,DF.R2m.diff,DF.R2c.diff,DF.AIC.diff,DF.delta.AIC,
DF.best.and.all.AIC <- do.call('rbind',lapply(unique(DF.results$id2),FUN=fun.AIC,DF.results))
DF.best.and.all.AICc <- do.call('rbind',lapply(unique(DF.results$id2),FUN=fun.AICc,DF.results))
table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$trait,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set)
## AIC weights
AIC.weights <- do.call('rbind',lapply(1:nrow(DF.best.and.all.AICc),FUN=function(i,DF) exp((min(DF[i,])-DF[i,])/2)/sum(exp((min(DF[i,])-DF[i,])/2)),DF.best.and.all.AIC[,9:14]))
AIC.weights <- do.call('rbind',
lapply(1:nrow(DF.best.and.all.AICc),
FUN=function(i,DF) exp((min(DF[i,])-DF[i,])/2)/sum(exp((min(DF[i,])-DF[i,])/2)),
DF.best.and.all.AIC[,9:15]))
DF.AIC.weights <- data.frame(DF.best.and.all.AICc[,1],AIC.weights)
names(DF.AIC.weights) <- c('id2',paste('AIC.weight',names(DF.AIC.weights)[-1],sep='.'))
DF.best.and.all.AIC <- merge(DF.best.and.all.AIC,DF.AIC.weights,by='id2')
#### compute percentage of vqariance explained by var
## compute a global AIC
fun.global.aic <- function(DF.results){
DF.results <- DF.results[DF.results$nobs>1000,] # select set ecocode with more than 1000 obs
DF.results <- DF.results[DF.results$id2 %in% names(table(DF.results$id2))[table(DF.results$id2)==7],] # select set ecocode with 7 model tested
## species
DF.results.sp <- DF.results[DF.results$filling=='species',] # select set ecocode with more than 1000 obs
test.same.n.model.ecoregion <- apply(table(DF.results.sp$trait,DF.results.sp$model),MARGIN=1,function(x) all(x==x[1]))
if(!all(test.same.n.model.ecoregion))
stop(paste('error not the same number of ecoregion for traits',
names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
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,DF.results.sp$model),FUN=length))
## genus
DF.results.ge <- DF.results[DF.results$filling=='genus',] # select set ecocode with more than 1000 obs
test.same.n.model.ecoregion <- apply(table(DF.results.ge$trait,DF.results.ge$model),MARGIN=1,function(x) all(x==x[1]))
if(!all(test.same.n.model.ecoregion))
stop(paste('error not the same number of ecoregion for traits',
names(test.same.n.model.ecoregion)[!test.same.n.model.ecoregion]))
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),FUN=length))
return(list(sp=list.sp,ge=list.ge))
}
global.AIC.list <- fun.global.aic(DF.results)
#### compute percentage of variance 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){
......@@ -76,8 +102,8 @@ aa <- as.data.frame.matrix( t(with(subset(DF.best.and.all.AIC,subset=trait %in%
table(best.model,set))))
AIC.pandoc.table <- data.frame(set=rownames(aa),aa[,c(4,6,1,5,2,3)],row.names=NULL)
names(AIC.pandoc.table) <- gsub("lmer.LOGLIN.","", names(AIC.pandoc.table))
print(pandoc.table(AIC.pandoc.table,caption=paste("Best models per data set for trait",trait.select)))
}
print(pandoc.table(AIC.pandoc.table,caption=paste("Best models per data set for trait",trait.select)))}
fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="SLA")
fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select="Leaf.N")
......@@ -91,13 +117,14 @@ fun.AIC.print.pandoc.table.trait(DF.best.and.all.AIC,trait.select=c("SLA","Leaf.
#############################################
### DO THE PLOT
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
var.y.l <- list('R2m.simplecomp','R2m.simplecomp')
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.x='MAP',var.y='R2m.simplecomp',ylim=c(-0.015,0.06))
var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.06))
dev.off()
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
......@@ -112,44 +139,29 @@ dev.off()
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)
pdf('figs/perc.var.relative.BATOT.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,100))
dev.off()
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='sumTnTfBn.abs.VAR',ylim=c(-0.015,0.2))
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='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='abs.perc.var',ylim=c(-0.015,100))
par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(2,2))
boxplot(effect.response.var~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
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)
boxplot(sumTfBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
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)
boxplot(sumTnBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
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)
boxplot(sumTnTfBn.abs.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(2,2))
boxplot(sumBn.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
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)
boxplot(Tf.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
par(mar=c(5.1,9.1,4.1,2.1))
boxplot(logD.VAR~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3,ylim=c(0,0.5))
par(mar=c(5.1,9.1,4.1,2.1),mfrow=c(1,2))
......@@ -157,18 +169,19 @@ boxplot(ER.perc.var~trait+set,data=DF.results[DF.results$nobs>1000,],horizontal=
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(R2m~trait+model+set,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
boxplot(R2m~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,3))
boxplot(sumTnBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
abline(v=0,col='red')
boxplot(sumTfBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
abline(v=0,col='red')
boxplot(sumBn~set+trait,data=DF.results[DF.results$nobs>1000,],horizontal=TRUE,las=2,cex=0.3)
abline(v=0,col='red')
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)
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()
......@@ -189,6 +202,18 @@ fun.plot.panel.lmer.parameters.c(models=models,
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.ER.Tf')
names(models) <- c('Effect/response effect','Effect/response response')
list.params <- list(c(Response='sumTnBn'),
c(Effect='sumTfBn'))
pdf('figs/parameters.MAT.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='MAT',
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.ER.Tf')
names(models) <- c('Effect/response effect','Effect/response response')
......@@ -241,16 +266,16 @@ dev.off()
## models <- c('lmer.LOGLIN.ER')
## names(models) <- c('Effect/response')
## list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'))
models <- c('lmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response')
list.params <- list('sumBn')
## 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()
pdf('figs/parameters.MAP.sumBn.ER.pdf',width=8,height=3)
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.4,0.13),threshold.delta.AIC=10000)
dev.off()
## models <- c('lmer.LOGLIN.AD')
## names(models) <- c('Absolute distance')
......
R/analysis/lmer.run.R
View file @ 67ee0df5
......@@ -98,9 +98,13 @@ output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
#============================================================
load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
base.dir = "output/processed/",std=TRUE){
### load data
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"), stringsAsFactors = FALSE)
if(std) { data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"),
stringsAsFactors = FALSE)}else{
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.csv"),
stringsAsFactors = FALSE)}
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
......
R/analysis/plots.resid.regression.lines.R
View file @ 67ee0df5
......@@ -148,7 +148,7 @@ lines(seq.from.to.quantile(df.lmer$sumTfBn),
fun.load.data.for.residual <- function(trait,set,ecoregion,type.filling){
df.lmer <- load.and.prepare.data.for.lmer(trait,set,ecoregion,
min.obs=10, sample.size=NA,
type.filling) # return a DF
type.filling,std=FALSE) # return a DF
simple <- readRDS(file.path("output/lmer", set,ecoregion,
trait,type.filling,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment