election88.stan

data {
  int<lower=0> N;
  int<lower=0> n_age;
  int<lower=0> n_age_edu;
  int<lower=0> n_edu;
  int<lower=0> n_region;
  int<lower=0> n_state;
  array[N] int<lower=1, upper=n_age> age;
  array[N] int<lower=1, upper=n_age_edu> age_edu;
  vector<lower=0, upper=1>[N] black;
  array[N] int<lower=1, upper=n_edu> edu;
  vector<lower=0, upper=1>[N] female;
  array[n_state] int<lower=1, upper=n_region> region;
  array[N] int<lower=1, upper=n_state> state;
  vector[n_state] v_prev;
  array[N] int<lower=0, upper=1> y;
}
parameters {
  vector[n_age] b_age;
  vector[n_age_edu] b_age_edu;
  vector[n_edu] b_edu;
  vector[n_region] b_region;
  vector[n_state] b_state;
  real b_v_prev;
  vector[4] beta;
  real mu;
  real mu_age;
  real mu_age_edu;
  real mu_edu;
  real mu_region;
  real<lower=0, upper=100> sigma_age;
  real<lower=0, upper=100> sigma_edu;
  real<lower=0, upper=100> sigma_age_edu;
  real<lower=0, upper=100> sigma_region;
  real<lower=0, upper=100> sigma_state;
}
transformed parameters {
  vector[n_age] b_age_adj;
  vector[n_age_edu] b_age_edu_adj;
  vector[n_edu] b_edu_adj;
  vector[n_region] b_region_adj;
  vector[n_state] b_state_hat;
  real mu_adj;
  vector[N] Xbeta;
  vector[N] p;
  vector[N] p_bound;
  
  for (i in 1 : N) {
    Xbeta[i] = beta[1] + beta[2] * female[i] + beta[3] * black[i]
               + beta[4] * female[i] * black[i] + b_age[age[i]]
               + b_edu[edu[i]] + b_age_edu[age_edu[i]] + b_state[state[i]];
  }
  
  mu_adj = beta[1] + mean(b_age) + mean(b_edu) + mean(b_age_edu)
           + mean(b_state);
  b_age_adj = b_age - mean(b_age);
  b_edu_adj = b_edu - mean(b_edu);
  b_age_edu_adj = b_age_edu - mean(b_age_edu);
  b_region_adj = b_region - mean(b_region);
  
  for (j in 1 : n_state) {
    b_state_hat[j] = b_region[region[j]] + 100 * b_v_prev * v_prev[j];
  }
}
model {
  mu_age ~ normal(0, 1);
  mu_edu ~ normal(0, 1);
  mu_age_edu ~ normal(0, 1);
  mu_region ~ normal(0, 1);
  mu ~ normal(0, 100);
  sigma_age ~ uniform(0, 100);
  sigma_edu ~ uniform(0, 100);
  sigma_age_edu ~ uniform(0, 100);
  sigma_region ~ uniform(0, 100);
  sigma_state ~ uniform(0, 100);
  
  beta ~ normal(0, 100);
  b_age ~ normal(100 * mu_age, sigma_age);
  b_edu ~ normal(100 * mu_edu, sigma_edu);
  b_age_edu ~ normal(100 * mu_age_edu, sigma_age_edu);
  b_state ~ normal(b_state_hat, sigma_state);
  
  b_v_prev ~ normal(0, 1);
  b_region ~ normal(100 * mu_region, sigma_region);
  
  y ~ bernoulli_logit(Xbeta);
}