data { int N_CONDITION; int N_TREATMENT; int N; int disease_condition_seg[N_CONDITION+1]; int treatment_seg[N_TREATMENT+1]; vector[N] concentration; vector[N] y; } parameters { real mu_overarch_0; real mu_overarch_1; vector[N_CONDITION] sigma_overarch_0; vector[N_CONDITION] sigma_overarch_1; // vector[N_CONDITION] nu_overarch_0; // vector[N_CONDITION] nu_overarch_1; vector[N_TREATMENT] sigma_treatment_0; vector[N_TREATMENT] sigma_treatment_1; // vector[N_TREATMENT] nu_treatment_0; // vector[N_TREATMENT] nu_treatment_1; vector[N_CONDITION] mu_condition_0; vector[N_CONDITION] mu_condition_1; vector[N_TREATMENT] mu_observation_0; vector[N_TREATMENT] mu_observation_1; } transformed parameters { vector[N_TREATMENT] mu_treat_0; vector[N_TREATMENT] mu_treat_1; vector[N] mu_obs_0; vector[N] mu_obs_1; vector[N] mu_obs; // vector[N] nu_obs; vector[N] sigma_obs; for (i in 2:N_CONDITION+1) { if ((disease_condition_seg[i]-disease_condition_seg[i-1]) == 1) { mu_treat_0[disease_condition_seg[i]] = mu_condition_0[i-1]; mu_treat_1[disease_condition_seg[i]] = mu_condition_1[i-1]; } else { mu_treat_0[(disease_condition_seg[i-1]+1):disease_condition_seg[i]] = rep_vector(mu_condition_0[i], disease_condition_seg[i]-disease_condition_seg[i-1]); mu_treat_1[(disease_condition_seg[i-1]+1):disease_condition_seg[i]] = rep_vector(mu_condition_1[i], disease_condition_seg[i]-disease_condition_seg[i-1]); } } // mu = beta_0 + beta_1 * concentration for (i in 2:N_TREATMENT+1) { // print(treatment_seg[i]); // print(treatment_seg[i-1]); if ((treatment_seg[i] - treatment_seg[i-1]) == 1) { mu_obs_0[treatment_seg[i]] = mu_observation_0[i-1]; mu_obs_1[treatment_seg[i]] = mu_observation_1[i-1]; // nu_obs[treatment_seg[i]] = nu_treatment_0[i-1]; sigma_obs[treatment_seg[i]] = sigma_treatment_0[i-1]; } else { mu_obs_0[treatment_seg[i-1]+1:treatment_seg[i]] = rep_vector(mu_observation_0[i-1], treatment_seg[i] - treatment_seg[i-1]); mu_obs_1[treatment_seg[i-1]+1:treatment_seg[i]] = rep_vector(mu_observation_1[i-1], treatment_seg[i] - treatment_seg[i-1]); // nu_obs[treatment_seg[i-1]+1:treatment_seg[i]] = rep_vector(nu_treatment_0[i-1], treatment_seg[i] - treatment_seg[i-1]); sigma_obs[treatment_seg[i-1]+1:treatment_seg[i]] = rep_vector(sigma_treatment_0[i-1], treatment_seg[i] - treatment_seg[i-1]); } } mu_obs = mu_obs_0 + mu_obs_1 .* concentration; } model { // mu_overarch_0 ~ normal(0, 1); target += normal_lpdf(mu_overarch_0 | 0, 1); // mu_overarch_1 ~ normal(0, 1); target += normal_lpdf(mu_overarch_1 | 0, 1); // //sigma_overarch_0 ~ inv_gamma(1, 1); // target += inv_gamma_lpdf(sigma_overarch_0 | 1, 1); // // sigma_overarch_1 ~ inv_gamma(1, 1); // target += inv_gamma_lpdf(sigma_overarch_1 | 1, 1); // // sigma_treatment_0 ~ inv_gamma(1, 1); // target += inv_gamma_lpdf(sigma_treatment_0 | 1, 1); // // sigma_treatment_1 ~ inv_gamma(1, 1); // target += inv_gamma_lpdf(sigma_treatment_1 | 1, 1); // beta_0 // mu_condition_0 ~ normal(mu_overarch_0, sigma_overarch_0); target += normal_lpdf(mu_condition_0 | mu_overarch_0, sigma_overarch_0); // mu_observation_0 ~ normal(mu_treat_0, sigma_treatment_0); target += normal_lpdf(mu_observation_0 | mu_treat_0, sigma_treatment_0); // beta_1 // mu_condition_1 ~ normal(mu_overarch_1, sigma_overarch_1); target += normal_lpdf(mu_condition_1 | mu_overarch_1, sigma_overarch_1); // mu_observation_1 ~ normal(mu_treat_1, sigma_treatment_1); target += normal_lpdf(mu_observation_1 | mu_treat_1, sigma_treatment_1); // beta_0 + beta_1 * concentration // y ~ normal(mu_obs, sigma_obs); target += normal_lpdf(y | mu_obs, sigma_obs); }