Commit 32c60103 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

include fixed effect biomes code ready to be sent for jags and comparison between stan and jags

parent 1ae052a9
...@@ -10,16 +10,18 @@ source('R/analysis/stan.run.R') ...@@ -10,16 +10,18 @@ source('R/analysis/stan.run.R')
model.files.jags.Tf.1 <- model.files.jags.Tf.1 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.biomes.R") c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.biomes.R")
model.files.jags.Tf.2 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.ecocode.R")
fun.call.jags.and.save <- function(jags.model, fun.call.jags.and.save <- function(jags.model,
list.jags, list.jags,
path.out, path.out,
var.sample, var.sample,
iter = 100, iter = 50000,
warmup = 50, warmup = 5000,
chains = 3, chains = 3,
thin = 5){ thin = 50){
start <- Sys.time() start <- Sys.time()
model.file <- model.file <-
cat(jags.model$bug, file = file.path(path.out,"model.file.bug") cat(jags.model$bug, file = file.path(path.out,"model.file.bug")
......
...@@ -27,6 +27,10 @@ run.model.for.set.one.trait <- function(model.file, fun.model, trait, ...@@ -27,6 +27,10 @@ run.model.for.set.one.trait <- function(model.file, fun.model, trait,
#===================================================================== #=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model # Run lmer() (in package lme4) for one ecoregion, one trait, one model
#===================================================================== #=====================================================================
model.files.lmer.Tf.0 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.r.biomes.species.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.r.biomes.species.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.Tf.r.biomes.species.R")
model.files.lmer.Tf.1 <- model.files.lmer.Tf.1 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.species.R") c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.species.R")
...@@ -36,6 +40,8 @@ model.files.lmer.Tf.2 <- ...@@ -36,6 +40,8 @@ model.files.lmer.Tf.2 <-
model.files.lmer.Tf.3 <- model.files.lmer.Tf.3 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.biomes.species.R") c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.biomes.species.R")
model.files.lmer.Tf.3b <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species.R")
model.files.lmer.Tf.4 <- model.files.lmer.Tf.4 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R") c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R")
...@@ -175,8 +181,7 @@ load.and.save.data.for.lmer <- function(trait, ...@@ -175,8 +181,7 @@ load.and.save.data.for.lmer <- function(trait,
min.obs= 10, min.obs= 10,
type.filling = species, type.filling = species,
fname = 'data.all.no.log.rds', fname = 'data.all.no.log.rds',
base.dir = "output/processed", base.dir = "output/processed"){
data.table.TF = FALSE){
data.tree.tot <- readRDS(file.path(base.dir, fname)) data.tree.tot <- readRDS(file.path(base.dir, fname))
df <- fun.select.data.for.lmer(data.tree.tot, trait, df <- fun.select.data.for.lmer(data.tree.tot, trait,
type.filling = type.filling) type.filling = type.filling)
...@@ -190,7 +195,7 @@ fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp, ...@@ -190,7 +195,7 @@ fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp,
perc.neigh.gs, BATOT, min.obs, perc.neigh.gs, BATOT, min.obs,
min.perc.neigh.sp = 0.90, min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95, min.perc.neigh.gs = 0.95,
min.BA.G = -60, min.BA.G = -100,
max.BA.G = 500){ max.BA.G = 500){
select.temp <-(1:nrow(data.tree))[!is.na(data.tree[["BA.G"]]) & select.temp <-(1:nrow(data.tree))[!is.na(data.tree[["BA.G"]]) &
...@@ -232,8 +237,8 @@ return(x/sd(x)) ...@@ -232,8 +237,8 @@ return(x/sd(x))
### get variables for lmer ### get variables for lmer
fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn, fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, abs.CWM.tntf, tf,
min.BA.G = 60){ min.BA.G = 100){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]] + min.BA.G)) logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]] + min.BA.G+1))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]])) logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]]) MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]]) MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
...@@ -273,8 +278,8 @@ return(data.frame(logG = logG, ...@@ -273,8 +278,8 @@ return(data.frame(logG = logG,
fun.get.the.variables.for.lmer.tree.id.cat <- function(data.tree, BATOT, CWM.tn, fun.get.the.variables.for.lmer.tree.id.cat <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, abs.CWM.tntf, tf,
min.BA.G = 60){ min.BA.G = 100){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G)) logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G+1))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]])) logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]]) MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]]) MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
......
load.model <- function(){
list(
name = "LOGLIN.ER.AD.Tf.r.ecocode",
## parameters to save
pars = c('intercept' , 'mean_logD', 'mean_Tf',
'mean_sumBn', 'mean_sumTfBn', 'mean_sumTnBn',
'mean_sumTnTfBn_abs', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_inter_tree', 'sigma_logD_species',
'sigma_Tf_ecocode', 'sigma_sumBn_ecocode',
'sigma_sumBn_species', 'sigma_sumTfBn_ecocode',
'sigma_sumTnBn_ecocode', 'sigma_sumTnTfBn_abs_ecocode',
'sigma', 'param_Tf_ecocode',
'param_sumBn_ecocode', 'param_sumTfBn_ecocode',
'param_sumTnBn_ecocode', 'param_sumTnTfBn_abs_ecocode'),
bug =
"#######################################################
######################## growth model with JAGS #######
model {
############ Likelihood ###################
for (i in 1:N_indiv) {
logG[i] ~ dnorm(theo_g[i],tau)
theo_g[i] <- inter[i] + DD[i] +TTf[i] + SUMBBN[i] +SUMTFBBN[i] + SUMTNBBN[i] + SUMABSBBN[i]
inter[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]]
DD[i] <- (mean_logD + param_logD_species[species_id[i]])*logD[i]
TTf[i] <- (mean_Tf + param_Tf_ecocode[ecocode_id[i]])*Tf[i]
SUMBBN[i] <-(mean_sumBn + param_sumBn_ecocode[ecocode_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (mean_sumTfBn + param_sumTfBn_ecocode[ecocode_id[i]])*sumTfBn[i]
SUMTNBBN[i] <- (mean_sumTnBn + param_sumTnBn_ecocode[ecocode_id[i]])*sumTnBn[i]
SUMABSBBN[i] <- (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_ecocode[ecocode_id[i]])*sumTnTfBn_abs[i]
}
################################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N_species)
{
param_logD_species[n] ~ dnorm(0,tau_logD_species)
param_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)
intercept_species[n] ~ dnorm(0,tau_inter_species)
}
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)
}
## tree effect
for (j in 1:N_tree)
{
intercept_tree[j] ~ dnorm(0,tau_inter_tree)
}
## tree effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
## ecocode
for (j in 1:N_ecocode)
{
param_Tf_ecocode[j] ~ dnorm(0, tau_Tf_ecocode)
param_sumBn_ecocode[j] ~ dnorm(0, tau_sumBn_ecocode)
param_sumTfBn_ecocode[j] ~ dnorm(0 ,tau_sumTfBn_ecocode)
param_sumTnBn_ecocode[j] ~ dnorm(0 ,tau_sumTnBn_ecocode)
param_sumTnTfBn_abs_ecocode[j] ~ dnorm(0 ,tau_sumTnTfBn_abs_ecocode)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
intercept ~ dnorm(0,tau0)T(-10, 10)
mean_logD ~ dnorm(0,tau0)T(-10, 10)
mean_Tf ~ dnorm(0,tau0)T(-10, 10)
mean_sumBn ~ dnorm(0,tau0)T(-10, 10)
mean_sumTfBn ~ dnorm(0,tau0)T(-10, 10)
mean_sumTnBn ~ dnorm(0,tau0)T(-10, 10)
mean_sumTnTfBn_abs ~ dnorm(0,tau0)T(-10, 10)
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau_inter_plot <- pow(sigma_inter_plot,-2)
sigma_inter_plot ~ dunif(0.00001,5)
tau_inter_tree <- pow(sigma_inter_tree,-2)
sigma_inter_tree ~ dunif(0.00001,5)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.00001,5)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.00001,5)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.00001,5)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.00001,5)
tau_Tf_ecocode <- pow(sigma_Tf_ecocode,-2)
sigma_Tf_ecocode ~ dunif(0.00001,5)
tau_sumBn_ecocode <- pow(sigma_sumBn_ecocode,-2)
sigma_sumBn_ecocode ~ dunif(0.00001,5)
tau_sumTfBn_ecocode <- pow(sigma_sumTfBn_ecocode,-2)
sigma_sumTfBn_ecocode ~ dunif(0.00001,5)
tau_sumTnBn_ecocode <- pow(sigma_sumTnBn_ecocode,-2)
sigma_sumTnBn_ecocode ~ dunif(0.00001,5)
tau_sumTnTfBn_abs_ecocode <- pow(sigma_sumTnTfBn_abs_ecocode,-2)
sigma_sumTnTfBn_abs_ecocode ~ dunif(0.00001,5)
} # End of the jags model
")
}
load.model <- function () {
list(name="lmer.LOGLIN.AD.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+Tf+Tf:biomes.id+logD+sumBn+sumBn:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)"))
}
load.model <- function(){ load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.biomes", list(name = "LOGLIN.ER.AD.Tf.r.biomes",
...@@ -68,14 +65,14 @@ parameters{ ...@@ -68,14 +65,14 @@ parameters{
vector[N_biomes] param_sumTnBn_biomes; vector[N_biomes] param_sumTnBn_biomes;
vector[N_biomes] param_sumTnTfBn_abs_biomes; vector[N_biomes] param_sumTnTfBn_abs_biomes;
vector[N_species] param_logD_species; vector[N_species] param_logD_species;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
vector[N_species] intercept_species; vector[N_species] intercept_species;
vector[N_plot] intercept_plot; vector[N_plot] intercept_plot;
vector[N_set] intercept_set; vector[N_set] intercept_set;
vector[N_tree] intercept_tree; vector[N_tree] intercept_tree;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){ for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_species[species_id[i]] theo_g[i] <- intercept + intercept_species[species_id[i]]
...@@ -83,8 +80,8 @@ transformed parameters { ...@@ -83,8 +80,8 @@ transformed parameters {
+ intercept_plot[plot_id[i]] + intercept_plot[plot_id[i]]
+ intercept_set[set_id[i]] + intercept_set[set_id[i]]
+ (mean_logD + param_logD_species[species_id[i]])*logD[i] + (mean_logD + param_logD_species[species_id[i]])*logD[i]
+ (mean_Tf + param_Tf_biomes[biomes_id[i]])*Tf[i] + (mean_Tf +param_Tf_biomes[biomes_id[i]])*Tf[i]
+ (mean_sumBn + param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i] + (mean_sumBn + param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (mean_sumTfBn + param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i] + (mean_sumTfBn + param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
+ (mean_sumTnBn + param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i] + (mean_sumTnBn + param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
+ (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i] + (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
...@@ -104,19 +101,19 @@ model { ...@@ -104,19 +101,19 @@ model {
### species random param ### species random param
param_logD_species ~ normal(0,sigma_logD_species); param_logD_species ~ normal(0,sigma_logD_species);
param_sumBn_species ~ normal(0,sigma_sumBn_species); param_sumBn_species ~ normal(0,sigma_sumBn_species);
intercept_species ~ normal(0,sigma_inter_species); intercept_species ~ normal(0, sigma_inter_species);
### plot random param ### plot random param
intercept_plot ~ normal(0,sigma_inter_plot); intercept_plot ~ normal(0,sigma_inter_plot);
### tree random effect ### tree random effect
intercept_tree ~ normal(0,sigma_inter_tree); intercept_tree ~ normal(0, sigma_inter_tree);
### biomes random param intercept_set ~ normal(0, sigma_inter_set);
param_Tf_biomes ~ normal(0, sigma_Tf_biomes); ### biomes random param
param_sumBn_biomes ~ normal(0, sigma_sumBn_biomes); param_Tf_biomes ~ normal(0, sigma_Tf_biomes);
param_sumBn_biomes ~ normal(0,sigma_sumBn_biomes);
param_sumTfBn_biomes ~ normal(0 ,sigma_sumTfBn_biomes); param_sumTfBn_biomes ~ normal(0 ,sigma_sumTfBn_biomes);
param_sumTnBn_biomes ~ normal(0 ,sigma_sumTnBn_biomes); param_sumTnBn_biomes ~ normal(0 ,sigma_sumTnBn_biomes);
param_sumTnTfBn_abs_biomes ~ normal(0 ,sigma_sumTnTfBn_abs_biomes); param_sumTnTfBn_abs_biomes ~ normal(0 , sigma_sumTnTfBn_abs_biomes);
intercept_set ~ normal(0,sigma_inter_set);
############################################### ###############################################
########### Non-Hierarchical parameters ######## ########### Non-Hierarchical parameters ########
# constant for prior # constant for prior
......
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.biomes",
pars = c('intercept' , 'mean_logD', 'mean_Tf',
'mean_sumBn', 'mean_sumTfBn', 'mean_sumTnBn',
'mean_sumTnTfBn_abs', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_inter_tree', 'sigma_logD_species',
'sigma_Tf_biomes', 'sigma_sumBn_biomes',
'sigma_sumBn_species', 'sigma_sumTfBn_biomes',
'sigma_sumTnBn_biomes', 'sigma_sumTnTfBn_abs_biomes',
'sigma', 'param_Tf_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumTnBn_biomes', 'param_sumTnTfBn_abs_biomes'),## parameters to save
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_biomes;
int<lower=0> N_set;
int<lower=0> N_tree;
int species_id[N_indiv];
int plot_id[N_indiv];
int biomes_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_biomes;
real<lower=0.0001,upper=5> sigma_sumBn_biomes;
real<lower=0.0001,upper=5> sigma_sumBn_species;
real<lower=0.0001,upper=5> sigma_sumTfBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_abs_biomes;
real<lower=0.0001,upper=5> sigma;
vector[N_biomes] param_Tf_biomes;
vector[N_biomes] param_sumBn_biomes;
vector[N_species] param_sumBn_species;
vector[N_biomes] param_sumTfBn_biomes;
vector[N_biomes] param_sumTnBn_biomes;
vector[N_biomes] param_sumTnTfBn_abs_biomes;
vector[N_species] param_logD_species;
vector[N_species] intercept_species;
vector[N_plot] intercept_plot;
vector[N_set] intercept_set;
vector[N_tree] intercept_tree;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + sigma_inter_species*intercept_species[species_id[i]]
+ sigma_inter_tree*intercept_tree[tree_id[i]]*tree_01[i]
+ sigma_inter_plot*intercept_plot[plot_id[i]]
+ sigma_inter_set*intercept_set[set_id[i]]
+ (mean_logD + sigma_logD_species*param_logD_species[species_id[i]])*logD[i]
+ (mean_Tf + sigma_Tf_biomes*param_Tf_biomes[biomes_id[i]])*Tf[i]
+ (mean_sumBn + sigma_sumBn_biomes*param_sumBn_biomes[biomes_id[i]] + sigma_sumBn_species*param_sumBn_species[species_id[i]])*sumBn[i]
+ (mean_sumTfBn + sigma_sumTfBn_biomes*param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
+ (mean_sumTnBn + sigma_sumTnBn_biomes*param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
+ (mean_sumTnTfBn_abs + sigma_sumTnTfBn_abs_biomes*param_sumTnTfBn_abs_biomes[biomes_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,1);
param_sumBn_species ~ normal(0,1);
intercept_species ~ normal(0,1);
### plot random param
intercept_plot ~ normal(0,1);
### tree random effect
intercept_tree ~ normal(0,1);
### biomes random param
param_Tf_biomes ~ normal(0, 1);
param_sumBn_biomes ~ normal(0,1);
param_sumTfBn_biomes ~ normal(0 ,1);
param_sumTnBn_biomes ~ normal(0 ,1);
param_sumTnTfBn_abs_biomes ~ normal(0 , 1);
intercept_set ~ normal(0, 1);
###############################################
########### 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_set~ 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_biomes~ uniform(0,6);
sigma_sumBn_biomes~ uniform(0,6);
sigma_sumTfBn_biomes~ uniform(0,6);
sigma_sumTnBn_biomes~ uniform(0,6);
sigma_sumTnTfBn_abs_biomes~ uniform(0,6);
sigma~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
}
...@@ -10,6 +10,7 @@ model.files.stan.Tf.1 <- ...@@ -10,6 +10,7 @@ model.files.stan.Tf.1 <-
c("R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.no.R", c("R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.no.R",
"R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.set.R", "R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.set.R",
"R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.biomes.R", "R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.biomes.R",
"R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.biomes.bis.R",
"R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.ecocode.R") "R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.ecocode.R")
......
### test jags
source('R/analysis/jags.run.R')
run.multiple.model.for.set.one.trait(model.files.jags.Tf.1, run.jags,
trait = 'Wood.density',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 200, warmup = 20, chains = 3,
init.TF = TRUE)
t1000 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.jags.Tf.1,
run.jags,
trait = 'Wood.density',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 1000, warmup = 20, chains = 3,
init.TF = FALSE))
t1000.4 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.jags.Tf.1,
run.jags,
trait = 'Wood.density',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 5000, warmup = 20, chains = 3,
init.TF = FALSE))
t1000.8 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.jags.Tf.1,
run.jags,
trait = 'Wood.density',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 10000, warmup = 20, chains = 3,
init.TF = FALSE, thin = 10))
t1000.4 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[3], run.stan,
trait = 'SLA',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 400, warmup = 20, chains = 2, chain_id = '1',
init.TF = FALSE))
t1000.8 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[3], run.stan,
trait = 'SLA',type.filling='species', sample.size = 1000,
var.sample = 'ecocode.id',
iter = 800, warmup = 20, chains = 2, chain_id = '1',
init.TF = FALSE))
sample.vec <- c(1000, 10000, 50000, 100000)
t <- c(t1000['elapsed'], t10000['elapsed'], t50000['elapsed'], t100000['elapsed'])
t4 <- c(t1000.4['elapsed'], t10000.4['elapsed'], t50000.4['elapsed'], t100000.4['elapsed'])
t8 <- c(t1000.8['elapsed'], t10000.8['elapsed'], t50000.8['elapsed'])
pdf('time.stan.pdf')
plot(sample.vec, t, ylim = range(t, t4, t.b, t4.b), type = 'b')
lines(sample.vec, t4, type = 'b', col = 'green')
lines(sample.vec[1:3], t8, type = 'b', col = 'pink')
dev.off()
###
### test stan ### test stan
## I changed the coding of the random effect with norm(0,10*sigma it seems more slow!
## PROBLEM WITH STANDARDIZED VARIABLE DO I NEED TO INCLUDE AN INTERCEPT??
## VERY BAD CONVERGENCE ???
## Need to compare the two code
source('R/analysis/stan.run.R') source('R/analysis/stan.run.R')
run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[3], run.stan, run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[3], run.stan,