diff --git a/R/analysis.tar.gz b/R/analysis.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..77bf3dabd28d0d9ede0c9c015b591c20137d630a
Binary files /dev/null and b/R/analysis.tar.gz differ
diff --git a/R/analysis/_region_.tex b/R/analysis/_region_.tex
new file mode 100644
index 0000000000000000000000000000000000000000..39ec4ccec2c79d48be50c396316e894facac8278
--- /dev/null
+++ b/R/analysis/_region_.tex
@@ -0,0 +1,5 @@
+\message{ !name(lmer.output.figs.R.tex)}
+\message{ !name(lmer.output.figs.R) !offset(305) }
+ddply(DF,'id',.fun=function(data) data$d[which.max(data$g)])
+
+\message{ !name(lmer.output.figs.R.tex) !offset(-5) }
diff --git a/R/analysis/glmer.run.R b/R/analysis/glmer.run.R
new file mode 100644
index 0000000000000000000000000000000000000000..b9739f15dae810e14565558d56d1250166bcabd4
--- /dev/null
+++ b/R/analysis/glmer.run.R
@@ -0,0 +1,192 @@
+###########################
+###########################
+### FUNCTION TO RUN GLMER ESTIMATION
+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, ...){
+ for(trait in traits)
+ run.multiple.model.for.set.one.trait(model.file,fun.model, trait, set, type.filling=type.filling, ...)
+}
+
+run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
+ for (m in model.files)
+ try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling, ...))
+}
+
+
+run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
+ fun.model <- match.fun(fun.model)
+ for (e in ecoregions)
+ try(fun.model(model.file, trait, set, e, type.filling=type.filling,...))
+}
+
+
+#=====================================================================
+# Run glmer() (in package lme4) for one ecoregion, one trait, one model
+#=====================================================================
+model.files.glmer.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.R.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.ER.R")
+model.files.glmer.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.AD.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.R")
+
+model.files.glmer.Tf.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R.",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
+model.files.glmer.Tf.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
+
+
+fun.test.if.multi.census <- function(data){
+return("tree.id" %in% names(data))
+}
+
+fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
+ Start <- Sys.time()
+ glmer.output <- glmer(formula=formula,data=df.lmer,family=binomial)
+ end <- Sys.time()
+ print(end - Start)
+ print(summary(glmer.output))
+ saveRDS(glmer.output,file=file.path(path.out, "glmer.results.rds"))
+}
+
+run.glmer <- function (model.file, trait, set, ecoregion,
+ min.obs=10, sample.size=NA,
+ type.filling) {
+ require(lme4)
+ source(model.file, local = TRUE)
+ model <- load.model()
+ #= Path for output
+ path.out <- output.dir.glmer(model$name, trait, set,
+ ecoregion,type.filling=type.filling)
+ print(path.out)
+ dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
+ cat("run glmer for model",model.file," for set",
+ set,"ecoregion",ecoregion,"trait",
+ trait,"\n")
+ df.glmer <- load.and.prepare.data.for.glmer(trait, set, ecoregion,
+ min.obs, sample.size,
+ type.filling=type.filling) # return a DF
+ cat("Ok data with Nobs",nrow(df.glmer),"\n")
+ #= Run model
+ fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
+}
+#========================================================================
+
+
+output.dir.glmer <- function (model, trait, set, ecoregion,type.filling) {
+ file.path("output/glmer", set,ecoregion,trait,type.filling,model)
+}
+
+
+#============================================================
+# Function to prepare data for lmer
+#============================================================
+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)
+ fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
+}
+
+fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
+ min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
+## remove tree with NA
+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 no missing traits
+data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
+ !is.na(data.tree[[BATOT]])),]
+# 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
+data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
+return(data.tree)
+}
+
+fun.center.and.standardized.var <- function(x){
+return((x-mean(x))/sd(x))
+}
+
+### get variables for lmer
+fun.get.the.variables.for.glmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
+dead <- data.tree[["dead"]]
+logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
+species.id <- unclass(factor(data.tree[["sp"]]))
+tree.id <- unclass(factor(data.tree[["tree.id"]]))
+plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
+#= multiply CWMs by BATOT
+sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
+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,
+ tree.id=tree.id,
+ plot.species.id=plot.species.id,
+ sumTnTfBn.diff=sumTnTfBn.diff,
+ sumTnTfBn.abs=sumTnTfBn.abs,
+ Tf=data.tree[[tf]],
+ sumTnBn=sumTnBn,
+ sumTfBn=sumTfBn,
+ sumBn=sumBn))
+}
+
+fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
+dead <- data.tree[["dead"]]
+logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
+species.id <- unclass(factor(data.tree[["sp"]]))
+tree.id <- unclass(factor(data.tree[["tree.id"]]))
+plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
+#= multiply CWMs by BATOT
+sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
+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,
+ plot.species.id=plot.species.id,
+ sumTnTfBn.diff=sumTnTfBn.diff,
+ sumTnTfBn.abs=sumTnTfBn.abs,
+ Tf=data.tree[[tf]],
+ sumTnBn=sumTnBn,
+ sumTfBn=sumTfBn,
+ sumBn=sumBn))
+}
+
+#============================================================
+# Function to prepare data for lmer with new CWM data
+# 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 ")
+# get vars names
+CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
+abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
+tf <- paste(trait,"focal",sep=".")
+BATOT <- "BATOT.log"
+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
+glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
+ return(glmer.data)
+}
+
+
+
diff --git a/R/analysis/lmer.output.figs.R b/R/analysis/lmer.output.figs.R
index f044784f531adc0e875ea0ef42884487d67fc701..dc5591adfe623499ae599761cc34a4b962148b48 100644
--- a/R/analysis/lmer.output.figs.R
+++ b/R/analysis/lmer.output.figs.R
@@ -506,7 +506,7 @@ 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,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
- var.y='R2c.simplecomp',ylim=c(-0.015,0.06),
+ 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)
dev.off()
@@ -537,9 +537,15 @@ 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='R2c.simplecomp',ylim=c(-0.015,0.065))
+ var.x='MAP',var.y='R2m.ratio',ylim=c(-0.5,3))
dev.off()
+
+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.06))
+
+
pdf('figs/R2.boxplot.pdf',width=8,height=6)
fun.plot.panel.lmer.res.boxplot(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
@@ -550,7 +556,7 @@ dev.off()
### plot parameters
-models <- c('lmer.LOGLIN.ER','lmer.LOGLIN.E','lmer.LOGLIN.R','lmer.LOGLIN.AD')
+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'),
diff --git a/R/analysis/lmer.run.R b/R/analysis/lmer.run.R
index 216a991fe0b8f194a6e8c5d31d8afbb7b6e8f040..15d94cc425f957d3b763cec3cbc5bbfd330dac1b 100644
--- a/R/analysis/lmer.run.R
+++ b/R/analysis/lmer.run.R
@@ -51,7 +51,7 @@ return("tree.id" %in% names(data))
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
- lmer.output <- lmer(formula=formula,data=df.lmer)
+ lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
end <- Sys.time()
print(end - Start)
print(summary(lmer.output))
@@ -67,6 +67,7 @@ run.lmer <- function (model.file, trait, set, ecoregion,
#= Path for output
path.out <- output.dir.lmer(model$name, trait, set,
ecoregion,type.filling=type.filling)
+ print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run lmer for model",model.file," for set",
set,"ecoregion",ecoregion,"trait",
@@ -86,10 +87,6 @@ run.lmer <- function (model.file, trait, set, ecoregion,
}
#========================================================================
-output.dir <- function (model, trait, set, ecoregion) {
- file.path("output/jags",set,ecoregion,trait, model)
-}
-
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
file.path("output/lmer", set,ecoregion,trait,type.filling,model)
}
@@ -106,8 +103,8 @@ load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
-fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs,
- min.perc.genus=0.90,min.BA.G=-50,max.BA.G=150){
+fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
+ min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
@@ -120,8 +117,8 @@ data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
# 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.genus
-data.tree <- subset(data.tree,subset=data.tree[[perc.genus]] > min.perc.genus & !is.na(data.tree[[perc.genus]]))
+# 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)
}
@@ -186,12 +183,12 @@ return(data.frame(logG=logG,
fun.data.for.lmer <- 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 ")
# get vars names
-CWM.tn <- paste(trait,"CWM",type.filling,"log",sep=".")
-abs.CWM.tntf <- paste(trait,"abs.CWM",type.filling,"log",sep=".")
+CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
+abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
-perc.genus <- paste(trait,"perc.genus",sep=".")
-data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs)
+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
if (length(table(table(data.tree[["tree.id"]])))>1){
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..42cf802629c60a4968d5048c04e9947744411fab
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.AD.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..efa145ff4fd2144d600ea7fc498bf85032ea7054
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.E.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..417358525ceec40d95ee72588943ae050d2eedee
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.ER.AD.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..04663a45c8e08dd74d06faf8d72dd4cb02604532
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.ER.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..f6429cb1d2fe8ecc46c09023a042261abe39da2e
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.R.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..e1a7144731d62211603e8eb450ff10131abe8471
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.simplecomp.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..7a12a8126cd5b3991564e438c01c5951cfc01a98
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R
@@ -0,0 +1,7 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.simplecomp.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn"))
+}
+
+
+
diff --git a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R
new file mode 100644
index 0000000000000000000000000000000000000000..ca20de257b13041d31da665becdeb23b6d40055c
--- /dev/null
+++ b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R
@@ -0,0 +1,5 @@
+load.model <- function () {
+ list(name="lmer.LOGLIN.ER.AD.Tf",
+ lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs"),
+ lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumTfBn+sumBn+sumTnBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
index 29e481a0793f5b33e1fce8ee5c569d23783ea01e..040d4bce6eecff865f26040ceef648fddf764dd0 100644
--- a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
+++ b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
@@ -1,7 +1,7 @@
load.model <- function () {
list(name="lmer.LOGLIN.ER.Tf",
- lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnBn+sumTfBn"),
- lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumBn+sumTnBn+sumTfBn"))
+ lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn"),
+ lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumTfBn+sumBn+sumTnBn"))
}
diff --git a/R/analysis/run.local.R b/R/analysis/run.local.R
index 2d541f55dfb3f0b39349426db35ed521ba4261eb..f17aefb33ff24518727d0d07927e21e86a590574 100644
--- a/R/analysis/run.local.R
+++ b/R/analysis/run.local.R
@@ -1,88 +1,65 @@
-#####
+##### SCRIPT TO TEST BEFORE TO SEDN ON CLUSTER
source("R/analysis/lmer.run.R")
-sets.inv <- c("NVS","US","Sweden","France","NSW","Canada","Spain","Swiss") #
-sets.trop <- c("BCI","Fushan","Paracou","Mbaiki","Luquillo")
-#==========
-# Test lmer
-model.files <- c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.R.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.E.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.HD.R")
-library(multicore)
-options(cores = 5)
-
-run.models.for.set.all.traits('France',model.files,run.lmer,type.filling="species")
-
#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-
-df.lmer <- load.and.prepare.data.for.lmer(trait='SLA',set='France',ecoregion='AB',
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='genus')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.HD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+
+df.lmer <- load.and.prepare.data.for.lmer(trait='SLA',set='France',ecoregion='HI',
min.obs=10, sample.size=NA,
- type.filling='fill') # return a DF
+ type.filling='species') # return a DF
-simple <- readRDS("output/lmer/lmer.LOGLIN.simplecomp/SLA/France/AB/fill/results.rds")
-nocomp <- readRDS("output/lmer/lmer.LOGLIN.nocomp/SLA/France/AB/fill/results.rds")
-ERcomp <- readRDS("output/lmer/lmer.LOGLIN.ER/SLA/France/AB/fill/results.rds")
+## simple <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp/results.rds")
+## nocomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp/results.rds")
+## ERcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER/results.rds")
+## ADcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD/results.rds")
+simple.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp.Tf/results.rds")
+nocomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp.Tf/results.rds")
+ERcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.Tf/results.rds")
+ADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD.Tf/results.rds")
+ERADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.AD.Tf/results.rds")
+
+
+anova(simple,nocomp,ERcomp,ADcomp,simple.Tf,nocomp.Tf,ERcomp.Tf,ADcomp.Tf,ERADcomp.Tf)
res.fix.simple <- fitted(simple)+resid(simple) -simple@pp$X %*% fixef(simple)
res.fix.no <- fitted(nocomp)+resid(nocomp) -nocomp@pp$X %*% fixef(nocomp)
res.fix.ER <- fitted(ERcomp)+resid(ERcomp) -ERcomp@pp$X %*% fixef(ERcomp)
- plot(df.lmer$sumBn,res.fix.no)
+par(mfrow=c(2,2))
+ plot(df.lmer$sumBn,res.fix.no,cex=0.1)
lines(0:15,0:15*fixef(simple)['sumBn'],col='red')
- plot(df.lmer$sumTnBn,res.fix.simple)
+ plot(df.lmer$sumTnBn,res.fix.simple,cex=0.1)
lines((-10):10,(-10):10*fixef(ERcomp)['sumTnBn'],col='red')
- plot(df.lmer$sumTfBn,res.fix.simple)
+ plot(df.lmer$sumTfBn,res.fix.simple,cex=0.1)
lines((-10):10,(-10):10*fixef(ERcomp)['sumTfBn'],col='red')
-
-lapply(c(sets.inv),run.models.for.set.all.traits,model.files,run.lmer,type.filling="species")
+ plot(df.lmer$sumTnTfBn.abs,res.fix.simple,cex=0.1)
+lines((-10):10,(-10):10*fixef(ADcomp)['sumTnTfBn.abs'],col='red')
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="species")
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="genus")
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="fill")
+ plot(df.lmer$Tf,res.fix.ER,cex=0.1)
+lines((-10):10,(-10):10*fixef(ERcomp.Tf)['Tf'],col='red')
-
-## # load results
-## mod.AD <- readRDS("output/lmer/lmer.LOGLIN.AD/SLA/France/AB/results.rds")
-## mod.HD <- readRDS("output/lmer/lmer.LOGLIN.HD/SLA/France/AB/results.rds")
-## mod.nocomp <- readRDS("output/lmer/lmer.LOGLIN.nocomp/SLA/France/AB/results.rds")
-## mod.simplecomp <- readRDS("output/lmer/lmer.LOGLIN.simplecomp/SLA/France/AB/results.rds")
-## # summary
-## summary(mod.AD)
-## summary(mod.HD)
-## summary(mod.nocomp)
-## summary(mod.simplecomp)
-## # predictions (back-transforming the log)
-## pred.AD <- exp(fitted(mod.AD)+var(residuals(mod.AD))/2)
-## pred.HD <- exp(fitted(mod.HD)+var(residuals(mod.HD))/2)
-## pred.nocomp <- exp(fitted(mod.nocomp)+var(residuals(mod.nocomp))/2)
-## pred.simplecomp <- exp(fitted(mod.simplecomp)+var(residuals(mod.simplecomp))/2)
+### GLMR
-## # observations
-## str(mod.AD)
-## obs <- exp(mod.AD@frame$logG)
+source("R/analysis/glmer.run.R")
-## #= R2 function
-## r2.calc <- function (obs,pred) {
-## SCR <- sum((obs-pred)^2,na.rm=TRUE)
-## SCT <- sum((obs-mean(obs))^2,na.rm=TRUE)
-## R.square <- 1-SCR/SCT
-## return(R.square)
-## }
-
-## # r2
-## r2.nocomp <- r2.calc(obs,pred.nocomp) # -0.48, Very bad model due to log and absence of competition
-## r2.AD <- r2.calc(obs,pred.AD) # 0.60 !
-## r2.HD <- r2.calc(obs,pred.HD) # -0.52
-## r2.simplecomp <- r2.calc(obs,pred.simplecomp) # -0.52
+#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
diff --git a/R/process.data/test.tree.CWM.R b/R/process.data/test.tree.CWM.R
index e6b350980e88fe98788570489d63b2f5567a30de..650ce4ddb333068201437cb1ad89ccbdaa4670d7 100644
--- a/R/process.data/test.tree.CWM.R
+++ b/R/process.data/test.tree.CWM.R
@@ -280,8 +280,10 @@ var.name.fill <- paste(traits.name,"abs.CWM.fill",sep=".")
var.name.perc.genus <- paste(traits.name,"perc.genus",sep=".")
var.name.perc.species <- paste(traits.name,"perc.species",sep=".")
perc.non.missing.sp <- lapply(1:length(var.name.sp),fun.perc.non.missing,var.name.sp,var.name.perc.species,data)
-num.non.missing.sp <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,var.name.perc.species,data)
-num.non.missing.genus <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,var.name.perc.genus,data)
+num.non.missing.sp <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,
+ var.name.perc.species,data)
+num.non.missing.genus <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,
+ var.name.perc.genus,data)
perc.non.missing.genus <- lapply(1:length(var.name.sp),fun.perc.non.missing,var.name.genus,var.name.perc.genus,data)
perc.genus <- lapply(var.name.perc.genus,function(i,data) mean(data[[i]],na.rm=TRUE),data)
perc.species <- lapply(var.name.perc.species,function(i,data) mean(data[[i]],na.rm=TRUE),data)
@@ -366,15 +368,21 @@ mat.perc <- data.frame(lapply(mat.perc, function(x) (unlist(x))))
write.csv(mat.perc,file=file.path(filedir, "all.sites.perc.traits.csv"), row.names=FALSE)
mat.perc <- read.csv(file=file.path(filedir, "all.sites.perc.traits.csv"))
-mat.num <- mat.perc[,c(1:3,9:13)]
-mat.num[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
+mat.num.g <- mat.perc[,c(1:3,9:13)]
+mat.num.g[,4:8] <- mat.perc[,9:13]/mat.perc[,3]
-names(mat.num) <- c('set','ecoregion','P obs total','P Leaf N','P Seed mass','P SLA','P Wood density','P Max height')
+names(mat.num.g) <- c('set','ecoregion','P obs total','P Leaf N','P Seed mass','P SLA','P Wood density','P Max height')
library('pander')
-pandoc.table(mat.num,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+pandoc.table(mat.num.g,caption="Number of tree radial growth observation per data sets and ecoregion.",split.tables='Inf')
mat.perc <- read.csv(file=file.path(filedir, "all.sites.perc.traits.csv"))
-mat.num <- mat.perc[,c(1:8)]
-mat.num[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
-pandoc.table(mat.num,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+mat.num.sp <- mat.perc[,c(1:8)]
+mat.num.sp[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
+pandoc.table(mat.num.sp,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+
+par(mfrow=c(2,2))
+plot(mat.num.sp[,4],mat.num.g[,4])
+plot(mat.num.sp[,5],mat.num.g[,5])
+plot(mat.num.sp[,6],mat.num.g[,6])
+plot(mat.num.sp[,7],mat.num.g[,7])
diff --git a/launch_all_glmer.bash b/launch_all_glmer.bash
new file mode 100644
index 0000000000000000000000000000000000000000..10ffefc12f76fec6abb24d4f8b92478659218c25
--- /dev/null
+++ b/launch_all_glmer.bash
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+export LD_LIBRARY_PATH=/usr/lib64/R/library
+
+mkdir -p trait.workshop
+
+for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies1$site.sh
+ qsub trait.workshop/glmerspecies1$site.sh -l nodes=1:ppn=1 -N "glmerspecies1$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies2$site.sh
+ qsub trait.workshop/glmerspecies2$site.sh -l nodes=1:ppn=1 -N "glmerspecies2$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus1$site.sh
+ qsub trait.workshop/glmergenus1$site.sh -l nodes=1:ppn=1 -N "glmergenus1$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus2$site.sh
+ qsub trait.workshop/glmergenus2$site.sh -l nodes=1:ppn=1 -N "glmergenus2$site" -q opt32G -j oe
+
+done
+
diff --git a/launch_all_lmer.bash b/launch_all_lmer.bash
index c809594d5eb00ad32b2a4b0a73be570a96a8a8e8..adbe3879809b83a3305a0ebd0a3269741db08949 100644
--- a/launch_all_lmer.bash
+++ b/launch_all_lmer.bash
@@ -6,17 +6,17 @@ mkdir -p trait.workshop
for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');sessionInfo();run.models.for.set.all.traits($site,model.files.lmer.1, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill1$site.sh
- qsub trait.workshop/fill1$site.sh -l nodes=1:ppn=1 -N "$site-1" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species');\n\"" > trait.workshop/species1$site.sh
+ qsub trait.workshop/species1$site.sh -l nodes=1:ppn=1 -N "lmerspecies1$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.2, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill2$site.sh
- qsub trait.workshop/fill2$site.sh -l nodes=1:ppn=1 -N "$site-2" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species');\n\"" > trait.workshop/species2$site.sh
+ qsub trait.workshop/species2$site.sh -l nodes=1:ppn=1 -N "lmerspecies2$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill3$site.sh
- qsub trait.workshop/fill3$site.sh -l nodes=1:ppn=1 -N "$site-3" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus1$site.sh
+ qsub trait.workshop/genus1$site.sh -l nodes=1:ppn=1 -N "lmergenus1$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill4$site.sh
- qsub trait.workshop/fill4$site.sh -l nodes=1:ppn=1 -N "$site-4" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus2$site.sh
+ qsub trait.workshop/genus2$site.sh -l nodes=1:ppn=1 -N "lmergenus2$site" -q opt32G -j oe
done