Current state of the model - vvandermeersch/climategrowthshifts GitHub Wiki

1. Baseline growth

Short description Name Mathematical form Domain Prior model
Global intercept mu_alpha Constant Global mu_alpha ~ normal(0, log(10)/2.32);
Stand intercepts alpha_stand Constant Stand sigma_alpha_stand ~ normal(0, log(10)/2.32); alpha_stand ~ normal(mu_alpha, sigma_alpha_stand);

2. GP component for short-term trajectories (canopy dynamics, microenvironment etc.)

Short description Name Mathematical form Domain Prior model
Length scale rho_short Constant Global (was Species before) rho_short ~ lognormal(log(2.5), 0.3);
Marginal variation log_gamma_short Constant Global (was Species before) log_gamma_short ~ normal(log(0.22), 0.58);

3. GP component for long-term trajectories (allometry, etc.)

Short description Name Mathematical form Domain Prior model
Length scale log_delta_rho Constant Species mu_log_delta_rho ~ normal(log(30), 0.5); sigma_log_delta_rho ~ normal(0, 0.5); log_delta_rho ~ normal(mu_log_delta_rho, sigma_log_delta_rho); rho_long = rho_short + exp(log_delta_rho);
Marginal variation log_gamma_long Constant Species mu_log_gamma_long ~ normal(log(0.42), 0.58); sigma_log_gamma_long ~ normal(0, 0.5); log_gamma_long ~ normal(mu_log_gamma_long, sigma_log_gamma_long);
Individual functional behavior f_ind GP Tree f_ind ~ multi_normal_cholesky(0, L_cov_ind[sp]);

The two GPs are not marginalized (because of missing rings). Each tree gets its own functional behavior, with a non-centered parametrization, i.e. something like:

f_tilde_ind ~ normal(0, 1);
array[N_species] matrix[N_all_years, N_all_years] L_cov_ind;
  {
    for(sp in 1:N_species){
      matrix[N_all_years, N_all_years] cov_ind
      = gp_exp_quad_cov(all_years, gamma_short, rho_short)
      + gp_exp_quad_cov(all_years, gamma_long[sp], rho_long[sp])
      + diag_matrix(rep_vector(1e-8, N_all_years));
      
      L_cov_ind[sp] = cholesky_decompose(cov_ind);
    }
  }
f_ind = L_cov_ind[sp] * f_tilde_ind;

4. Climate regression

Short description Name Mathematical form Domain Prior model
GDD slope beta_gdd Constant Species mu_beta_gdd ~ normal(0, log(1.3) / 2.57); sigma_beta_gdd ~ normal(0, log(1.3) / 2.57); beta_gdd ~ normal(mu_beta_gdd, sigma_beta_gdd);
Winter prec. slope beta_pre Constant Species mu_beta_pre ~ normal(0, log(1.3) / 2.57); sigma_beta_pre ~ normal(0, log(1.3) / 2.57); beta_pre ~ normal(mu_beta_pre, sigma_beta_pre);
VPD slope beta_vpd Constant Species mu_beta_vpd ~ normal(0, log(1.3) / 2.57); sigma_beta_vpd ~ normal(0, log(1.3) / 2.57); beta_vpd ~ normal(mu_beta_vpd, sigma_beta_vpd);
Climate regression ... Linear combination Species x Stand (scale of climate) x Year beta_gdd[sp] * (gdd_obs[stand,year] - gdd0) + beta_pre[sp] * (pre_obs[stand,year] - pre0) + beta_vpd[sp] * (vpd_obs[stand,year] - vpd0)

5. "Missing" mesoenvironment

Short description Name Mathematical form Domain Prior model
Amplitude tau_clim Constant Global tau_clim ~ normal(0, log(2)/2.57);
Autoregressive parameter rho_clim Constant Global rho_prim_clim ~ beta(4, 2); rho_clim = rho_prim_clim*2-1;
Intercepts delta_clim Constant Stand x Years delta_clim[s][1] = 0; delta_clim[s][i] = rho_clim * delta_clim[s][i-1] + tau_clim * delta_tilde_clim[s][i-1];
Species sensitivity kappa_clim Constant Species log_kappa_clim ~ normal(0, 0.3 / 2.57); kappa_clim = exp(log_kappa_clim);

The stand x years intercepts are scaled by kappa_clim[sp], as some stands have several species that might not respond in a similar way to these variations. I also defined sum_to_zero_vector[N_species] log_kappa_clim; as a sum_to_zero vector, so the baseline is set to a virtual species of reference with kappa_clim = 1. Maybe it was a bad idea, and I could also pick one of the species of the dataset as a reference and set kappa_clim = 1 for this species (i.e. estimating only [N_species-1] log_kappa_clim;).

5. Idiosyncratic abrupt deviations (and missing rings)

Short description Name Mathematical form Domain Prior model
Amplitude tau_idio Constant Global tau_idio ~ normal(0, log(5) / 2.57);
Simplex thetas_idio Constant Tree thetas_baseline = [100, 20, 1]'; omega_thetas = 10; alphas = thetas_baseline / omega_thetas + ones; thetas_idio[t] ~ dirichlet(alphas);

thetas_idio[t][2] is the probability of an abrupt deviation ("idiosyncratic shock"), thetas_idio[t][3] is the probability of a missing ring ("idiosyncratic shutdown").

5. Concordant shocks (and shutdowns)

5.1. Stand-level extreme event
Short description Name Mathematical form Domain Prior model
Baseline probability of stand-level event (intercept) logit_phi_sck Constant Stand mu_phi ~ normal(-2.17, 0.35); sigma_phi ~ normal(0, 0.7); logit_phi_sck ~ normal(mu_phi, sigma_phi);
VPD slope beta_phi_vpd Constant Global beta_phi_vpd ~ normal(0, 0.3);
Winter precipitation slope beta_phi_pre Constant Global beta_phi_pre ~ normal(0, 0.3);
Stand-level event probability phi_sck Linear combination Stand x Year phi_sck = inv_logit(mu_phi + logit_phi_sck[st] + beta_phi_vpd * (vpd_obs[stand,year] - vpd0) + beta_phi_pre * (pre_obs[stand,year] - pre0));

Oct. 8: I spotted a huge mistake. mu_phi is applied both in logit_phi_sck ~ normal(mu_phi, sigma_phi); AND in phi_sck = inv_logit(mu_phi + ...)... I will run the correct model again and finish this...

5.2. Tree-level concordant response
Short description Name Mathematical form Domain Prior model
Shock amplitude tau_conc Constant Species mu_log_tau_conc ~ normal(log(1), 0.5); sigma_log_tau_conc ~ normal(0, 0.5); tau_conc ~ lognormal(mu_log_tau_conc, sigma_log_tau_conc);
Baseline concordance mu_omega_conc Constant Global mu_omega_conc ~ normal(-1.5, 0.28);
Species-level concordance alpha_omega_species Constant Species sigma_omega_species ~ normal(0, 0.45); alpha_omega_species ~ normal(0, sigma_omega_species);
Stand-level concordance alpha_omega_stand Constant Stand sigma_omega_stand ~ normal(0, 0.45); alpha_omega_stand ~ normal(0, sigma_omega_stand);
Probability of concordant response (given a stand-level event) omega_conc_sck Constant Stand x Species omega_conc_sck = inv_logit(mu_omega_conc + alpha_omega_species[sp] + alpha_omega_stand[st]);