# Specifying multilevel model structure when random effects exhaust the population

I have been working with a dataset featuring observations at the county level for about 1300 of the ~3100 or so counties of the United States. These 1300 counties are drawn from every state in the union plus DC.

In my initial models, I estimated relatively simple linear regressions:

lm(y ~ x)

I eventually worked my way up to multilevel models where I allow the intercept and slope to vary across the states. In lmer-style syntax, this looks like this:

lmer(y ~ x + (x|state))

Using a cross-validation model comparison procedure, the latter model is preferred quite strongly. However, I've begun to rethink the wisdom of even estimating this model in the first place. I can think of a number of competing reasons for and against specifying this model:

In it's favor, the observations within states are correlated. The ICC for a random-intercept only version of the model described above is .23. So, the observations are not independent, given the basic linear model.

Also in it's favor, the cross-validation procedure prefers it! In some ways, I think this is the best evidence in favor of using this model.

The main point against it is that by specifying varying intercepts/slopes across states, we are making the assumption that these states are drawn from some population of states. The fixed effects are the expected association across this population (i.e. population effect). The problem is that there is no such population of states. There are 51 (including DC), and we have sampled all of them.

I would appreciate anyone pointing out either where my thinking is incorrect, or indicating some way of resolving these apparent inconsistencies.