factor-analytic/reduced-rank models

[Ben Bolker]

Just floating this in case there are others who have tried this
before ...

glmmTMB has an option to fit models with
reduced-rank/factor-analytic covariance structures: that is, we assume
the latent variables for each group/block are drawn from a multivariate
normal distribution with a rank-deficient covariance matrix, Sigma =
Lambda Lambda', where Lambda is n*d (n is the number of variables per
block, d is the rank of Sigma).

The way we've done this so far is to use the *spherical* (iid Normal)
latent variables, u, with u ~ MVN(0, sigma^2 I) and then computing b =
Lambda*u to get the predicted values for the conditional log-likelihood.
This works fine, but things blow up when we try to do some computations
to calculate the hat matrix/leverage on large data sets (ridiculously
large memory requirements).

My guess is that this happens because of the pattern of the Hessian
matrix of the latent variables. If we work with spherical latent
variables u, each latent variable u_i affects all of the non-spherical
latent variables (b) for that block; if we work with non-spherical
variables, there should be less overlap in effect (maybe).

There are some technical problems with computing a multivariate
normal distribution with a non-full-rank covariance matrix
(<https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case>);
I thought I could deal with these problems, and some others, by working
in the *precision* space instead, i.e. parameterizing the model with the
elements of the (pseudo)*inverse* factor loading matrix Phi =
pseudo-inv(Lambda) (i.e. Phi*Phi' = pseudo-inv(Sigma).

I have banging my head on this off and on for a few days; I will
presumably work it out eventually (or give up), but was curious if
anyone had tried anything like this before/had pointers or references. A
quick look at projects like Jim Thorson's spatial factor-analytic models
<https://github.com/James-Thorson/spatial_factor_analysis> suggests that
they do things the straightforward way (i.e. using Laplace approx on
spherical latent variables and transforming to b = Lambda*u for
computing predictions ...

cheers
Ben Bolker

[james...@gmail.com]

Two quick thoughts:

1. I agree that the space-space version (i.e., with projected factors as random effects) is often more computationally efficient, but also doesn't apply in the rank-deficient case, so in dsem/tinyVAST etc I allow the user to instead select the MVN on unprojected factors and then projecting them afterwards, e.g., using `dsem_control(  "gmrf_parameterization" )` ;
2. Zheng et al. In review here (https://arxiv.org/abs/2411.14185) has demo code for computing the diagonal of the hat matrix for GLMMs with continuous responses which it then uses to compute effective degrees of freedom (I think it's Eq-4 and associated code) ... am I right in guessing that you're only interested in the diagonal of the hat matrix, or the full matrix?  I haven't actually compared the calculation with conventional leverage metrics though.  And Nan would be the expert to address follow up questions on that

Jim