// version: OCT-07, 2019
// Programmer: Yuxin Chen (chenyux9@gmail.com)

data {
int<lower = 0> N; // number of observations (plots across years)
int<lower = 0> N_plot;
int<lower = 0> N_comm; // number of community composition
int<lower = 0> N_div; // number of diversity levels
int<lower = 0> N_age;
vector<lower = 0>[N] y; // stand volume
vector[N] site;
vector[N] log_SR;
int comm[N]; // ID of community composition
int plot[N]; // ID of plot
int age_level[N]; // ID of forest age
vector[N_age] age_level_c; // standardized forest age
}

parameters {
real<lower = 0> sd_0;
real<lower = 0> sd_2;
real<lower = 0> sd_comm; // SD for random effects of community composition
real<lower = 0> sd_plot; // SD for random effects of plot
vector<lower = 0>[N_age] sd_y;
real alpha1;
real<lower = 0> beta0_0;
real beta0_1;
real beta2_0;
real beta2_1;
vector[N_comm] comm_raw;
vector[N_plot] plot_raw;
vector<lower = 0>[N_age] alpha0;
vector[N_age] alpha2;
}

transformed parameters {
vector[N_plot] r_plot; // random effects of plot
vector[N_comm] r_comm; // random effects of community composition
for (i in 1:N_plot){
r_plot[i] = sd_plot * plot_raw[i];
}
for (j in 1:N_comm){
r_comm[j] = sd_comm * comm_raw[j];
}
}

model {
alpha1 ~ normal(0, 5);
beta0_0 ~ normal(0, 5);
beta0_1 ~ normal(0, 5);
beta2_0 ~ normal(0, 5);
beta2_1 ~ normal(0, 5);
comm_raw ~ normal(0, 1);
plot_raw ~ normal(0, 1);
sd_0 ~ cauchy(0, 5);
sd_2 ~ cauchy(0, 5);
sd_comm ~ cauchy(0, 5);
sd_plot ~ cauchy(0, 5);
sd_y ~ cauchy(0, 5);

alpha0 ~ normal(beta0_0 + beta0_1 * age_level_c, sd_0);
alpha2 ~ normal(beta2_0 + beta2_1 * age_level_c, sd_2);

// main-model 
for (i in 1:N){
y[i] ~ normal(alpha0[age_level[i]] + alpha1 * site[i] + alpha2[age_level[i]] * log_SR[i] +
              r_comm[comm[i]] + r_plot[plot[i]], sd_y[age_level[i]]);
}
}