Population-wide county-based data: reasonable to report temporal trends using hierarchical modelling?

I have a population-wide data, including county-level information. Subjects are unequally distributed between the 10 counties in the dataset, resulting in multifold differences. The problem is that crude estimates would depend mostly on large counties. However, would it be reasonable to report country-wide temporal trends using hierarchical modelling y ~ predictor + (1|county)?

I think that using a multilevel model for this task makes a lot of sense. A critical issue is how time enters the model. Most typically it would enter as a "fixed" predictor treated as a continuous variable to estimate a linear association between time passed and the outcome:

In lmer:

m1 <- lmer(y ~ time + (1|county), df)


Depending on your goal for the analysis, you might be interested in whether the time trend varies across counties, in which case you can augment the model to allow for county variation in the linear relation between time and y:

m2 <- lmer(y ~ time + (time|county), df)


m1 is nested within m2 and you can use a likelihood ratio test to determine whether the added complexity of m2 (a random slope for time and the random covariance between the time slopes and the county intercepts) provides a better fit to the data than just a single random intercept for county in m1:

anova(m2, m1)


A completely different direction would be to think of the time effect as being crossed with county such that all counties are affected similarly by some event or characteristics that are tracked in the occasions of measurement. This is called a two-way error components model by economists because there are two random intercepts for different clustering units. Psychologists and others call this a cross-classified model:

m3 <- lmer(y ~ 1 + (1|county) + (1|time), df)


The residual from this model ($$e_{ij}$$) captures any interaction between occasion and county as well as other county effect specific to county$$_i$$ on occasion$$_j$$. This model is less common, but is just as valid, especially if you expect the occasion effect to have similar influences on all counties. Note that this model is not nested within either m1 or m2 so you cannot use likelihood ratio testing to compare it to either of them.

• Thank you so much! Am I understand it correctly, that hierarchical modelling gives equal weight for all counties in the analysis. To put it in other words, the temporal trends of small counties will be taken into account similarly as the trends of bigger counties? Thus, hierarchical modelling somewhat makes a country, where all counties have equally sized? – st4co4 Oct 8 '20 at 7:27
• And one more thing about m2. May we also model it without fixed effects: m2 <- lmer(y ~ (time|county), df)? How does it differ? – st4co4 Oct 8 '20 at 7:53
• Hierarchical modeling gives different weights based primarily on sample size. This is particularly true when you allow a parameter (slope or intercept) to be random or varying. Thus these models are sometimes referred to as partial pooling. – Erik Ruzek Oct 8 '20 at 14:33
• I would not generally recommend fitting a model with a random slope but not the accompanying fixed slope. The fixed slope is the what the random slope is referring to, in terms of deviations for each group off that fixed slope estimate. Without it, they are somewhat rudderless. You can always use anova to test these two models against each other. – Erik Ruzek Oct 8 '20 at 14:34
• Thank you so much! Very useful answers! – st4co4 Oct 8 '20 at 15:36