Mixed distribution columns

[spinkney]

In the Collins paper they write

Would it be possible to add the ability for each column in Y to be distributed by a different exponential family distribution?

[willtownes]

Yes, absolutely that is possible in principle, however it's not currently implemented. One difference between this implementation and that of Collins is we do not require the factorized parameters to be the natural parameter- rather it should just be the mean. This is particularly relevant for the negative binomial likelihood since the log link is not canonical. Can you describe what kind of data would have this characteristic?

[willtownes]

By the way here is a similar idea as what you are suggesting applied to unsupervised clustering: https://arxiv.org/abs/1510.05491

[spinkney]

I've seen your package when it was released and I think this is a great idea. I'm a Stan developer and would love to write a probabilistic glmpca (just like bayesian gaussian pca has been done). I mentioned your package on the Stan forums as a method to use. The person responded that they need mixed type PCA. I went sleuthing through the Collins paper to see more about the method and was surprised to see that paragraph.

I don't have a specific use case outside of that post on the Stan forums but I can see this being useful for my work in the future. I just completed a bunch of bayesian copula stuff and just helped with a mixed continuous-discrete margin case. Having a nice latent space representation of mixed outcomes would be a really nice tool to add.

[rmcminds]

I am that user on the Stan forums, and without having dug too deeply into your package yet, my use case is a multi-omics study that has 4 HTS count datasets, a number of continuous datasets such as temperature measurements and telomere lengths, and a number of categorical/binary datasets with observations such as 'is the organism in contact with this other species in a photo'. Obviously the categorical stuff in particular, and also much of the continuous data, could be used as 'independent variables' in a glm analyzing the count data, but I think there are theoretical benefits to doing something very unsupervised with even those, given that they are not directly manipulated experimental treatments and may contain error themselves. There are also quite a few of them, without clear precedence for selecting just a few and testing them.

The linked preprint looks interesting but a bit dense for me on a first pass!

[willtownes]

Sounds like y'all are doing some really interesting work! Unfortunately the bottom line is I will probably not be able to add this feature anytime soon. However, I'll keep the issue open for if/when I can find a scrap of free time to try it out.

One package I would also recommend checking out is Jeff Miller's generalized bilinear models in Julia. It is more theoretically sophisticated than glmpca and I think it supports interactions better, but I can't remember if it can handle mixed likelihood outcomes.

I glanced at the Stan forum discussion and one suggestion (and what is implemented in glmpca) is to avoid trying to implement constraints within Stan, instead, just do a post-processing step similar to what I described in section 5 of the glmpca technical report. Unfortunately this will only work with point estimation (variational inference or maximum a posteriori). If you are doing MCMC, the latent factors will probably be unidentifiable and you can't get the posterior mean by averaging across the samples as usual. Instead, it is usually best to compute an identifiable quantity (such as the outer product of the loadings matrix VV') that is invariant to rotations then average that across the samples. In my experience uncertainty quantification is less essential for these dimension reduction models anyway, at least in comparison to regression models, and the speed and scalability concerns are more important.

This is something you probably already know how to do, but here's a primitive implementation of Bayesian factor analysis in Stan that could potentially be generalized to mixed likelihood outcomes.

[rmcminds]

Thanks for the suggestions and for keeping the issue open!

In my experience uncertainty quantification is less essential for these dimension reduction models anyway, at least in comparison to regression models, and the speed and scalability concerns are more important.

Yes, you are absolutely right! I am stubbornly still thinking it'd be nice to have my cake and eat it, too, but for all practical purposes we should focus on the point estimates. I'm happy to see that your package does something similar to my suggested post-processing, and that may be all I need to justify it, cite it, and move on with my analysis for the moment. Thanks!

[spinkney]

Thanks for that! I checked out the generalized bilinear model package but it only has negative binomial. Though this is a really cool package. It reminds me of doing a regression using matrix variate normals but now the generalized matrix variate model.

Also that stan code you have can be made more performant by declaring the parameters as matrices. Since they're all getting standard normal iid priors you can vectorize the call. The std_normal() is faster than normal(). Lastly you can vectorize the observation call by using vectorized indices. The rows_dot_product() does the dot product row wise. It's about 33% faster for me. But that's only a 1s difference for this small code.

parameters {
  matrix[N, L] u;
  matrix[G, L] v;
}
model {
  // prior
  to_vector(u) ~ std_normal();
  to_vector(v) ~ std_normal();
  
  // likelihood
  y ~ normal(rows_dot_product(u[nn], v[gg]), 1);
}