WAIC for probit model confirmation or troubleshooting?

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) );
	
	}
}
"

Code seems fine, but computing WAIC is not. Use LOOIC from loo() instead.

1 Like

Thanks for the tip, Ben. I went ahead and gave this a try and for other users wishing to do likewise, I found the help documentation to be helpful: help(loo)

It shows how to extract the log likelihoods from the stan model.

log_likb <- extract_log_lik(mb)
loo(log_likb)

My discipline is just getting familiar with WAIC, so LOOIC will be something new for most. Is there a canonical reference that details the advantages of LOOIC?

Same as the one in the documentation for the loo package:
http://www.stat.columbia.edu/~gelman/research/published/loo_stan.pdf

1 Like

Andrew’s preprint mentions “To appear”, but in case you want to cite it, it has appeared in Statistics and Computing, 27(5):1413–1432.

Aki