Log-linear models applied to contingency tables where the cells are modeled as independent Poisson variables has a brief mention in Ch16 p428-431 BDA3, but I haven’t been able to find a Stan model anywhere illustrating the core ideas: especially featuring hierarchical models with interactions and main effects. Does anyone know of any examples?
@martinmodrak, since it’s widely used, if there isn’t any examples, would the Stan team be interested in one for the documentation section?
I am unfortunately completely unfamiliar with these types of model as well as with their applications (but I guess I understand what you mean). So it is hard for me to judge the importance/wide use.
I think a general rule about examples is that we love ourselves a good example :-). Now, aiming directly for the Stan docs might not always be the best course of action, especially since the devs may take a while to react (although AFAIK @Bob_Carpenter tends to respond quickly there and knows way more details about the process of inclusion in the docs than I do). I would personally recommend people to just write an example (possibly in a similar style as the Stan docs are written), put it somewhere on the Interwebs and advertise here on Discourse (or even post the whole example directly as a Discourse topic). This way people can use it immediately. Inclusion in Stan docs can then happen on its own pace.
Does that make sense?
Totally, thanks for the advise. I’ll do accordingly.
@emiruz have you made any progress on this? I’m starting a problem that is essentially a hierarchical contingency table problem (think of it as the 8-schools problem, but where we get a 2x2 (or larger) contingency table from each school, and I think a Bayesian log-linear model would be useful. I haven’t seen any good documentation on appropriate priors for this type of problem, although the actual model itself in Stan shouldn’t be hard to implement.
Not yet; I’ve drowning in other things at the moment, but I intend to. It’s a really interesting class of models.
Dr A Overstall wrote the conting package (https://cran.r-project.org/web/packages/conting/index.html). Check out the vignettes because it may do what you need out of the box.
The really interesting part of Overstall’s work comes when you assume that a contingency table has interaction terms but you don’t know what they are. He and others worked out an elegant prior on the space of possible designs and an efficient Gibbs sampler for it. That part would be difficult to do in Stan since it’d have to be discrete.
I was wondering whether you have made any progress on this? im looking to implement some log-linear models for contingency tables in Stan
As long as you can follow our style, we love examples in the User’s Guide. The way to do that is create a pull request in the stan-dev/docs repo.
We’d just marginalize out discrete parameters, which many people find difficult to code. There are examples in the user’s guide chapter on latent discrete parameters. It’s usually way
more efficient to marginalize out discrete parameters than to sample them with Gibbs unless you can be very clever with blockwise Gibbs sampler (i.e., roll your own rather than using BUGS/JAGS).
These models are simple enough that you can probably just code them as written mathematically in Stan. If someone can be more specific about which model they want, I’d be happy to code it into a simple Stan program. For example, here’s a 2-way “independence” Poisson model as described in BDA3, section 16.7 and also in Alan Agresti’s An Introduction to Categorical Data Analysis in section 7.1 (log linear models for two-way and three-way tables):
data {
int<lower=0> M;
int<lower=0> N;
array[M, N] int<lower=0> y;
}
parameters {
real alpha;
vector[M] beta;
vector[N] gamma;
}
model {
alpha ~ ... prior ...
beta ~ ... prior ...
gamma ~ ... prior ...
for (m in 1:M) {
for (n in 1:N) {
y[m, n] ~ poisson_log(alpha + beta[m] + gamma[n]);
}
}
}
Here, poisson_log(lambda) = poisson(exp(lambda)), so it’s still the log linear model being suggested.
I don’t know much about these models, but the prior needs to be sensitive to the scale of the counts or there needs to be something like a constant rate term plus logit link.
Then if you want to build the “saturated” model, you add another predictor
matrix[M, N] delta;
and give it a prior and add it to the linear model
y[m, n] ~ poisson_log(alpha + beta[m] + gamma[n] + delta[m, n]);
These models can be hard to identify with weakly informative independent priors because of the additive non-identifiabilities. You can try to deal with the non-identifiabilities as @andrewgelman et al. suggest and pin some of the values in beta, gamma, and delta to zero for identifiability. Or set them up so that they sum to zero.
P.S. I didn’t understand @andrewgelman et al.'s discussion in the Prior Distributions section because it suggests putting a prior on the derived quantity mu[m, n] = alpha[m] + beta[n]. We can do that here directly because the Jacobian determinant is constant, but it’s going to be overparameterized in that we get m * n prior terms for m + n parameters.
[edit: escaped operators]
It’s ok to have m*n prior terms for m+n parameters. I mean, it could be ok, or it could be terrible, depending on the details. I’m speaking here in general terms without thinking too hard about this particular example. The trick to thinking about this sort of prior information is to think of each term in the prior as a separate piece of information. So if you can think of these as m*n pieces of info, you’re cool; if you’re using them to represent a population distribution, then, yeah, that won’t work because the parameters are entangled in the model.
How about the model from this paper https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.5695
specifically the model at the top of page 2575 for 3 tests (P=3)
The paper’s paywalled and I don’t have access. Can you summarize your question here?