Hi everyone,

I’m analyzing some dyadic data with a probit model and thinking about calculating WAIC for model selection using built-in functions in *loo*. I think that the following code in “generated quantities” will produce the **log_lik** values that feed into *loo*, but it’d be helpful to get some confirmation of that because my first attempt returned a moderately higher value than expected (compared against a simpler model). I’m debugging other aspects of the code now, but thinking that it’d be good to confirm the **log_lik** portion of the code as well. Suggestions are appreciated. Thanks, Jeremy

```
model_code_b <- "
data{
int N;
int N_i__j_ID;
int N_ij_dyads;
int y[N];
int i_ID[N];
int j_ID[N];
int ij_dyad_ID[N];
}
parameters{
real a;
// Node level effects cov matrix
matrix[2,N_i__j_ID] z_i__j_ID;
vector<lower=0>[2] sigma_i__j_ID;
cholesky_factor_corr[2] L_Rho_i__j_ID;
vector[N_ij_dyads] v_ij_dyad_ID;
real<lower=0> sigma_ij_dyad_ID;
}
transformed parameters{
// Define all matrices
matrix [N_i__j_ID,2] v_i__j_ID;
{
v_i__j_ID = (diag_pre_multiply ( sigma_i__j_ID, L_Rho_i__j_ID ) * z_i__j_ID)';
}
}
model{
vector[N] p;
//priors
a ~ normal( 0 , 5 );
to_vector (z_i__j_ID) ~ normal (0,1);
v_ij_dyad_ID ~ normal( 0 , 1 );
sigma_ij_dyad_ID ~ exponential (1);
sigma_i__j_ID ~ exponential (1);
L_Rho_i__j_ID ~ lkj_corr_cholesky(4);
// LIKELIHOOD
for ( i in 1:N ) {
p[i] =
a
+ v_i__j_ID [i_ID[i],1]
+ v_i__j_ID [j_ID[i],2]
+ v_ij_dyad_ID [ ij_dyad_ID[i]] * sigma_ij_dyad_ID
;
}
// update target
y ~ bernoulli (Phi(p));
}
generated quantities{
matrix[2,2] Rho_i__j_ID;
vector[N] log_lik;
Rho_i__j_ID = L_Rho_i__j_ID * L_Rho_i__j_ID';
for ( i in 1:N ) {
real p;
p =
a
+ v_i__j_ID [i_ID[i],1]
+ v_i__j_ID [j_ID[i],2]
+ v_ij_dyad_ID [ ij_dyad_ID[i]] * sigma_ij_dyad_ID
;
log_lik[i] = bernoulli_lpmf ( y[i] | Phi(p) );
}
}
"
```