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

progress on slide

parent 4516d3ce
......@@ -37,45 +37,38 @@ include mk/cible.tree.traits.mk
TRY: $(D2)/TRY/data.TRY.std.rds
$(D2)/TRY/data.TRY.std.rds:
Rscript R/format.data/TRY.R
Rscript scripts/format.data/TRY.R
#-------------------------------------------------------
GLOBAL: $(D3)/Done.txt
$(D3)/Done.txt: R/process.data/remove.out.R $(D3)/Done.t.txt
$(D3)/Done.txt: scripts/process.data/remove.out.R $(D3)/Done.t.txt
Rscript $<
$(D3)/Done.t.txt: R/process.data/merge.all.processed.data.R
$(D3)/Done.t.txt: scripts/process.data/merge.all.processed.data.R
Rscript $<
#-------------------------------------------------------
TEST.TREE: $(D4)/Done.txt tree.all.sites
$(D4)/Done.txt: R/format.data/test.tree.R $(D3tree)
$(D4)/Done.txt: scripts/format.data/test.tree.R $(D3tree)
Rscript $<
#-------------------------------------------------------
TEST.TRAITS: $(D5)/Done.txt #traits.all.sites
$(D5)/Done.txt: R/find.trait/test.traits.R $(D3traits)
$(D5)/Done.txt: scripts/find.trait/test.traits.R $(D3traits)
Rscript $<
#-------------------------------------------------------
TEST.CWM: $(D6)/Done.txt
$(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D3Done)
$(D6)/Done.txt: scripts/process.data/test.tree.CWM.R $(D3Done)
Rscript $<
#-------------------------------------------------------
# #-------------------------------------------------------
# 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?
......
......@@ -6,20 +6,20 @@ source("R/analysis/glmer.run.all.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
var.sample <- 'ecocode.id'
files <- c()
for (trait in traits){
for(model in c(model.files.glmer.Tf.1, model.files.glmer.Tf.2,
model.files.glmer.Tf.3[-3], model.files.glmer.Tf.4)){
for(model in c(model.files.glmer.Tf.1)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.glmer(model.obj$name, trait,
pathout <- output.dir.glmer('all.no.log', model.obj$name, trait,
type.filling=type.filling)
files <- c(files,file.path(pathout,"glmer.results.global.rds"))
files <- c(files,file.path(pathout,paste(var.sample, "glmer.results.global.rds")))
}
}
out <- lapply(files, summarise.glmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/glmer.global.list.out.rds')
out <- lapply(files, summarise.glmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),function(x) paste(as.vector(x[-5]),collapse="_"))
saveRDS(out,file='output/glmer.global.list.out.ecocode.id.rds')
......@@ -22,7 +22,7 @@ summarise.glmer.output <- function(x){
fixed.var=variance.fixed.glmm.lmer(x))
}
summarise.glmer.all.output <- function(x){
summarise.glmer.all.output <- function(x, random.name = 'biomes.id'){
list( nobs = nobs(x),
R2m =Rsquared.glmm.lmer(x)$Marginal,
R2c =Rsquared.glmm.lmer(x)$Conditional,
......@@ -33,7 +33,7 @@ summarise.glmer.all.output <- function(x){
fixed.coeff.E=fixef(x),
fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
fixed.var=variance.fixed.glmm.lmer(x),
set.BLUP=coef(x)$set.id)
set.BLUP=coef(x)[[random.name]])
}
......@@ -48,11 +48,11 @@ summarise.glmer.output.list <- function(f ){
return(res)
}
summarise.glmer.output.list <- function(f ){
summarise.glmer.output.all.list <- function(f, random.name = 'biomes.id' ){
output.glmer <- read.glmer.output(f)
if(!is.null(output.glmer)){
res <- list(files.details=files.details.all(f),
glmer.summary=summarise.glmer.all.output( output.glmer))
glmer.summary=summarise.glmer.all.output( output.glmer, random.name = random.name))
}else{
res <- NULL
}
......@@ -69,7 +69,7 @@ files.details <- function(x){
files.details.all <- function(x){
s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
names(s) <- c("d1", "d2", "set", "trait", "filling", "model", "file" )
names(s) <- c("d1", "d2", "trait", "filling", "model", "file" )
s[-(1:2)]
}
......
......@@ -8,7 +8,7 @@ library(lme4)
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
type.filling,
fname = 'data.all.no.log.csv',
fname = 'data.all.no.log.rds',
cat.TF = FALSE, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model, trait,
......@@ -32,6 +32,18 @@ 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.ER.AD.Tf.r.species.R")
model.files.glmer.Tf.2 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.r.biomes.species.R")
model.files.glmer.Tf.3 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R")
model.files.glmer.Tf.4 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.biomes.species.R")
model.files.glmer.Tf.5 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species.R")
### TODO ER QND QD TO CHQNGE NEED TO ADD CAT
fun.call.glmer.and.save <- function(formula, df.glmer, path.out, var.sample){
......@@ -67,8 +79,6 @@ logexp <- function(exposure = 1)
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
# add starting value ? start = list(theta = 1, fixef = c(1, 0))
summary(mfit_glmer)
saveRDS(glmer.output,
file = file.path(path.out,
......@@ -84,9 +94,10 @@ run.glmer <- function (model.file, trait,
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
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,
if(fname == 'data.all.no.log.rds') dir <- 'all.no.log'
if(fname == 'data.all.rds') dir <- 'all.log'
path.out <- output.dir.glmer(dir, model = model$name,
trait = trait,
type.filling = type.filling)
print(path.out)
dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
......@@ -130,8 +141,13 @@ if (fname == 'data.all.no.log.rds') type <- 'no.log'
if (fname == 'data.all.log.rds') 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 = '.')))
if(!is.na(sample.size)){
df$species.id[is.na(df$species.id)] <- 'missing.sp'
......
......@@ -8,8 +8,11 @@ traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.4b,
model.files.lmer.Tf.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)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
......@@ -30,53 +33,53 @@ names(out) <- lapply(lapply(files,files.details.all),
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.ecocode.id.no.log.rds')
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.4b,
model.files.lmer.Tf.3)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
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.species.id.no.log.rds')
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.4b,
model.files.lmer.Tf.3)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"obs.id.results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
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.obs.id.no.log.rds')
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
## traits <- c("SLA", "Wood.density","Max.height")
## type.filling <- 'species'
## files <- c()
## for (trait in traits){
## for(model in c(model.files.lmer.Tf.4b,
## model.files.lmer.Tf.3)){
## source(model, local = TRUE)
## model.obj <- load.model()
## pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
## type.filling=type.filling)
## files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
## }
## }
## out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
## 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.species.id.no.log.rds')
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
## traits <- c("SLA", "Wood.density","Max.height")
## type.filling <- 'species'
## files <- c()
## for (trait in traits){
## for(model in c(model.files.lmer.Tf.4b,
## model.files.lmer.Tf.3)){
## source(model, local = TRUE)
## model.obj <- load.model()
## pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
## type.filling=type.filling)
## files <- c(files,file.path(pathout,"obs.id.results.nolog.all.rds"))
## }
## }
## out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
## 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.obs.id.no.log.rds')
......@@ -716,6 +716,7 @@ list.new.data <- list(MAT = new.data.MAT, MAP = new.data.MAP)
return(list.new.data)
}
plot.fun.ci.param.var.climate <- function(list.output, trait, param.vec, var.vec,
n.length = 100, add.mean.effect = FALSE, ...){
list.new.data <- fun.generate.pred.dat(n.length = 100)
......@@ -765,10 +766,6 @@ plot.param.BLUP <- function(list.res,
"sumTfBn", "sumTnTfBn.abs"),
col.vec,
pch.vec){
## m <- rbind(c(0.1, 0.5, 0.55, 1), c(0.5, 0.9, 0.55, 1),
## c(0.1, 0.5, 0.1, 0.45), c(0.5, 0.9, 0.1, 0.45),
## c(0.90, 1, 0.1, 1))
## split.screen(m)
m <- matrix(c(1:4,5,5), 2, 3)
layout(m, widths=c(4, 4, 1),
heights=c(4, 4))
......@@ -784,7 +781,7 @@ for (i in traits){
param.BLUP <- list.temp$set.BLUP
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-0.6, 0.6))
main = i, ylim = c(-0.45, 0.25))
abline(h = 0)
axis(1, 1:length(param.vec), labels = param.vec, las = 3)
fun.plot.error.bar(1:length(param.vec),
......@@ -793,16 +790,54 @@ for (i in traits){
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)])
col = col.vec[rownames(param.BLUP)],
cex= 2)
}
}
## par(mar=c(0, 0, 0, 0))
plot(1,1,pch = NA, xaxt = 'n', yaxt = 'n')
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names(col.vec), col = col.vec, pch = pch.vec,
bty = 'n')
}
##
plot.param <- function(list.res,
model = 'lmer.LOGLIN.ER.AD.Tf',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("logD", "Tf","sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs"),
col.vec,
pch.vec){
m <- matrix(c(1:4), 2, 2)
layout(m, widths=c(4, 4),
heights=c(4, 4))
for (i in traits){
## screen((1:length(traits))[traits == i])
list.temp <- list.res[[paste("all.no.log_", i ,
"_species_", model,
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]
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-0.45, 0.25))
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)
}
}
##
fun.get.conv <- function(list.res.lme4){
......
load.model <- function () {
list(name="glmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.biomes.species",
glmer.formula=formula("dead~1+Tf+Tf:MAT+Tf:MAP+logD+MAT+MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"),
glmer.formula.offset=formula("dead~1+offset(logyear)+Tf+Tf:MAT+Tf:MAP+logD+MAT+MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)")
)
}
load.model <- function () {
list(name="glmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species",
glmer.formula=formula("dead~1+Tf+Tf:MAT+Tf:MAP+logD+MAT+MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"),
glmer.formula.offset=formula("dead~1+offset(logyear)+Tf+Tf:MAT+Tf:MAP+logD+MAT+MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)")
)
}
load.model <- function () {
list(name="glmer.LOGLIN.ER.AD.Tf.r.biomes.species",
glmer.formula=formula("dead~1+Tf+logD+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"),
glmer.formula.offset=formula("dead~1+offset(logyear)+Tf+logD+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)")
)
}
load.model <- function () {
list(name="glmer.LOGLIN.ER.AD.Tf.r.ecocode.species",
glmer.formula=formula("dead~1+Tf+logD+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"),
glmer.formula.offset=formula("dead~1+offset(logyear)+Tf+logD+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)")
)
}
library(rstan) #http://www.mc-stan.org/
source("R/analysis/lmer.run.R")
stan.model <- "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_ecocode;
int<lower=0> N_set;
int<lower=0> N_tree;
int species_id[N_indiv];
int plot_id[N_indiv];
int ecocode_id[N_indiv];
int set_id[N_indiv];
int tree_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
real Tf[N_indiv];
real sumTfBn[N_indiv];
real sumTnBn[N_indiv];
real sumTnTfBn_abs[N_indiv];
real tree_01[N_indiv];
}
parameters{
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean_logD;
real<lower=-10,upper=10> mean_Tf;
real<lower=-10,upper=10> mean_sumBn;
real<lower=-10,upper=10> mean_sumTfBn;
real<lower=-10,upper=10> mean_sumTnBn;
real<lower=-10,upper=10> mean_sumTnTfBn_abs;
real<lower=0.0001,upper=5> sigma_inter_species;
real<lower=0.0001,upper=5> sigma_inter_set;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_inter_tree;
real<lower=0.0001,upper=5> sigma_logD_species;
real<lower=0.0001,upper=5> sigma_Tf_ecocode;
real<lower=0.0001,upper=5> sigma_sumBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumBn_species;
real<lower=0.0001,upper=5> sigma_sumTfBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTnBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_abs_ecocode;
real<lower=0.0001,upper=5> sigma;
vector[N_ecocode] param_Tf_ecocode;
vector[N_ecocode] param_sumBn_ecocode;
vector[N_species] param_sumBn_species;
vector[N_ecocode] param_sumTfBn_ecocode;
vector[N_ecocode] param_sumTnBn_ecocode;
vector[N_ecocode] param_sumTnTfBn_abs_ecocode;
}
transformed parameters {
# vector for prediction
vector[N_species] param_logD_species;
vector[N_species] intercept_species;
vector[N_plot] intercept_plot;
vector[N_set] intercept_set;
vector[N_indiv] theo_g;
vector[N_tree] intercept_tree;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_species[species_id[i]]
+ intercept_tree[tree_id[i]]*tree_01[i]
+ intercept_plot[plot_id[i]]
+ intercept_set[set_id[i]]
+ (mean_logD + param_logD_sp[species_id[i]])*logD[i]
+ (mean_Tf + param_Tf_ecocode[ecocode_id[i]])*Tf[i]
+ (mean_sumBn + param_sumBn_ecocode[ecocode_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (mean_sumTfBn + param_sumTfBn_ecocode[ecocode_id[i]])*sumTfBn[i]
+ (mean_sumTnBn + param_sumTnBn_ecocode[ecocode_id[i]])*sumTnBn[i]
+ (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_ecocode[ecocode_id[i]])*sumTnTfBn_abs[i]
;
}
}
model {
# constants for prior
real sigma0;
################################################################################
######################## growth model with STAN #############################
###############################################
########### Hierarchical parameters ########
### species random param
param_logD_species ~ normal(0,sigma_logD_species);
param_sumBn_species ~ normal(0,sigma_sumBn_species);
intercept_species ~ normal(0,sigma_inter_species);
### plot random param
intercept_plot ~ normal(0,sigma_inter_plot);
### tree random effect
intercept_tree ~ normal(0,sigma_inter_tree);
### ecocode random param
param_Tf_ecocode ~ normal(0, sigma_Tf_ecocode);
param_sumBn_ecocode ~ normal(0, sigma_sumBn_ecocode);
param_sumTfBn_ecocode ~ normal(0 ,sigma_sumTfBn_ecocode);
param_sumTnBn_ecocode ~ normal(0 ,sigma_sumTnBn_ecocode);
param_sumTnTfBn_abs_ecocode ~ normal(0 ,sigma_sumTnTfBn_abs_ecocode);
intercept_set ~ normal(0,sigma_inter_ecocode);
###############################################
########### Non-Hierarchical parameters ########
# constant for prior
sigma0 <- 10;
# slope and intercept
intercept ~ normal(0,sigma0);
mean_logD ~ normal(0,sigma0);
mean_Tf ~ normal(0,sigma0);
mean_sumBn ~ normal(0,sigma0);
mean_sumTfBn ~ normal(0,sigma0);
mean_sumTnBn ~ normal(0,sigma0);
# sigma
sigma_inter_species~ uniform(0,6);
sigma_inter_ecocode~ uniform(0,6);
sigma_inter_plot~ uniform(0,6);
sigma_inter_tree~ uniform(0,6);
sigma_logD_species~ uniform(0,6);
sigma_sumBn_species~ uniform(0,6);
sigma_Tf_ecocode~ uniform(0,6);
sigma_sumBn_ecocode~ uniform(0,6);
sigma_sumTfBn_ecocode~ uniform(0,6);
sigma_sumTnBn_ecocode~ uniform(0,6);
sigma_sumTnTfBn_abs_ecocode~ uniform(0,6);
sigma~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
"
# load data
source("R/analysis/lmer.run.R")
trait <- "SLA"
type.filling <- 'species'
df.lmer <- load.data.for.lmer(trait,
cat.TF = FALSE,
sample.size = 10000) # return a DF
df.lmer <- df.lmer[!is.na(df.lmer$species.id), ]
df.lmer.t <- df.lmer
stan.dat <- list(N_indiv=nrow(df.lmer.t),
N_species=nlevels(factor(df.lmer.t$species.id)),
N_plot=nlevels(factor(df.lmer.t$plot.id)),
N_ecocode=nlevels(factor(df.lmer.t$set.id)),
species_id=unclass(factor(df.lmer.t$species.id)),
plot_id=unclass(factor(df.lmer.t$plot.id)),
ecocode_id=unclass(factor(df.lmer.t$set.id)),
logG=df.lmer.t$logG,
logD=df.lmer.t$logD,
Tf=df.lmer.t$Tf,
sumBn=df.lmer.t$sumBn,
sumTfBn=df.lmer.t$sumTfBn,
sumTnBn=df.lmer.t$sumTnBn)
# param to save
param.to.save <- c('intercept','mean_logD','mean_Tf','mean_sumBn',
'mean_sumTfBn','mean_sumTnBn','sigma_inter_sp',
'sigma_inter_ecocode','sigma_inter_plot','sigma_logD_sp',
'sigma_Tf_ecocode','sigma_sumBn_ecocode',
'sigma_sumTfBn_ecocode','sigma_sumTnBn_ecocode',
'sigma',
'param_Tf_ecocode',
'param_sumBn_ecocode','param_sumTfBn_ecocode',
'param_sumTnBn_ecocode','intercept_ecocode')
#'intercept_sp','param_logD_sp',
#warmup = burn-in
system.time(stan.out <- stan(
model_code=stan.model,
data=stan.dat,
pars=param.to.save,