In the Data Analysis Using Regression and Multilevel/Hierarchical Models (ARM) book by Gelman & Hill, it discusses Five ways to write the same model in section 12.5.

The final example is a Large regression with correlated errors. It defines the covariance matrix as:

For any unit i: var_y+var_alpha
For any units i,k within the same group j: var_alpha
For any units i,k in different groups: cov(error_i, error_k)=0

What if instead, these are adjusted to:

For any unit i: var_y+var_alpha+var_base
For any units i,k within the same group j: var_alpha+var_base
For any units i,k in different groups: cov(error_i, error_k)=var_base

where var_base is some base level of correlation among the variables.

I would like to do this in Stan without having to write the explicit covariance matrix. I tried some naive attempts to adopt a model similar to the radon_no_pool.stan example model to account for it, but didn’t have much luck (divergences, tree depth issues, Bayesian fraction of information issues…the works). I suspect it’s because I was basically introducing a latent variable that couldn’t be identified properly.

So I suppose the big question is: am I wasting my time thinking about this?

I’m fitting a cross-sectional model for simplicity, but the full dataset is a panel. I was just curious if I can incorporate the correlation structure that is evident when looking at the full panel in the cross-section.

I don’t think I am following you here, could you elaborate? In particular I don’t understand how a multilevel model is related to covariance matrices (I didn’t read the ARM book, and am not an expert so I might be missing something basic). Maybe the thing you are trying is simply adding one more level to the model?

In any case, it would be useful if you posted your model code - my gut feeling is that the model should be identifiable, so it might be that the problem is somewhere else in the model, or your data actually do not fit the model very well.

The point the book makes is that the typical hierarchical model is equivalent to a large regression with correlated errors. In my example, I was looking at a cross-sectional slice of a panel. I know that the panel had some correlation between the variables on a time-series basis and I figured I might test to see if it also showed up if I just did a cross-sectional slice of the data.

One way to write the original model was:
y_i ~ N(X_i*B + n_j[i], sigma_y)
n_j ~ N(0, sigma_a)

So I think what I’m talking about is actually
y_i ~ N(X_i*B + n_j[i] + v, std_y)
n_j ~ N(0, std_alpha)
v ~ N(0, std_base)

which is why it was hard to identify.

Anyway, after your post, I was thinking about this a little more today… The issue I was having trouble wrapping my head around was that I wanted to start with a cross-sectional slice and move to a panel. I knew correlation was an issue in the panel, so I figured that I should think about it in the cross-section as well. However, in the cross-section, the intercept should be capturing the average level at that period, which is one of the major drivers of the effect within the panel. So I think I’ve convinced myself that it’s not something to worry about.