Commit 92b3ce98 authored by ghislainv's avatar ghislainv
Browse files

resolve conflict Paracou

parents c560d2f8 db81af37
......@@ -4,7 +4,7 @@ D3 := output/processed
D4 := figs/test.format.tree
D5 := figs/test.traits
D5 := figs/test.CWM
sites:= Fushan Paracou BCI Mbaiki Luquillo Canada France Spain US Sweden Swiss NSW NVS Japan
sites:= Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS US
D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) ))
D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) ))
D2tree := $(addsuffix /tree.csv,$(addprefix $(D2)/, $(sites) ))
......@@ -26,10 +26,10 @@ $(D2)/TRY/data.TRY.std.rds:
BCI: $(D3)/BCI/Done.txt
$(D3)/BCI/Done.txt: R/process.data/process.fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15); process_bigplot_dataset('BCI', Rlim=15,std.traits=FALSE);"
$(D3)/BCI/Done.txt: R/process.data/process-fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local'); process_bigplot_dataset('BCI', Rlim=15,std.traits='global'); process_bigplot_dataset('BCI', Rlim=15,std.traits='no');"
$(D2)/BCI/traits.csv: R/find.trait/BCI.R R/find.trait/trait.fun.R $(D2)/BCI/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/BCI/traits.csv: R/find.trait/BCI.R R/find.trait/trait-fun.R $(D2)/BCI/tree.csv
Rscript $<
$(D2)/BCI/tree.csv: R/format.data/BCI.R $(shell find $(D1)/BCI -type f)
......@@ -39,10 +39,10 @@ $(D2)/BCI/tree.csv: R/format.data/BCI.R $(shell find $(D1)/BCI -type f)
Japan: $(D3)/Japan/Done.txt
$(D3)/Japan/Done.txt: R/process.data/process.fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15); process_bigplot_dataset('Japan', Rlim=15,std.traits=FALSE);"
$(D3)/Japan/Done.txt: R/process.data/process-fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local'); process_bigplot_dataset('Japan', Rlim=15,std.traits='global'); process_bigplot_dataset('Japan', Rlim=15,std.traits='no');"
$(D2)/Japan/traits.csv: R/find.trait/Japan.R R/find.trait/trait.fun.R $(D2)/Japan/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Japan/traits.csv: R/find.trait/Japan.R R/find.trait/trait-fun.R $(D2)/Japan/tree.csv
Rscript $<
$(D2)/Japan/tree.csv: R/format.data/Japan.R $(shell find $(D1)/Japan -type f)
......@@ -51,10 +51,10 @@ $(D2)/Japan/tree.csv: R/format.data/Japan.R $(shell find $(D1)/Japan -type f)
Luquillo: $(D3)/Luquillo/Done.txt
$(D3)/Luquillo/Done.txt: R/process.data/process.fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15); process_bigplot_dataset('Luquillo', Rlim=15,std.traits=FALSE);"
$(D3)/Luquillo/Done.txt: R/process.data/process-fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='no');process_bigplot_dataset('Luquillo', Rlim=15,std.traits='global');"
$(D2)/Luquillo/traits.csv: R/find.trait/Luquillo.R R/find.trait/trait.fun.R $(D2)/Luquillo/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Luquillo/traits.csv: R/find.trait/Luquillo.R R/find.trait/trait-fun.R $(D2)/Luquillo/tree.csv
Rscript $<
$(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -type f)
......@@ -64,10 +64,10 @@ $(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -t
Mbaiki: $(D3)/Mbaiki/Done.txt
$(D3)/Mbaiki/Done.txt: R/process.data/process.fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits=FALSE);"
$(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='local'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='no');"
$(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait.fun.R $(D2)/Mbaiki/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait-fun.R $(D2)/Mbaiki/tree.csv
Rscript $<
$(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f)
......@@ -77,10 +77,10 @@ $(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f)
Canada: $(D3)/Canada/Done.txt
$(D3)/Canada/Done.txt: R/process.data/process.fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada'); process_inventory_dataset('Canada',std.traits=FALSE);"
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local'); process_inventory_dataset('Canada',std.traits='global');process_inventory_dataset('Canada',std.traits='no');"
$(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait.fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait-fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f)
......@@ -89,10 +89,10 @@ $(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f)
#-------------------------------------------------------
France: $(D3)/France/Done.txt
$(D3)/France/Done.txt: R/process.data/process.fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France'); process_inventory_dataset('France',std.traits=FALSE);"
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local'); process_inventory_dataset('France',std.traits='global');process_inventory_dataset('France',std.traits='no');"
$(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait.fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait-fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f)
......@@ -102,11 +102,11 @@ $(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f)
Fushan: $(D3)/Fushan/Done.txt
$(D3)/Fushan/Done.txt: R/process.data/process.fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15); process_bigplot_dataset('Fushan', Rlim=15,std.traits=FALSE);"
$(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='local'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');process_bigplot_dataset('Fushan', Rlim=15,std.traits='no');"
$(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait.fun.R $(D2)/Fushan/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait-fun.R $(D2)/Fushan/tree.csv
Rscript $<
$(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
......@@ -116,10 +116,10 @@ $(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
NSW: $(D3)/NSW/Done.txt
$(D3)/NSW/Done.txt: R/process.data/process.fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW'); process_inventory_dataset('NSW',std.traits=FALSE);"
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local'); process_inventory_dataset('NSW',std.traits='global');process_inventory_dataset('NSW',std.traits='no');"
$(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait.fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait-fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f)
......@@ -129,10 +129,10 @@ $(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f)
NVS: $(D3)/NVS/Done.txt
$(D3)/NVS/Done.txt: R/process.data/process.fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS'); process_inventory_dataset('NVS',std.traits=FALSE);"
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local'); process_inventory_dataset('NVS',std.traits='global'); process_inventory_dataset('NVS',std.traits='no');"
$(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait.fun.R $(D2)/NVS/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait-fun.R $(D2)/NVS/tree.csv
Rscript $<
$(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f)
......@@ -142,10 +142,10 @@ $(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f)
Paracou: $(D3)/Paracou/Done.txt
$(D3)/Paracou/Done.txt: R/process.data/process.fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15); process_bigplot_dataset('Paracou', Rlim=15,std.traits=FALSE);"
$(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='local'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='no');"
$(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait.fun.R $(D2)/Paracou/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait-fun.R $(D2)/Paracou/tree.csv
Rscript $<
$(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type f)
......@@ -155,10 +155,10 @@ $(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type
Spain: $(D3)/Spain/Done.txt
$(D3)/Spain/Done.txt: R/process.data/process.fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain'); process_inventory_dataset('Spain',std.traits=FALSE);"
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local'); process_inventory_dataset('Spain',std.traits='global'); process_inventory_dataset('Spain',std.traits='no');"
$(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait.fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait-fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f)
......@@ -168,10 +168,10 @@ $(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f)
Sweden: $(D3)/Sweden/Done.txt
$(D3)/Sweden/Done.txt: R/process.data/process.fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden'); process_inventory_dataset('Sweden',std.traits=FALSE);"
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local'); process_inventory_dataset('Sweden',std.traits='global'); process_inventory_dataset('Sweden',std.traits='no');"
$(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait.fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait-fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f)
......@@ -181,10 +181,10 @@ $(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f)
Swiss: $(D3)/Swiss/Done.txt
$(D3)/Swiss/Done.txt: R/process.data/process.fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss'); process_inventory_dataset('Swiss',std.traits=FALSE);"
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local'); process_inventory_dataset('Swiss',std.traits='global'); process_inventory_dataset('Swiss',std.traits='no');"
$(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait.fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait-fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f)
......@@ -194,10 +194,10 @@ $(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f)
US: $(D3)/US/Done.txt
$(D3)/US/Done.txt: R/process.data/process.fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US'); process_inventory_dataset('US',std.traits=FALSE);"
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local'); process_inventory_dataset('US',std.traits='global'); process_inventory_dataset('US',std.traits='no');"
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait.fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/US/tree.csv: R/format.data/US.R $(shell find $(D1)/US -type f)
......@@ -223,6 +223,15 @@ $(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D5)/Done.txt $(D3Done)
#-------------------------------------------------------
#-------------------------------------------------------
MERGE.ALL: output/processed/data.all.csv
output/processed/data.all.csv: R/process.data/merge.all.processed.data.R
Rscript $<
#-------------------------------------------------------
# This susbtitution rule should work as rule, but not, why not?
# docs/output/formatted/%/tree.csv: $(D1)/%/* %.R
# Rscript %.R
......
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.R")
## TODO NEED TO CHANGE THAT TO 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.nolog.rds"))
}
}
}
}
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/lmer.list.out.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.nolog.DF.rds')
#!/usr/bin/env Rscript
rm(list=ls())
source("R/analysis/lmer.output-fun.R")
source("R/utils/plot.R")
## load results
DF.results <- readRDS('output/lmer.nolog.DF.rds')
## 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=".")
## par(mfrow=c(1,2),mar=c(5, 9, 4, 1) +0.1)
## boxplot(site.clim.all$clim.all.MAT~site.clim.all$id,horizontal=TRUE,las=2,cex.axis=0.8)
## boxplot(site.clim.all$clim.all.MAP~site.clim.all$id,horizontal=TRUE,las=2,cex.axis=0.8)
mean.MAT <- tapply(site.clim.all$clim.all.MAT,INDEX=site.clim.all$id,FUN=mean,na.rm=TRUE)
mean.MAP <- tapply(site.clim.all$clim.all.MAP,INDEX=site.clim.all$id,FUN=mean,na.rm=TRUE)
data.clim.ecocode <- data.frame(id=names(mean.MAT),MAT=mean.MAT,MAP=mean.MAP)
### merge climate and lmer results
DF.results <- merge(DF.results,data.clim.ecocode,by="id")
DF.results$id2 <- paste(DF.results$id,DF.results$trait,DF.results$filling)
### save
saveRDS(DF.results,file='output/lmer.nolog.DF.results.merged.rds')
DF.results <- readRDS(file='output/lmer.nolog.DF.results.merged.rds')
### Analysis of the results
DF.R2m.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"R2m"))
DF.R2c.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"R2c"))
DF.AIC.diff <- do.call("rbind",lapply(1:nrow(DF.results),fun.compute.criteria.diff,DF.results,"AIC"))
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)
### report best model based on 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)
t(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',]$set))/(apply(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',]$set),
MARGIN=2,sum))
## 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: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 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
#############################################
#############################################
### DO THE PLOT
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,threshold.delta.AIC=10000)
dev.off()
......@@ -283,10 +283,13 @@ segments( df.selected[[var.x]]-small.bar,df.selected[[param]] + 1.96*df.selected
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,...){
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}
plot(df.selected[[var.x]],df.selected[[params[1]]],col=col.vec2,...)
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)
}
......@@ -393,9 +396,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('topright',legend=levels(DF.results$set),pch=16,
## col=col.vec,bty='n',ncol=2)
}
}
......
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
files <- list.files("output/lmer", recursive=TRUE, full.names =TRUE, pattern = "rds")
source("R/analysis/lmer.run.R")
## TODO NEED TO CHANGE THAT TO 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.no.std.rds"))
}
}
}
}
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/lmer.list.out.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))
......
......@@ -49,10 +49,18 @@ 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)
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)
t(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',]$set))/(apply(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',]$set),
MARGIN=2,sum))
## AIC weights
AIC.weights <- do.call('rbind',
......@@ -63,6 +71,20 @@ 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')
## assume ER best model if either E or R best model
f <- function(i,DF){
d.AIC <- as.vector(DF$delta.AIC[DF$id2==i] )
best <- DF$model[DF$id2==i][DF$delta.AIC[DF$id2==i]==0]
if(best %in% c('lmer.LOGLIN.E','lmer.LOGLIN.R')) {
d.AIC[ DF$model[DF$id2==i]=='lmer.LOGLIN.ER'] <- 0}
return(d.AIC)
}
d.AIC <- do.call('c',lapply(unique(DF.results$id2),FUN=f,DF.results))
DF.results$delta.AIC <- d.AIC
## 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
......@@ -136,14 +158,14 @@ fun.plot.panel.lmer.res.x.y(models=models,
var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.17))
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.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.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.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()
......@@ -199,7 +221,8 @@ 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(-.15,.25),threshold.delta.AIC=10000)
list.params=list.params,small.bar=0.02,
threshold.delta.AIC=10000)
dev.off()
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.ER.Tf')
......
......@@ -50,13 +50,16 @@ fun.test.if.multi.census <- function(data){
return("tree.id" %in% names(data))
}
fun.call.lmer <- function(formula,df.lmer){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
return(lmer.output)
}
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
end <- Sys.time()
print(end - Start)
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.no.std.rds"))
saveRDS(lmer.output,file=file.path(path.out, "results.rds"))
}
run.lmer <- function (model.file, trait, set, ecoregion,
......@@ -82,10 +85,37 @@ run.lmer <- function (model.file, trait, set, ecoregion,
#= Run model
if(test.if.multi.census){
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out)
fun.ci.boot(df.lmer,formula=model$lmer.formula.tree.id,path.out,level=0.95,nsim=500)
}else{
fun.call.lmer.and.save(formula=model$lmer.formula,df.lmer,path.out)
fun.ci.boot(df.lmer,formula=model$lmer.formula,path.out,level=0.95,nsim=500)
}
}
## new function to compute boot ci
fun.ci.boot <- function(df.lmer,formula,path.out,level=0.95,nsim=500){
require(boot)
require(multicore)
bb <- boot(data=df.lmer, statistic=boot.fun,R= nsim,formula=formula)
bci <- lapply(seq_along(bb$t0), boot.out = bb, boot::boot.ci,
type = "perc", conf = level)
citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- paste('CI',round(a, 3),sep='.')
dimnames(citab) <- list(names(bb[["t0"]]), pct)
saveRDS(citab,file=file.path(path.out, "results.ci.rds"))
}
boot.fun <- function(data, indices, formula){
df.lmer <- data[indices,] # select obs. in bootstrap sample
res <- fun.call.lmer(formula=formula,df.lmer)
fixef(res)
}
#========================================================================
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
......@@ -98,7 +128,7 @@ 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/",std=TRUE){
base.dir = "output/processed/",std=FALSE){
### load data
if(std) { data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"),
stringsAsFactors = FALSE)}else{
......
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA
library(lme4)
get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
}
run.models.for.set.all.traits <- function(set,model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),type.filling, std,...){
for(trait in traits)
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, set, type.filling=type.filling, std=std,...)
}
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), std, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling,std=std, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), std, ...){
fun.model <- match.fun(fun.model)
for (e in ecoregions)
try(fun.model(model.file, trait, set, e, type.filling=type.filling,std=std,...))
}
#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.lmer.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.R.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.R")
model.files.lmer.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R")
model.files.lmer.Tf.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.R.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R")
model.files.lmer.Tf.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.HD.Tf.R",