Skip to content

GitLab

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

progress on stan with random species f n

parent 6c1cc9a2
Hide whitespace changes
Inline Side-by-side
Makefile
View file @ edf742ae
......@@ -5,7 +5,7 @@ D4 := figs/test.format.tree
D5 := figs/test.traits
D6 := figs/test.CWM
sites:= US Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS
D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) ))
D3Done := $(addsuffix /Done.no.txt,$(addprefix $(D3)/, $(sites) ))
D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) ))
D2tree := $(addsuffix /tree.csv,$(addprefix $(D2)/, $(sites) ))
......@@ -42,7 +42,11 @@ $(D2)/TRY/data.TRY.std.rds:
#-------------------------------------------------------
GLOBAL: $(D3)/Done.txt
$(D3)/Done.txt: R/process.data/merge.all.processed.data.R #sites R/process.data/process-fun.R # $(D5)/Done.txt
$(D3)/Done.txt: R/process.data/remove.out.R $(D3)/Done.t.txt
Rscript $<
$(D3)/Done.t.txt: R/process.data/merge.all.processed.data.R
Rscript $<
#-------------------------------------------------------
......@@ -65,11 +69,11 @@ $(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D3Done)
#-------------------------------------------------------
#-------------------------------------------------------
MERGE.ALL: output/processed/data.all.csv
# #-------------------------------------------------------
# MERGE.ALL: output/processed/data.all.csv
output/processed/data.all.csv: R/process.data/merge.all.processed.data.R
Rscript $<
# output/processed/data.all.csv: R/process.data/merge.all.processed.data.R
# Rscript $<
#-------------------------------------------------------
......
R/analysis/glmer.run.R
View file @ edf742ae
......@@ -6,19 +6,22 @@ library(lme4)
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait,
type.filling, fname = 'data.all.no.log.csv',
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
type.filling,
fname = 'data.all.no.log.csv',
cat.TF = FALSE, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait,
type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
try(run.model.for.set.one.trait (m, fun.model, trait,
type.filling = type.filling,
fname = fname, cat.TF = cat.TF, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait,
run.model.for.set.one.trait <- function(model.file, fun.model, trait,
type.filling, fname, cat.TF, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
try(fun.model(model.file, trait, type.filling = type.filling,
fname = fname, cat.TF = cat.TF, ...))
}
......@@ -27,42 +30,65 @@ run.model.for.set.one.trait <- function(model.file,fun.model, trait,
#=====================================================================
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.HD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
c("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.AD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
model.files.glmer.Tf.3 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.MAT.MAP.R")
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.MAT.MAP.R")
model.files.glmer.Tf.4 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.MAT.MAP.R",
## "R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.MAT.MAP.R")
fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
### TODO ER QND QD TO CHQNGE NEED TO ADD CAT
fun.call.glmer.and.save <- function(formula, df.lmer, path.out){
library(MASS)
logexp <- function(exposure = 1)
{
linkfun <- function(mu) qlogis(mu^(1/exposure))
## FIXME: is there some trick we can play here to allow
## evaluation in the context of the 'data' argument?
linkinv <- function(eta) plogis(eta)^exposure
logit_mu_eta <- function(eta) {
ifelse(abs(eta)>30, .Machine$double.eps,
exp(eta)/(1+exp(eta))^2)
## OR .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
}
mu.eta <- function(eta) {
exposure * plogis(eta)^(exposure-1) *
logit_mu_eta(eta)
}
valideta <- function(eta) TRUE
link <- paste("logexp(", deparse(substitute(exposure)), ")",
sep = "")
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta,
name = link),
class = "link-glm")
}
Start <- Sys.time()
glmer.output <- glmer(formula=formula,data=df.lmer,
family=binomial(link = 'cloglog'))
## glmer.output <- glmer(formula = formula, data = df.lmer,
## family = binomial(link = 'cloglog'))
glmer.output <- glmer(formula = formula, data = df.lmer,
family = binomial(link = logexp(exp(df.lmer$logyear))))
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.global.rds"))
# add starting value ? start = list(theta = 1, fixef = c(1, 0))
summary(mfit_glmer)
saveRDS(glmer.output, file = file.path(path.out, "glmer.results.global.rds"))
}
run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA,
min.obs = 5, sample.size = NA,
type.filling, fname , cat.TF) {
require(lme4)
source(model.file, local = TRUE)
......@@ -71,36 +97,34 @@ run.glmer <- function (model.file, trait,
if(fname == 'data.all.no.log.csv') dir <- 'all.no.log'
if(fname == 'data.all.csv') dir <- 'all.log'
path.out <- output.dir.glmer(dir, model$name, trait,
type.filling=type.filling)
type.filling = type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run glmer for model",model.file," for set",
'all',"trait",
trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait,
min.obs, sample.size,
type.filling=type.filling,
fname = fname,
cat.TF = cat.TF) # return a DF
cat("Ok data with Nobs",nrow(df.glmer),"\n")
dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
cat("run glmer for model", model.file, " for set",
'all', "trait",
trait, "\n")
df.glmer <- load.data.for.glmer(trait,
fname = fname,
cat.TF = cat.TF) # 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)
fun.call.glmer.and.save(formula = model$glmer.formula, df.glmer, path.out)
}
#========================================================================
output.dir.glmer <- function (dir, model, trait, type.filling) {
file.path("output/glmer", dir,trait,type.filling,model)
file.path("output/glmer", dir, trait, type.filling, model)
}
fun.load.data.all <- function(base.dir,fname){
fun.load.data.all <- function(base.dir, fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
nrows <- as.numeric(system(paste('wc -l < ', file.path(base.dir, fname)),
intern = TRUE))
data.tree.tot <- read.csv(file = file.path(base.dir, fname),
......@@ -114,62 +138,67 @@ return(data.tree.tot)
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait,
min.obs, sample.size, type.filling,
load.data.for.glmer <- function(trait, cat.TF = FALSE,
fname = 'data.all.no.log.csv',
base.dir = "output/processed/",
cat.TF = FALSE){
base.dir = "output/processed/"){
### load data
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- readRDS(file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
df <- readRDS(file = file.path(base.dir,
paste('data.glmer', trait, type,
cat,'rds', sep = '.')))
return(df)
}
load.and.save.data.for.glmer <- function(trait,
min.obs = 10, sample.size = NA, type.filling = 'species',
fname = 'data.all.no.log.csv',
base.dir = "output/processed/", cat.TF = FALSE){
min.obs = 10, sample.size = NA,
type.filling = 'species',
fname = 'data.all.no.log.csv',
base.dir = "output/processed/",
cat.TF = FALSE){
### load data
data.tree.tot <- fun.load.data.all(base.dir,fname)
data.tree.tot <- fun.load.data.all(base.dir, fname)
if (fname == 'data.all.no.log.csv') type <- 'no.log'
if (fname == 'data.all.csv') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling, cat.TF = cat.TF)
saveRDS(df,file = file.path(base.dir,paste('data.glmer', trait, type, cat,'rds', sep = '.')))
df <- fun.data.for.glmer(data.tree.tot,trait, type.filling = type.filling,
cat.TF = cat.TF)
saveRDS(df, file = file.path(base.dir, paste('data.glmer', trait,
type, cat,'rds',
sep = '.')))
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,
fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs,
BATOT, min.obs,
min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) &
(!is.na(data.tree[['year']])))
data.tree <- subset(data.tree, subset = (!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) &
(!is.na(data.tree[['year']])))
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
data.tree <- subset(data.tree, subset = data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
!is.na(data.tree[[BATOT]])), ]
# select obs abov min perc.neigh species
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
min.perc.neigh.sp &
!is.na(data.tree[[perc.neigh.sp]]))
# select obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs & !is.na(data.tree[[perc.neigh.gs]]))
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
data.tree <- subset(data.tree, subset = data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(
data.tree[["sp"]]))>min.obs])
return(data.tree)
......@@ -184,8 +213,9 @@ return(x/sd(x))
}
### get variables for lmer
fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,
abs.CWM.tntf,tf,min.BA.G=50){
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"]]))
logyear <- log(data.tree[["year"]])
......@@ -202,28 +232,31 @@ 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,
logyear=logyear,
MAT = MAT,
MAP = MAP,
species.id=species.id,
tree.id=tree.id,
plot.id=plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.diff= fun.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs= fun.standardized.var(sumTnTfBn.abs),
Tf= fun.standardized.var(data.tree[[tf]]),
sumTnBn=fun.standardized.var(sumTnBn),
sumTfBn=fun.standardized.var(sumTfBn),
sumBn=fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
return(
data.frame(dead = dead,
logD = logD,
logyear = logyear,
MAT = MAT,
MAP = MAP,
species.id = species.id,
tree.id = tree.id,
plot.id = plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.diff = fun.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs = fun.standardized.var(sumTnTfBn.abs),
Tf = fun.standardized.var(data.tree[[tf]]),
sumTnBn = fun.standardized.var(sumTnBn),
sumTfBn = fun.standardized.var(sumTfBn),
sumBn = fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
}
### get variables for lmer
fun.get.the.variables.for.glmer.no.tree.id.cat <- function(data.tree,BATOT,CWM.tn,
abs.CWM.tntf,tf,min.BA.G=50){
fun.get.the.variables.for.glmer.no.tree.id.cat <- 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"]]))
logyear <- log(data.tree[["year"]])
......@@ -252,34 +285,34 @@ sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *sumTfBn
sumBn <- data.tree[[BATOT]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
return(data.frame(dead=dead,
logD=logD,
logyear=logyear,
return(data.frame(dead = dead,
logD = logD,
logyear = logyear,
MAT = MAT,
MAP = MAP,
species.id=species.id,
tree.id=tree.id,
plot.id=plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.abs.A_EV= fun.standardized.var(sumTnTfBn.abs.A_EV),
sumTnTfBn.abs.A_D= fun.standardized.var(sumTnTfBn.abs.A_D),
sumTnTfBn.abs.C= fun.standardized.var(sumTnTfBn.abs.C),
Tf.A_EV= fun.standardized.var(Tf.A_EV),
Tf.A_D= fun.standardized.var(Tf.A_D),
Tf.C= fun.standardized.var(Tf.C),
sumTnBn.A_EV=fun.standardized.var(sumTnBn.A_EV),
sumTnBn.A_D=fun.standardized.var(sumTnBn.A_D),
sumTnBn.C=fun.standardized.var(sumTnBn.C),
sumTfBn.A_EV=fun.standardized.var(sumTfBn.A_EV),
sumTfBn.A_D=fun.standardized.var(sumTfBn.A_D),
sumTfBn.C=fun.standardized.var(sumTfBn.C),
sumBn=fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
species.id = species.id,
tree.id = tree.id,
plot.id = plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.abs.A_EV = fun.standardized.var(sumTnTfBn.abs.A_EV),
sumTnTfBn.abs.A_D = fun.standardized.var(sumTnTfBn.abs.A_D),
sumTnTfBn.abs.C = fun.standardized.var(sumTnTfBn.abs.C),
Tf.A_EV = fun.standardized.var(Tf.A_EV),
Tf.A_D = fun.standardized.var(Tf.A_D),
Tf.C = fun.standardized.var(Tf.C),
sumTnBn.A_EV = fun.standardized.var(sumTnBn.A_EV),
sumTnBn.A_D = fun.standardized.var(sumTnBn.A_D),
sumTnBn.C = fun.standardized.var(sumTnBn.C),
sumTfBn.A_EV = fun.standardized.var(sumTfBn.A_EV),
sumTfBn.A_D = fun.standardized.var(sumTfBn.A_D),
sumTfBn.C = fun.standardized.var(sumTfBn.C),
sumBn = fun.standardized.var(sumBn),
stringsAsFactors = FALSE))
}
......@@ -287,16 +320,16 @@ return(data.frame(dead=dead,
# 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', cat.TF = FALSE) {
fun.data.for.glmer <- function(data.tree, trait, min.obs = 10,
type.filling = 'species', cat.TF = FALSE) {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA",
"Wood.density","Max.height"))
"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',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".")
CWM.tn <- paste(trait, "CWM", 'fill', sep = ".")
abs.CWM.tntf <- paste(trait, "abs.CWM", 'fill', sep = ".")
tf <- paste(trait, "focal", sep = ".")
BATOT <- "BATOT"
perc.neigh.sp <- paste(trait, "perc", 'species', sep = ".")
perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
......@@ -304,7 +337,7 @@ data.tree <- fun.select.data.for.analysis(data.tree,
abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs)
BATOT, min.obs)
#= DATA LIST FOR JAGS
if(!cat.TF) {
print('no cat')
......
R/analysis/lmer.all.log.output.R 0 → 100644
View file @ edf742ae
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.log.rds')
R/analysis/lmer.all.no.log.output.R
View file @ edf742ae
......@@ -8,8 +8,8 @@ traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
for(model in c(model.files.lmer.Tf.1, model.files.lmer.Tf.2,
model.files.lmer.Tf.3, model.files.lmer.Tf.4,
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
source(model, local = TRUE)
model.obj <- load.model()
......@@ -17,10 +17,10 @@ for (trait in traits){
type.filling=type.filling)
files <- c(files,file.path(pathout,"results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list)
......
R/analysis/lmer.all.output.figs.R 0 → 100644
View file @ edf742ae
#!/usr/bin/env Rscript
rm(list = ls())
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.R")
source("R/utils/plot.R")
source('R/utils/mem.R')
## load results all
list.all.results <- readRDS('output/list.lmer.out.all.no.log.rds')
print(names(list.all.results))
## list.all.results <- readRDS('output/list.lmer.out.all.rds')
## names.param <- unique(unlist(lapply(list.all.results,
## function(x) names(x$lmer.summary$fixed.coeff.E))))
## names.param2 <- names.param
## names.param2[names.param2 == "sumBn:MAT"] <- "MAT:sumBn"
## names.param2[names.param2 == "sumBn:MAP"] <- "MAP:sumBn"
## names(names.param2) <-names.param
## DF.global <- do.call( 'rbind', lapply(list.all.results,
## fun.merge.results.global,
## names.param = names.param2))
## if(sum(DF.global$conv != 0) >0) stop( 'bad convergence')
## ## AIC table
## for (i in c('Wood.density', 'SLA', 'Leaf.N', 'Max.height')){
## DF.t <- DF.global[DF.global$trait == i, ]
## print(i)
## print(DF.t$model[DF.t$AIC == min(DF.t$AIC)])
## res <- data.frame(model = DF.t$model,delta.aic = DF.t$AIC - min(DF.t$AIC))
## print(pander(res))
## }
## plot the parameters
## add set random effect estim
pdf('figs/param.global.model.ER.AD.3typesb.pdf')
#x11()
opar <- par()
m <- matrix(c(1:4,5,5), 2, 3)
layout(m, widths=c(4, 4, 1),
heights=c(4, 4))
for (i in c('Wood.density', 'SLA', 'Leaf.N', 'Max.height')){
param.vec <- c("logD", "Tf","sumBn", "sumTnBn","sumTfBn", "sumTnTfBn.abs")
list.temp <- list.all.results[[paste("all.no.log_", i ,
"_species_lmer.LOGLIN.ER.AD.Tf",
sep = '')]]$lmer.summary
param.mean <- list.temp$fixed.coeff.E[param.vec]
param.std <- list.temp$fixed.coeff.Std.Error
names(param.std) <- names(list.temp$fixed.coeff.E)
param.std <- param.std[param.vec]
param.BLUP <- list.temp$set.BLUP
print(i)
print(param.mean)
print(param.std)
print(list.temp$fixed.coeff.Std.Error)
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-0.6, 0.6))
abline(h = 0)
axis(1, 1:length(param.vec), labels = param.vec, las = 3)
fun.plot.error.bar(1:length(param.vec),
param.mean,
param.std)
for (j in 1:length(param.vec)){
points(rep(j,nrow(param.BLUP)), param.BLUP[,param.vec[j]],
pch = pch.vec[rownames(param.BLUP)],
col = col.vec[rownames(param.BLUP)])
}
}
par(mar=c(0, 0, 0, 0))
plot(1,1,pch = NA, xaxt = 'n', yaxt = 'n')
legend("center",legend = names(col.vec), col = col.vec, pch = pch.vec,
bty = 'n')
par(opar)