MrP with CCES example code

mrp_cces.R

vote_model <- brm(mvbind(vote_trump,vote_clinton) ~ 1 +
                    # pop level
                    state_clinton + state_black_pct + state_hispanic_pct + state_white_protestant + 
                    # group effects
                    (1|sex) + (1|race) + (1|age) + (1|edu) + (1|state_name) + (1|region) +
                    # interactions b/t group effects
                    (1|race:edu) + (1|age:edu) + (1|race:region),
                  data = cces,
                  family = bernoulli(link = "logit"),
                  # stan stuff
                  iter=1000, chains=1, cores=1, warmup=500,
                  algorithm='sampling',
                  save_model = 'output/stan_model_code.stan', 
                  # priors (https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations)
                  prior = c(set_prior("normal(0, 1)", class = "Intercept"),
                            set_prior("student_t(3, 0, 4)", class = "sd",resp=c('votetrump','voteclinton'))),
                  # help sampling algorithm converge (http://mc-stan.org/misc/warnings.html)
                  control = list(adapt_delta = 0.95,
                                 max_treedepth = 10),
                  seed = 20190605,
                  refresh = 1)


### generated model code
// generated with brms 2.9.0
functions {
}
data {
  int<lower=1> N;  // number of observations
  int<lower=1> N_votetrump;  // number of observations
  int Y_votetrump[N_votetrump];  // response variable
  int<lower=1> K_votetrump;  // number of population-level effects
  matrix[N_votetrump, K_votetrump] X_votetrump;  // population-level design matrix
  int<lower=1> N_voteclinton;  // number of observations
  int Y_voteclinton[N_voteclinton];  // response variable
  int<lower=1> K_voteclinton;  // number of population-level effects
  matrix[N_voteclinton, K_voteclinton] X_voteclinton;  // population-level design matrix
  // data for group-level effects of ID 1
  int<lower=1> N_1;
  int<lower=1> M_1;
  int<lower=1> J_1_votetrump[N_votetrump];
  vector[N_votetrump] Z_1_votetrump_1;
  // data for group-level effects of ID 2
  int<lower=1> N_2;
  int<lower=1> M_2;
  int<lower=1> J_2_votetrump[N_votetrump];
  vector[N_votetrump] Z_2_votetrump_1;
  // data for group-level effects of ID 3
  int<lower=1> N_3;
  int<lower=1> M_3;
  int<lower=1> J_3_votetrump[N_votetrump];
  vector[N_votetrump] Z_3_votetrump_1;
  // data for group-level effects of ID 4
  int<lower=1> N_4;
  int<lower=1> M_4;
  int<lower=1> J_4_votetrump[N_votetrump];
  vector[N_votetrump] Z_4_votetrump_1;
  // data for group-level effects of ID 5
  int<lower=1> N_5;
  int<lower=1> M_5;
  int<lower=1> J_5_votetrump[N_votetrump];
  vector[N_votetrump] Z_5_votetrump_1;
  // data for group-level effects of ID 6
  int<lower=1> N_6;
  int<lower=1> M_6;
  int<lower=1> J_6_votetrump[N_votetrump];
  vector[N_votetrump] Z_6_votetrump_1;
  // data for group-level effects of ID 7
  int<lower=1> N_7;
  int<lower=1> M_7;
  int<lower=1> J_7_votetrump[N_votetrump];
  vector[N_votetrump] Z_7_votetrump_1;
  // data for group-level effects of ID 8
  int<lower=1> N_8;
  int<lower=1> M_8;
  int<lower=1> J_8_votetrump[N_votetrump];
  vector[N_votetrump] Z_8_votetrump_1;
  // data for group-level effects of ID 9
  int<lower=1> N_9;
  int<lower=1> M_9;
  int<lower=1> J_9_votetrump[N_votetrump];
  vector[N_votetrump] Z_9_votetrump_1;
  // data for group-level effects of ID 10
  int<lower=1> N_10;
  int<lower=1> M_10;
  int<lower=1> J_10_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_10_voteclinton_1;
  // data for group-level effects of ID 11
  int<lower=1> N_11;
  int<lower=1> M_11;
  int<lower=1> J_11_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_11_voteclinton_1;
  // data for group-level effects of ID 12
  int<lower=1> N_12;
  int<lower=1> M_12;
  int<lower=1> J_12_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_12_voteclinton_1;
  // data for group-level effects of ID 13
  int<lower=1> N_13;
  int<lower=1> M_13;
  int<lower=1> J_13_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_13_voteclinton_1;
  // data for group-level effects of ID 14
  int<lower=1> N_14;
  int<lower=1> M_14;
  int<lower=1> J_14_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_14_voteclinton_1;
  // data for group-level effects of ID 15
  int<lower=1> N_15;
  int<lower=1> M_15;
  int<lower=1> J_15_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_15_voteclinton_1;
  // data for group-level effects of ID 16
  int<lower=1> N_16;
  int<lower=1> M_16;
  int<lower=1> J_16_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_16_voteclinton_1;
  // data for group-level effects of ID 17
  int<lower=1> N_17;
  int<lower=1> M_17;
  int<lower=1> J_17_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_17_voteclinton_1;
  // data for group-level effects of ID 18
  int<lower=1> N_18;
  int<lower=1> M_18;
  int<lower=1> J_18_voteclinton[N_voteclinton];
  vector[N_voteclinton] Z_18_voteclinton_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc_votetrump = K_votetrump - 1;
  matrix[N_votetrump, Kc_votetrump] Xc_votetrump;  // centered version of X_votetrump
  vector[Kc_votetrump] means_X_votetrump;  // column means of X_votetrump before centering
  int Kc_voteclinton = K_voteclinton - 1;
  matrix[N_voteclinton, Kc_voteclinton] Xc_voteclinton;  // centered version of X_voteclinton
  vector[Kc_voteclinton] means_X_voteclinton;  // column means of X_voteclinton before centering
  for (i in 2:K_votetrump) {
    means_X_votetrump[i - 1] = mean(X_votetrump[, i]);
    Xc_votetrump[, i - 1] = X_votetrump[, i] - means_X_votetrump[i - 1];
  }
  for (i in 2:K_voteclinton) {
    means_X_voteclinton[i - 1] = mean(X_voteclinton[, i]);
    Xc_voteclinton[, i - 1] = X_voteclinton[, i] - means_X_voteclinton[i - 1];
  }
}
parameters {
  vector[Kc_votetrump] b_votetrump;  // population-level effects
  real temp_votetrump_Intercept;  // temporary intercept
  vector[Kc_voteclinton] b_voteclinton;  // population-level effects
  real temp_voteclinton_Intercept;  // temporary intercept
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  vector[N_1] z_1[M_1];  // unscaled group-level effects
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations
  vector[N_2] z_2[M_2];  // unscaled group-level effects
  vector<lower=0>[M_3] sd_3;  // group-level standard deviations
  vector[N_3] z_3[M_3];  // unscaled group-level effects
  vector<lower=0>[M_4] sd_4;  // group-level standard deviations
  vector[N_4] z_4[M_4];  // unscaled group-level effects
  vector<lower=0>[M_5] sd_5;  // group-level standard deviations
  vector[N_5] z_5[M_5];  // unscaled group-level effects
  vector<lower=0>[M_6] sd_6;  // group-level standard deviations
  vector[N_6] z_6[M_6];  // unscaled group-level effects
  vector<lower=0>[M_7] sd_7;  // group-level standard deviations
  vector[N_7] z_7[M_7];  // unscaled group-level effects
  vector<lower=0>[M_8] sd_8;  // group-level standard deviations
  vector[N_8] z_8[M_8];  // unscaled group-level effects
  vector<lower=0>[M_9] sd_9;  // group-level standard deviations
  vector[N_9] z_9[M_9];  // unscaled group-level effects
  vector<lower=0>[M_10] sd_10;  // group-level standard deviations
  vector[N_10] z_10[M_10];  // unscaled group-level effects
  vector<lower=0>[M_11] sd_11;  // group-level standard deviations
  vector[N_11] z_11[M_11];  // unscaled group-level effects
  vector<lower=0>[M_12] sd_12;  // group-level standard deviations
  vector[N_12] z_12[M_12];  // unscaled group-level effects
  vector<lower=0>[M_13] sd_13;  // group-level standard deviations
  vector[N_13] z_13[M_13];  // unscaled group-level effects
  vector<lower=0>[M_14] sd_14;  // group-level standard deviations
  vector[N_14] z_14[M_14];  // unscaled group-level effects
  vector<lower=0>[M_15] sd_15;  // group-level standard deviations
  vector[N_15] z_15[M_15];  // unscaled group-level effects
  vector<lower=0>[M_16] sd_16;  // group-level standard deviations
  vector[N_16] z_16[M_16];  // unscaled group-level effects
  vector<lower=0>[M_17] sd_17;  // group-level standard deviations
  vector[N_17] z_17[M_17];  // unscaled group-level effects
  vector<lower=0>[M_18] sd_18;  // group-level standard deviations
  vector[N_18] z_18[M_18];  // unscaled group-level effects
}
transformed parameters {
  // group-level effects
  vector[N_1] r_1_votetrump_1 = (sd_1[1] * (z_1[1]));
  // group-level effects
  vector[N_2] r_2_votetrump_1 = (sd_2[1] * (z_2[1]));
  // group-level effects
  vector[N_3] r_3_votetrump_1 = (sd_3[1] * (z_3[1]));
  // group-level effects
  vector[N_4] r_4_votetrump_1 = (sd_4[1] * (z_4[1]));
  // group-level effects
  vector[N_5] r_5_votetrump_1 = (sd_5[1] * (z_5[1]));
  // group-level effects
  vector[N_6] r_6_votetrump_1 = (sd_6[1] * (z_6[1]));
  // group-level effects
  vector[N_7] r_7_votetrump_1 = (sd_7[1] * (z_7[1]));
  // group-level effects
  vector[N_8] r_8_votetrump_1 = (sd_8[1] * (z_8[1]));
  // group-level effects
  vector[N_9] r_9_votetrump_1 = (sd_9[1] * (z_9[1]));
  // group-level effects
  vector[N_10] r_10_voteclinton_1 = (sd_10[1] * (z_10[1]));
  // group-level effects
  vector[N_11] r_11_voteclinton_1 = (sd_11[1] * (z_11[1]));
  // group-level effects
  vector[N_12] r_12_voteclinton_1 = (sd_12[1] * (z_12[1]));
  // group-level effects
  vector[N_13] r_13_voteclinton_1 = (sd_13[1] * (z_13[1]));
  // group-level effects
  vector[N_14] r_14_voteclinton_1 = (sd_14[1] * (z_14[1]));
  // group-level effects
  vector[N_15] r_15_voteclinton_1 = (sd_15[1] * (z_15[1]));
  // group-level effects
  vector[N_16] r_16_voteclinton_1 = (sd_16[1] * (z_16[1]));
  // group-level effects
  vector[N_17] r_17_voteclinton_1 = (sd_17[1] * (z_17[1]));
  // group-level effects
  vector[N_18] r_18_voteclinton_1 = (sd_18[1] * (z_18[1]));
}
model {
  vector[N_votetrump] mu_votetrump = temp_votetrump_Intercept + Xc_votetrump * b_votetrump;
  vector[N_voteclinton] mu_voteclinton = temp_voteclinton_Intercept + Xc_voteclinton * b_voteclinton;
  for (n in 1:N_votetrump) {
    mu_votetrump[n] += r_1_votetrump_1[J_1_votetrump[n]] * Z_1_votetrump_1[n] + r_2_votetrump_1[J_2_votetrump[n]] * Z_2_votetrump_1[n] + r_3_votetrump_1[J_3_votetrump[n]] * Z_3_votetrump_1[n] + r_4_votetrump_1[J_4_votetrump[n]] * Z_4_votetrump_1[n] + r_5_votetrump_1[J_5_votetrump[n]] * Z_5_votetrump_1[n] + r_6_votetrump_1[J_6_votetrump[n]] * Z_6_votetrump_1[n] + r_7_votetrump_1[J_7_votetrump[n]] * Z_7_votetrump_1[n] + r_8_votetrump_1[J_8_votetrump[n]] * Z_8_votetrump_1[n] + r_9_votetrump_1[J_9_votetrump[n]] * Z_9_votetrump_1[n];
  }
  for (n in 1:N_voteclinton) {
    mu_voteclinton[n] += r_10_voteclinton_1[J_10_voteclinton[n]] * Z_10_voteclinton_1[n] + r_11_voteclinton_1[J_11_voteclinton[n]] * Z_11_voteclinton_1[n] + r_12_voteclinton_1[J_12_voteclinton[n]] * Z_12_voteclinton_1[n] + r_13_voteclinton_1[J_13_voteclinton[n]] * Z_13_voteclinton_1[n] + r_14_voteclinton_1[J_14_voteclinton[n]] * Z_14_voteclinton_1[n] + r_15_voteclinton_1[J_15_voteclinton[n]] * Z_15_voteclinton_1[n] + r_16_voteclinton_1[J_16_voteclinton[n]] * Z_16_voteclinton_1[n] + r_17_voteclinton_1[J_17_voteclinton[n]] * Z_17_voteclinton_1[n] + r_18_voteclinton_1[J_18_voteclinton[n]] * Z_18_voteclinton_1[n];
  }
  // priors including all constants
  target += normal_lpdf(temp_votetrump_Intercept | 0, 1);
  target += normal_lpdf(temp_voteclinton_Intercept | 0, 1);
  target += student_t_lpdf(sd_1 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_1[1] | 0, 1);
  target += student_t_lpdf(sd_2 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_2[1] | 0, 1);
  target += student_t_lpdf(sd_3 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_3[1] | 0, 1);
  target += student_t_lpdf(sd_4 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_4[1] | 0, 1);
  target += student_t_lpdf(sd_5 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_5[1] | 0, 1);
  target += student_t_lpdf(sd_6 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_6[1] | 0, 1);
  target += student_t_lpdf(sd_7 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_7[1] | 0, 1);
  target += student_t_lpdf(sd_8 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_8[1] | 0, 1);
  target += student_t_lpdf(sd_9 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_9[1] | 0, 1);
  target += student_t_lpdf(sd_10 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_10[1] | 0, 1);
  target += student_t_lpdf(sd_11 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_11[1] | 0, 1);
  target += student_t_lpdf(sd_12 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_12[1] | 0, 1);
  target += student_t_lpdf(sd_13 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_13[1] | 0, 1);
  target += student_t_lpdf(sd_14 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_14[1] | 0, 1);
  target += student_t_lpdf(sd_15 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_15[1] | 0, 1);
  target += student_t_lpdf(sd_16 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_16[1] | 0, 1);
  target += student_t_lpdf(sd_17 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_17[1] | 0, 1);
  target += student_t_lpdf(sd_18 | 3, 0, 4)
    - 1 * student_t_lccdf(0 | 3, 0, 4);
  target += normal_lpdf(z_18[1] | 0, 1);
  // likelihood including all constants
  if (!prior_only) {
    target += bernoulli_logit_lpmf(Y_votetrump | mu_votetrump);
    target += bernoulli_logit_lpmf(Y_voteclinton | mu_voteclinton);
  }
}
generated quantities {
  // actual population-level intercept
  real b_votetrump_Intercept = temp_votetrump_Intercept - dot_product(means_X_votetrump, b_votetrump);
  // actual population-level intercept
  real b_voteclinton_Intercept = temp_voteclinton_Intercept - dot_product(means_X_voteclinton, b_voteclinton);
}