Can Generalized Least Squares be used to fit a hierarchical model?

Question

I'm used to approaching the idea of a hierarchical model from the Bayesian perspective, such as the one developed in Andrew Gelman's book on regression. I find that I can understand these models from the generative model perspective, but I'm having trouble porting that understanding to the classical Generalized Least Squares framework. My question is about whether or not these hierarchical models can be implemented with GLS.

The simplest example I can think of is the following intercept-only hierarchical model. Imagine that we have $J$ groups. We measure a number of units for each group, taking note of some measured quantity, $X$, for each unit.

We then suppose the following model for our observations:

$X_i \sim N(\alpha_{[j]i}, \sigma)$

$\alpha_j \sim N(\mu_a, \sigma_a)$

We can assume some hyperprior over $\mu_a, \sigma_a$. The formulation above with the hyperprior allows us to write down the full posterior over the parameters. A Bayesian approach would be to compute the posteriors over all the parameters, conditioned on the data, using a tool like MCMC or Variational Inference. This sort of analysis would take in my observed $X_i$ and spit out posterior samples of the parameters, from which I could make my inference. This is an attractive method because if a unit is part of two groups, it is not hard to add a second set of cluster variables, priors and hyperpriors to the above model. As far as I understand, no OLS model specification will allow us to estimate the observed model, nor will and WLS model; there is something about the between-unit relationships in the above model which cannot be represented in the OLS or WLS framework.

I'm trying to understand whether this kind of model can be implemented at all in the GLS framework, or if there's some other tool I should be looking for. Of course, there will be no hyperpriors, which is fine. But it seems plausible to me that the above hierarchical model can be specified somehow by the correlation structure between the errors that we feed into a GLS model with our data. I know that GLS and random effects models (which seem to be very similar to the Bayesian Hierarchical methods) are closely related, but I'm having trouble understanding how. Is there some way to specify the covariance matrix, such as the parameter in statsmodels, that lets us estimate the above model? How about Hierarchical models more generally?

Answer 1

But it seems plausible to me that the above hierarchical model can be specified somehow by the correlation structure between the errors that we feed into a GLS model with our data.

That's exactly how it's done. The residual variance $\sigma^2$ and the higher level variance $\sigma^2_a$ are uknown parameters, however the covariance matrix of the errors of the model is completely determined by those two parameters:

$$ \text{Cov}[y_i, y_j]= \cases{ \sigma^2 + \sigma^2_a & if i = j\\ \sigma^2_a & if i and j belong to the same cluster\\ 0 & otherwise } $$

Then, GLS does the trick.

When the model is more complex, and takes also random slopes into account, the covariance matrix gets quite more complicate, but still, it can be computed starting from the matrix of random covariates.

I tried to find more thorough instructions in my university book, but I couldn't find them, so the best document that comes to my mind is this quite off topic paper that shows how to compute the covariance matrix of residuals in a completely different context, but that is still accurate.

Answer 2

I don't have enough reputation to create a comment but I simply want to suggest you look into Rob Hyndman's 2011 paper in which he proposes an approach to creating optimal hierarchical forecasts using generalized least squares. I imagine it is possible to apply a similar technique to non-time-series data. His 2017 paper mentions the shortcomings of this approach and discusses other approaches.