Stand-alone predictions in the context of a hierarchical random effect model

Interfaces RStan
,
2

Estimating the group-level parameters for new observations is a bit more of a complicated endeavour, and was discussed more over in this thread: Posterior Predictive for hierarchical model without further hyper-parameter inference - #8 by maxbiostat

Note that you can pretty significantly optimise your model by using the non-centered parameterisation for u_i and using the compound bernoulli_logit_glm distribution for your likelihood, like so:

data {
  int<lower=1> N;                                     // total number of observations
  int<lower=0,upper=1> outcome[N];                    // outcome
  int<lower=1> N_people;                              // number of unique individuals
  int<lower=1,upper=N_people> seq_person_id[N];       // SEQUENTIAL person id -- important for indexing
  matrix[N, 1]    predictor;                             // predictor that goes with random slope
}
parameters {
  // variance components
  vector<lower = 0>[2]  tau_u;
  // 'Raw' N(0,1) variates for the non-centered param
  matrix[2, N_people] u_i_raw;
  // correlation matrix
  cholesky_factor_corr[2] rho_matrix_cholesky;

}
transformed parameters {
  // random effects for people
  // Implies u_i ~ multi_normal( [0, 0]', Sigma )
  // Specified with column per group (person) since it's
  // faster to extract a column than a row
  matrix[2, N_people] u_i = diag_pre_multiply(tau_u, rho_matrix_cholesky) * u_i_raw;
}
model {
  // correlation
  rho_matrix_cholesky ~ lkj_corr_cholesky(1);
  // variances
  tau_u ~ cauchy(0,1);

  // 'Raw' N(0,1) variates for the non-centered param
  to_vector(u_i_raw) ~ std_normal();

  // outcome
  outcome ~ bernoulli_logit_glm(predictor,
                                u_i[1, seq_person_id]',
                                u_i[2, seq_person_id]');
}
1 Like