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

survival and growth model with lme4

parent cdf19fdf
\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) }
###########################
###########################
### 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)
}
......@@ -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'),
......
......@@ -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)
......
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"))
}
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"))
}
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"))
}
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"))
}
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"))
}
load.model <- function () {
list(name="glmer.LOGLIN.simplecomp.Tf",
glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)"))
}
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"))
}
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"))
}
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"))
}
......
#####
##### 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')
......@@ -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])
#!/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
......@@ -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