Background
Gelman and Hill's Data Analysis Using Regression and Multilevel/Hierarchical Models includes an example in section 13.5 of how to model non-nested data. The second example in this section is a regression of log earnings, $y_i$, on height, $z_i$ (mean-adjusted), and ethnicity and age categories. Using the notation of the book, there are $J = 4$ ethnic groups and $K=3$ age categories.
The multilevel model there is written
$$y_i \sim N(\alpha_{j[i],k[i]} + \beta_{j[i],k[i]}z_i, \sigma_y^2), \text{ for all $i=1,\dots,n$} ,$$ where
$$ \begin{pmatrix} \alpha_{j,k} \\ \beta_{j,k} \end{pmatrix} = \begin{pmatrix} u_0 \\ u_1 \end{pmatrix} + \begin{pmatrix} \gamma_{0j}^{eth} \\ \gamma_{1j}^{eth} \end{pmatrix} + \begin{pmatrix} \gamma_{0j}^{age} \\ \gamma_{1j}^{age} \end{pmatrix} + \begin{pmatrix} \gamma_{0j}^{eth\times age} \\ \gamma_{1j}^{eth\times age} \end{pmatrix}, $$ with priors
$\begin{pmatrix} \gamma_{0j}^{eth} \\ \gamma_{1j}^{eth} \end{pmatrix}\sim N(\mathbf0, \mathbf{\Sigma}^{eth})$
$\begin{pmatrix} \gamma_{0j}^{age} \\ \gamma_{1j}^{age} \end{pmatrix}\sim N(\mathbf0, \mathbf{\Sigma}^{age})$
$\begin{pmatrix} \gamma_{0j}^{eth\times age} \\ \gamma_{1j}^{eth\times age} \end{pmatrix}\sim N(\mathbf0, \mathbf{\Sigma}^{eth\times age})$.
Using lmer() this is modeled as
lmer (y ~ x.centered + (1 + x.centered | eth) + (1 + x.centered | age) +
(1 + x.centered | eth:age))
Question
The above brilliantly, to my mind, allows varying slopes and intercepts to be modeled in a non-nested way. Neither age nor ethnicity is "contained" with one another, yet one can still leverage the partial-pooling power of multilevel models.
In the setting I care about (with categories having nothing to do with race and age) $J=100$ and $K=20$. But in my setting the variable corresponding to $J$ has a nested structure (with 3 or 4 levels). Suppose that $J$ is the number of cities, which are clustered in counties, states, and countries (assuming for simplicity each country has the same hierarchical categories as the U.S.).
My question is, can one usefully specify a multilevel-model with this "partially-nested, partially non-nested" structure? This need not be able to be fit using lmer().
Vague idea #1: use the same specification as in the non-nested background example, but impose some sort of prior on the covariance of the 4-level-nested-category which makes the otherwise-noisy covariance matrix more robustly specified.
Vague idea #2: Compute the group level averages for each level of the hierarchy and include them as fixed effects.
I would prefer a specification that allows country/state/county and city effects to all be random, but it each of the lower levels would seem to require a centering on the level above it, which seems to complicate things. I also don't have a uniform amount of data in each $j\times k$ "bucket", if that's relevant.
Is there a way to model such data in an elegant way in the spirit of Gelman and Hill?
Edited: [I radically re-wrote this question to improve it - the specific question has changed but is exactly in the spirit of the original question.]
