Information criteria for parameterizations as multivariate Poisson

Modeling
1

This post is especially for @avehtari and others who are well-versed in questions about information criteria.

I have a Poisson regression model that, in principle, has a single response variable. However, because of a few somewhat complicated random effects structures (not show in this example), it made sense to split up the estimation into multiple groups/functions. For simplicity, imagine we are estimating the respective outcomes for males and females via separate likelihoods. Assuming different samples sizes for the two (N_A1 and N_A2), the model then looks something like this:

model {
	//group 1 -- males
	for ( i in 1:N_A1 ) {
	lambda_A1[i] = a
	;
	}
	   
	//group 2 -- females
	for ( i in 1:N_A2 ) {
	lambda_A2[i] = a
	;
	}
	
	y_A1 ~ poisson_log( lambda_A1 );
	y_A2 ~ poisson_log( lambda_A2 );
}

Meanwhile, my colleagues are interested in having information criteria (LOOIC, WAIC, etc.) for the models that we estimate. In general, I understand that the log_likelihood is not well defined for multivariate Poisson models. However, question number 1 – my assumption is that since this isn’t truly a multivariate response, it would be okay to estimate the log_lik for each of the functions. Is that right?

Assuming that it would be okay, I’m not quite sure how to generate log_lik quantities that are compatible with some of the convenience functions that take the log_lik from Stan models to generate summary calculations. In other words, the following represents my first pass at doing this, with the evident need to generate a separate log_lik quantity for each of the functions:

generated quantities{
    vector[N_A1] log_lik_A1;
    vector[N_A2] log_lik_A2;
    real lambda_A1;
    real lambda_A2;
    
    lambda_A1 = a;
    lambda_A1 = exp(lambda_A1);
    
    lambda_A2 = a;
    lambda_A2 = exp(lambda_A2);
    
    for ( i in 1:N_A1 ) log_lik_A1[i] = poisson_lpmf( y_A1[i] | lambda_A1 );
    
    for ( i in 1:N_A2 ) log_lik_A2[i] = poisson_lpmf( y_A2[i] | lambda_A2 );    
}

However, it seems that the log_lik for the respective observations would ideally end up as a single vector of log_liks rather than the two separate vectors that are defined here. Question number 2 – is there a way to accomplish that?

(The alternative would be for me to recode the model as a single likelihood function, but the indexing will get somewhat intense because there are actually closer to 30 sets of correlated random effects, not 2, but with common variances/covariances and fixed effect parameters across the functions.)

1 Like
2

In case of group or hierarchical structure, I recommend to think in terms of cross-validation as it makes more clear the prediction task or focus on certain model level (while information criteria language is hiding this)

Why it wouldn’t be defined?

It’s unclear what is function here. Do you want to predict or focus on new groups, individuals or observations? (It’s not clear what your observation process is)

Assuming your focus is on level of i, you can just concatenate log_lik_A1 and log_lik_A2 (see 5.10 Matrix Concatenation | Stan Functions Reference)

Use poisson_log_lpmf instead of poisson_lpmf to avoid the exponentation

2 Likes
3

Thanks, Aki. It sounds like it’s conceptually straightforward, and a fairly straightforward append_row[] in the generated quantities would create the vector of log_lik that we need. Something like the following should work, correct?

generated quantities{
    vector[N_A1] log_lik_A1;
    vector[N_A2] log_lik_A2;
    real lambda_A1;
    real lambda_A2;
    vector[N] log_lik;

    lambda_A1 = a;
    lambda_A1 = exp(lambda_A1);
    
    lambda_A2 = a;
    lambda_A2 = exp(lambda_A2);
    
    for ( i in 1:N_A1 ) log_lik_A1[i] = poisson_lpmf( y_A1[i] | lambda_A1 );
    
    for ( i in 1:N_A2 ) log_lik_A2[i] = poisson_lpmf( y_A2[i] | lambda_A2 );

    log_lik = append_row [log_lik_A1, log_lik_A2];
}
4

Yes (without checking whether it actually compiles)

1 Like