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

progress on optimizing JAGS and STAN

parent 32c60103
......@@ -8,9 +8,13 @@ source('R/analysis/stan.run.R')
### FUNCTION TO GENERATE THE GOOD FORMAT OF DATA
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.fixed.biomes.R")
model.files.jags.Tf.2 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.biomes.R")
model.files.jags.Tf.3 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.ecocode.R")
......@@ -23,16 +27,16 @@ fun.call.jags.and.save <- function(jags.model,
chains = 3,
thin = 50){
start <- Sys.time()
model.file <-
model.file <-
cat(jags.model$bug, file = file.path(path.out,"model.file.bug")
, sep=" ", fill = FALSE, labels = NULL, append = FALSE)
### SEND to jags
jags.output <- jags(data=list.jags,
jags.output <- jags(data=list.jags,
model.file = file.path(path.out,"model.file.bug"),
parameters.to.save = jags.model$pars,
parameters.to.save = jags.model$pars,
n.chains = chains ,
DIC = FALSE,
DIC = FALSE,
n.iter = iter,
n.burnin = warmup,
n.thin = thin)
......@@ -47,7 +51,7 @@ run.jags <- function (model.file, trait, init.TF,
min.obs = 10, sample.size = NA,
type.filling='species', cat.TF,
fname = 'data.all.no.log.rds',
var.sample = 'ecocode.id',
var.sample = 'ecocode.id',
...) {
require(R2jags)
source(model.file, local = TRUE)
......
load.model <- function(){
list(
name = "LOGLIN.ER.AD.Tf",
## 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'),
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*Tf[i]
SUMBBN[i] <-(mean_sumBn + param_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- mean_sumTfBn *sumTfBn[i]
SUMTNBBN[i] <- mean_sumTnBn *sumTnBn[i]
SUMABSBBN[i] <- mean_sumTnTfBn_abs *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)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
mean_Tf ~ dnorm(0,tau0)T(-5, 5)
mean_sumBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTfBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnTfBn_abs ~ dnorm(0,tau0)T(-5, 5)
# 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,3)
tau_inter_tree <- pow(sigma_inter_tree,-2)
sigma_inter_tree ~ dunif(0.00001,3)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.00001,3)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.00001,3)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.00001,3)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.00001,3)
} # End of the jags model
")
}
load.model <- function(){
list(
name = "LOGLIN.ER.AD.Tf.fixed.biomes",
## 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_sumBn_species',
'sigma', 'param_Tf_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumTnBn_biomes', 'param_sumTnTfBn_abs_biomes'),
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] <- (param_Tf_biomes[biomes_id[i]])*Tf[i]
SUMBBN[i] <-(param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
SUMTNBBN[i] <- (param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
SUMABSBBN[i] <- (param_sumTnTfBn_abs_biomes[biomes_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)T(-10, 10)
param_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)T(-10, 10)
intercept_species[n] ~ dnorm(0,tau_inter_species)T(-10, 10)
}
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)T(-10, 10)
}
## tree effect
for (j in 1:N_tree)
{
intercept_tree[j] ~ dnorm(0,tau_inter_tree)T(-10, 10)
}
## set effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)T(-10, 10)
}
## biomes
for (j in 1:N_biomes)
{
param_Tf_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
param_sumBn_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
param_sumTfBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
param_sumTnBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
param_sumTnTfBn_abs_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau_inter_plot <- pow(sigma_inter_plot,-2)
sigma_inter_plot ~ dunif(0.0001,5)
tau_inter_tree <- pow(sigma_inter_tree,-2)
sigma_inter_tree ~ dunif(0.0001,5)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.0001,5)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.0001,5)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.0001,5)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.0001,5)
#transformed variables
mean_Tf <- mean(param_Tf_biomes)
mean_sumBn <- mean(param_sumBn_biomes)
mean_sumTfBn <- mean(param_sumTfBn_biomes)
mean_sumTnBn <- mean(param_sumTnBn_biomes)
mean_sumTnTfBn_abs <- mean(param_sumTnTfBn_abs_biomes)
} # End of the jags model
")
}
......@@ -14,7 +14,7 @@ load.model <- function(){
'sigma_sumBn_species', 'sigma_sumTfBn_biomes',
'sigma_sumTnBn_biomes', 'sigma_sumTnTfBn_abs_biomes',
'sigma', 'param_Tf_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumTnBn_biomes', 'param_sumTnTfBn_abs_biomes'),
bug =
"#######################################################
......@@ -31,7 +31,7 @@ TTf[i] <- (mean_Tf + param_Tf_biomes[biomes_id[i]])*Tf[i]
SUMBBN[i] <-(mean_sumBn + param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (mean_sumTfBn + param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
SUMTNBBN[i] <- (mean_sumTnBn + param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
SUMABSBBN[i] <- (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
SUMABSBBN[i] <- (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
}
################################################
......@@ -48,28 +48,28 @@ for (n in 1:N_species)
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_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)
intercept_tree[j] ~ dnorm(0,tau_inter_tree)
}
## tree effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
## biomes
## biomes
for (j in 1:N_biomes)
{
param_Tf_biomes[j] ~ dnorm(0, tau_Tf_biomes)
param_sumBn_biomes[j] ~ dnorm(0, tau_sumBn_biomes)
param_sumTfBn_biomes[j] ~ dnorm(0 ,tau_sumTfBn_biomes)
param_sumTnBn_biomes[j] ~ dnorm(0 ,tau_sumTnBn_biomes)
param_sumTnTfBn_abs_biomes[j] ~ dnorm(0 ,tau_sumTnTfBn_abs_biomes)
param_Tf_biomes[j] ~ dnorm(0, tau_Tf_biomes)T(-5, 5)
param_sumBn_biomes[j] ~ dnorm(0, tau_sumBn_biomes)T(-5, 5)
param_sumTfBn_biomes[j] ~ dnorm(0 ,tau_sumTfBn_biomes)T(-5, 5)
param_sumTnBn_biomes[j] ~ dnorm(0 ,tau_sumTnBn_biomes)T(-5, 5)
param_sumTnTfBn_abs_biomes[j] ~ dnorm(0 ,tau_sumTnTfBn_abs_biomes)T(-5, 5)
}
###############################################
########### Non-hierarchical parameters ########
......@@ -77,40 +77,40 @@ for (j in 1:N_biomes)
# 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)
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
mean_Tf ~ dnorm(0,tau0)T(-5, 5)
mean_sumBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTfBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnTfBn_abs ~ dnorm(0,tau0)T(-5, 5)
# 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)
sigma_inter_plot ~ dunif(0.00001,3)
tau_inter_tree <- pow(sigma_inter_tree,-2)
sigma_inter_tree ~ dunif(0.00001,5)
sigma_inter_tree ~ dunif(0.00001,3)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.00001,5)
sigma_inter_set ~ dunif(0.00001,3)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.00001,5)
sigma_inter_species ~ dunif(0.00001,3)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.00001,5)
sigma_logD_species ~ dunif(0.00001,3)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.00001,5)
sigma_sumBn_species ~ dunif(0.00001,3)
tau_Tf_biomes <- pow(sigma_Tf_biomes,-2)
sigma_Tf_biomes ~ dunif(0.00001,5)
sigma_Tf_biomes ~ dunif(0.00001,3)
tau_sumBn_biomes <- pow(sigma_sumBn_biomes,-2)
sigma_sumBn_biomes ~ dunif(0.00001,5)
sigma_sumBn_biomes ~ dunif(0.00001,3)
tau_sumTfBn_biomes <- pow(sigma_sumTfBn_biomes,-2)
sigma_sumTfBn_biomes ~ dunif(0.00001,5)
sigma_sumTfBn_biomes ~ dunif(0.00001,3)
tau_sumTnBn_biomes <- pow(sigma_sumTnBn_biomes,-2)
sigma_sumTnBn_biomes ~ dunif(0.00001,5)
sigma_sumTnBn_biomes ~ dunif(0.00001,3)
tau_sumTnTfBn_abs_biomes <- pow(sigma_sumTnTfBn_abs_biomes,-2)
sigma_sumTnTfBn_abs_biomes ~ dunif(0.00001,5)
sigma_sumTnTfBn_abs_biomes ~ dunif(0.00001,3)
} # End of the jags model
......
......@@ -14,7 +14,7 @@ load.model <- function(){
'sigma_sumBn_species', 'sigma_sumTfBn_ecocode',
'sigma_sumTnBn_ecocode', 'sigma_sumTnTfBn_abs_ecocode',
'sigma', 'param_Tf_ecocode',
'param_sumBn_ecocode', 'param_sumTfBn_ecocode',
'param_sumBn_ecocode', 'param_sumTfBn_ecocode',
'param_sumTnBn_ecocode', 'param_sumTnTfBn_abs_ecocode'),
bug =
"#######################################################
......@@ -31,7 +31,7 @@ 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]
SUMABSBBN[i] <- (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_ecocode[ecocode_id[i]])*sumTnTfBn_abs[i]
}
################################################
......@@ -48,28 +48,28 @@ for (n in 1:N_species)
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_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)
intercept_tree[j] ~ dnorm(0,tau_inter_tree)
}
## tree effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
## ecocode
## 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)
param_Tf_ecocode[j] ~ dnorm(0, tau_Tf_ecocode)T(-5, 5)
param_sumBn_ecocode[j] ~ dnorm(0, tau_sumBn_ecocode)T(-5, 5)
param_sumTfBn_ecocode[j] ~ dnorm(0 ,tau_sumTfBn_ecocode)T(-5, 5)
param_sumTnBn_ecocode[j] ~ dnorm(0 ,tau_sumTnBn_ecocode)T(-5, 5)
param_sumTnTfBn_abs_ecocode[j] ~ dnorm(0 ,tau_sumTnTfBn_abs_ecocode)T(-5, 5)
}
###############################################
########### Non-hierarchical parameters ########
......@@ -77,40 +77,40 @@ for (j in 1:N_ecocode)
# 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)
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
mean_Tf ~ dnorm(0,tau0)T(-5, 5)
mean_sumBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTfBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnBn ~ dnorm(0,tau0)T(-5, 5)
mean_sumTnTfBn_abs ~ dnorm(0,tau0)T(-5, 5)
# 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)
sigma_inter_plot ~ dunif(0.00001,4)
tau_inter_tree <- pow(sigma_inter_tree,-2)
sigma_inter_tree ~ dunif(0.00001,5)
sigma_inter_tree ~ dunif(0.00001,4)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.00001,5)
sigma_inter_set ~ dunif(0.00001,4)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.00001,5)
sigma_inter_species ~ dunif(0.00001,4)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.00001,5)
sigma_logD_species ~ dunif(0.00001,4)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.00001,5)
sigma_sumBn_species ~ dunif(0.00001,4)
tau_Tf_ecocode <- pow(sigma_Tf_ecocode,-2)
sigma_Tf_ecocode ~ dunif(0.00001,5)
sigma_Tf_ecocode ~ dunif(0.00001,4)
tau_sumBn_ecocode <- pow(sigma_sumBn_ecocode,-2)
sigma_sumBn_ecocode ~ dunif(0.00001,5)
sigma_sumBn_ecocode ~ dunif(0.00001,4)
tau_sumTfBn_ecocode <- pow(sigma_sumTfBn_ecocode,-2)
sigma_sumTfBn_ecocode ~ dunif(0.00001,5)
sigma_sumTfBn_ecocode ~ dunif(0.00001,4)
tau_sumTnBn_ecocode <- pow(sigma_sumTnBn_ecocode,-2)
sigma_sumTnBn_ecocode ~ dunif(0.00001,5)
sigma_sumTnBn_ecocode ~ dunif(0.00001,4)
tau_sumTnTfBn_abs_ecocode <- pow(sigma_sumTnTfBn_abs_ecocode,-2)
sigma_sumTnTfBn_abs_ecocode ~ dunif(0.00001,5)
sigma_sumTnTfBn_abs_ecocode ~ dunif(0.00001,4)
} # End of the jags model
......
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.fixed.biomes",
pars = c('intercept' , 'mean_logD', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_inter_tree', 'sigma_logD_species',
'sigma_sumBn_species', '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=-5,upper=5> intercept;
real<lower=-5,upper=5> mean_logD;
real<lower=0.0001,upper=3> sigma_inter_species;
real<lower=0.0001,upper=3> sigma_inter_set;
real<lower=0.0001,upper=3> sigma_inter_plot;
real<lower=0.0001,upper=3> sigma_inter_tree;
real<lower=0.0001,upper=3> sigma_logD_species;
real<lower=0.0001,upper=3> sigma_sumBn_species;
real<lower=0.0001,upper=3> sigma;
real<lower=-5,upper=5> param_Tf_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_species[N_species];
real<lower=-5,upper=5> param_sumTfBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnTfBn_abs_biomes[N_biomes];
real<lower=-5,upper=5> param_logD_species[N_species];
real<lower=-5,upper=5> intercept_species[N_species];
real<lower=-5,upper=5> intercept_plot[N_plot];
real<lower=-5,upper=5> intercept_set[N_set];
real<lower=-5,upper=5> intercept_tree[N_tree];
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
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_species[species_id[i]])*logD[i]
+ (param_Tf_biomes[biomes_id[i]])*Tf[i]
+ ( param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
+ (param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
+ (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,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);
intercept_set ~ normal(0, sigma_inter_set);
### biomes fixed param
param_Tf_biomes ~ normal(0, sigma0);
param_sumBn_biomes ~ normal(0,sigma0);
param_sumTfBn_biomes ~ normal(0 ,sigma0);
param_sumTnBn_biomes ~ normal(0 ,sigma0);
param_sumTnTfBn_abs_biomes ~ normal(0 , sigma0);
###############################################
########### 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~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
}
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.biomes",
pars = c('intercept' , 'mean_logD', 'mean_Tf',
'mean_sumBn', 'mean_sumTfBn', 'mean_sumTnBn',