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