1
$\begingroup$

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.]

$\endgroup$
  • $\begingroup$ Do you want something like this: stats.stackexchange.com/questions/134995/… ? Can explain a bit more the issue that you have with this? $\endgroup$ – Martijn Weterings Jan 18 at 11:52
  • $\begingroup$ That captures the nested random effects, which I was already aware of how to do. I am interested in modeling, using random effects or a more Bayesian setting, nested effects + an effect which is not-nested in the nested-variables. I gave the example of trying to model the random effects of geographical location (nested as country, state, county and city effects) plus some other random effect, like age. I'm looking first and foremost for an intuitive and practical approach in this setting. $\endgroup$ – Lepidopterist Jan 18 at 14:30
  • $\begingroup$ "nested effects + an effect which is not-nested in the nested-variables." where is the problem with this? $\endgroup$ – Martijn Weterings Jan 18 at 17:34
  • $\begingroup$ What sort of Bayesian regression/modelling are you doing? This is unclear because you refer to lmer. $\endgroup$ – Martijn Weterings Jan 18 at 17:36
  • $\begingroup$ What lmer is doing can be specified in a fully Bayesian way if you want. Indeed that is how I am able to understand what lmer is doing given Gelman and Hill's exposition. I'm indicating that if lmer cannot do this, I am open to a Bayesian specification that makes sense. With regard to "where is the problem with this", it's that I don't know how to specify a model with partially nested and non-nested random effects. For example, naively I would guess lmer (y ~ x.centered + (1 + x.centered | eth1/eth2/eth3) + (1 + x.centered | age) + (1 + x.centered | eth1/eth2/eth3:age)). $\endgroup$ – Lepidopterist Jan 18 at 17:44
4
+200
$\begingroup$

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().

Yes, you can specify a model with a nested / partially nested structure. Provided that you encode the nested/crossed/partially-crossed factors uniquely, then lmer will handle the structure correctly if you specify the grouping variables as

(1 | eth1) + (1 | eth2) + (1 | eth3) + (1|eth3:age) + (1|age) 

Of course you can add random slopes for x-centered as necessary (It may not be necessary or warranted to fit random slopes to all the factors of the nesting hierarchy).

A nice little trick, if you want to estimate the same model in a Bayesian framework, and/or you are having problems with convergence, is simply to replace the call to lmer with a call to stan_lmer from the rstanarm package, which was developed by Andrew Gelman and colleagues. This is designed as a plug-and-play replacement for lmer which fits a Bayesian model with default priors and estimates it using Hamilton Monte Carlo. If you seek to set you own priors, then this can also be done quite easily. There is a lot of documentation at https://mc-stan.org


Edit (to address the point raised in the first comment below):

This random effects structure:

(1 | eth1/eth2/eth3)

Expands to:

(1 | eth1) + (1 | eth1:eth2) + (1 | eth1:eth2:eth3)

However, provided that the nesting is explicit in the coding of factors, this is also equivalent to :

(1 | eth1) + (1 | eth2) + (1 | eth3)

A simple simulation will demonstrate this:

> set.seed(15)
> dtA <- expand.grid(eth1 = c(1,2,3,4,5), eth2 = c(1,2,3,4,5), eth3 = c(1,2,3,4,5,6), measure = c(1,2))

> dtA$y <- dtA$eth1 + dtA$eth2 + dtA$eth3 + dtA$measure + rnorm(nrow(dtA),0,1)

> fm01 <- lmer(y ~ 1 + (1|eth1/eth2/eth3), data = dtA)
> summary(fm01)  

Random effects:
 Groups           Name        Variance Std.Dev.
 eth3:(eth2:eth1) (Intercept) 3.038    1.743   
 eth2:eth1        (Intercept) 1.718    1.311   
 eth1             (Intercept) 1.948    1.396   
 Residual                     1.377    1.173   
Number of obs: 300, groups:  eth3:(eth2:eth1), 150; eth2:eth1, 25; eth1, 5

However, the "nesting" here is not explicit. Each level of eth2 occurs in every level of eth1:

> xtabs(~ eth1 + eth2, dtA)  
    eth2
eth1  1  2  3  4  5
   1 12 12 12 12 12
   2 12 12 12 12 12
   3 12 12 12 12 12
   4 12 12 12 12 12
   5 12 12 12 12 12

Thus it is ambiguous whether these data are crossed or nested. Nesting is a property of the experimental design. So, if the data are nested, then eth2 = 1 in eth1= 1 is not the same unit of measurement as eth2 = 1 in eth1= 2. If the data are nested and the factors are not coded uniquely between clusters, then it is necessary to write the random effect structure as (1|eth1/eth2/eth3) or (1|eth1) + (1|eth1:eth2) + (1|eth1:eth2:eth3) which tells lmer that the data are nested.

If we fit this model:

> fm02 <- lmer(y ~ 1 + (1|eth1) + (1|eth2) + (1|eth3), data = dtA)

...then we obtain different estimates for the variance components, because unless we specify the nesting explicitly, it means that the random effects are crossed, and not nested:

> summary(fm02)  
Random effects:
 Groups   Name        Variance Std.Dev.
 eth3     (Intercept) 3.132    1.770   
 eth2     (Intercept) 2.232    1.494   
 eth1     (Intercept) 2.395    1.547   
 Residual             1.291    1.136   
Number of obs: 300, groups:  eth3, 6; eth2, 5; eth1, 5

Now, let us encode the factors uniquely:

> dtA$eth2.u <- paste(dtA$eth1, dtA$eth2, sep=".")
> dtA$eth3.u <- paste(dtA$eth1, dtA$eth2, dtA$eth3, sep=".")
> xtabs(~ eth1 + eth2.u, dtA)
    eth2.u
eth1 1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 4.1 4.2 4.3 4.4 4.5 5.1 5.2 5.3 5.4 5.5
   1  12  12  12  12  12   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   2   0   0   0   0   0  12  12  12  12  12   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   3   0   0   0   0   0   0   0   0   0   0  12  12  12  12  12   0   0   0   0   0   0   0   0   0   0
   4   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  12  12  12  12  12   0   0   0   0   0
   5   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  12  12  12  12  12

and now we see that the nesting is explicit.

So to complete the demonstration, we fit:

> fm03 <- lmer(y ~ 1 + (1|eth1/eth2.u/eth3.u), data = dtA)  

or equivalently:

> fm04 <- lmer(y ~ 1 + (1|eth1) + (1|eth2.u:eth1) + (1|eth3.u:eth2.u:eth1), data = dtA)

then we obtain:

> summary(fm03)
Random effects:
 Groups               Name        Variance Std.Dev.
 eth3.u:(eth2.u:eth1) (Intercept) 3.038    1.743   
 eth2.u:eth1          (Intercept) 1.718    1.311   
 eth1                 (Intercept) 1.948    1.396   
 Residual                         1.377    1.173   
Number of obs: 300, groups:  eth3.u:(eth2.u:eth1), 150; eth2.u:eth1, 25; eth1, 5

And finally, if we fit:

> fm05 <- lmer(y ~ 1 + (1|eth1) + (1|eth2.u) + (1|eth3.u), data = dtA)

we obtain:

> summary(fm03)
Random effects:
 Groups   Name        Variance Std.Dev.
 eth3.u   (Intercept) 3.038    1.743   
 eth2.u   (Intercept) 1.718    1.311   
 eth1     (Intercept) 1.948    1.396   
 Residual             1.377    1.173   
Number of obs: 300, groups:  eth3.u, 150; eth2.u, 25; eth1, 5

..which indeed is exactly the same as fm01, fm03 and fm04.

The main take-home point here is that when factors are coded uniquely between clusters, lmer will take care of the nesting, and the model can be specified in several ways, depending on what is most convenient, but if the factors are not unique then the nesting must either be made explicit by changing the coding of the factors, as we did above, or the nesting must be specified using / or : when writing the grouping variables in the call to lmer.

$\endgroup$
  • 1
    $\begingroup$ Well that was rather simple. I think I was hung up on the interaction between the nested variables and the non-nested variables, but I suppose it makes sense to only consider eth3:age. What confuses me is that people commonly code nested variables as (1 | eth1/eth2/eth3) = (1 | eth1) + (1 | eth1:eth2) + (1 | eth1:eth2:eth3), no? Yet in your specification there is no interaction. Can you clarify why that is? Your simpler formulation makes more sense to me, actually, but why the different practices? $\endgroup$ – Lepidopterist Jan 23 at 1:12
  • 1
    $\begingroup$ @Lepidopterist I have updated the answer to address this. $\endgroup$ – Robert Long Jan 23 at 11:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.