Computing pointwise log-likelihood for a binomial normal hierarchical model

Dear Stan experts,

I want to compute pointwise log-likelihood for a binomial-normal hierarchical model. Then, I want to do leave-one-out cross validation. Precisely, my model is "Blocker: “random effects meta-analysis of clinical trials” (https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/blocker/blocker.stan)

I followed the explanations given by “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC” (Vehtari et al.). Although there is Stan code in that paper which is actually very similar to my model, still I don’t get how to do it. Because, there is two likelihood functions in my model, one for control arm and the other is for the treatment arm:

...
model{
...
// likelihood
rt ~ binomial_logit(nt, mu + delta);
rc ~ binomial_logit(nc, mu);
}

Firstly, I thought I should write as follows:

generated quantities {
  vector[N] log_lik_c;
  vector[N] log_lik_t;
  for (n in 1:N)
    log_lik_c[n] = bernoulli_logit_lpmf(rc[i] | nc[i], mu[i]);
    log_lik_t[n] = bernoulli_logit_lpmf(rt[i] | nt[i], mu[i] + delta[i]);
}

But, now there are two “log_lik” vectors? How should I combine them? Because, afterwards I want to use “loo” package and compute “elpd” etc.

Then, I thought maybe I should create a new variable for treatment (trt = 0 for control and trt = 1 for treatment), and then compute pointwise log-likelihood as follows:

generated quantities {
vector[2*N] log_lik;

for (i in 1:(2*N)) {
  if (trt[i] == 0) {
    log_lik[i] = binomial_logit_lpmf(rc[i] | nc[i], mu[i]);
}  else {
     log_lik[i] = binomial_logit_lpmf(rt[i] | nt[i], mu[i] + delta[i]);
} 

}
}

Is one of the two approaches correct, or am I missing something?

Thank you very much for your help,

Best regards,

Burak Kürsad Günhan

Your second approach is on the right track. You end up with one one vector for LOO/etc…

Yeah, if you want to compute LOO the way we describe in our paper you reference then you’ll want to have one vector of length N+N (assuming both rt and rc are length N in this case, although they need not be equal in general). Then when you run Stan you’ll end up with a S by 2N matrix of draws of the point-wise log-likelihood. You can then use that matrix with the functions in the loo R package.

After creating log_lik_c and log_lik_t, you could also have
log_lik = append_row(log_lik_c, log_lik_t);

Aki

just a follow up to this discussion. Let’s say you have multiple terms in the likelihood (as is presented in the above example), but you want the marginal contribution of one of them for calculating LOO. It doesn’t make sense in the example above, but I’m thinking of a case where one sampling statement is defined for a latent variable (which is typically integrated out of the likelihood for predictive assessment). Would it be as simple as just ignoring the contribution of the latent variable term when generating the log_lik vector in generated quantities?

Yes.

Aki

Thank you very much for all three answers. Both suggested approaches works. But of course Aki’s approach is simpler.

Burak Kürsad Günhan