# Add extra level to multilevel model that was not part of the sampling process

Consider a population of students, clustered within schools. We are interested on explaining results of a math test at the student level. Assume we use a multistage sampling process in order to capture this two-level structure in our dataset. Thus, in this case, we follow a two-stage random sampling (school first, and then students among these schools). The resulting dataset can be used in a multilevel model, or mixed model like the following:

$$y_{ij} = \beta_{0} + \beta_{1}X_{ij} + \zeta_{i} + \epsilon_{ij}$$

where $i$ is school id and $j$ is student id. $X_{ij}$ includes regressors that vary across schools, students, or both.

Say I have data on the region in which the school is located. Can I add region as an extra, top level to my model? This is, if we index region as $k$, the new model is:

$$y_{ijk} = \beta_{0} + \beta_{1}X_{ijk} + \zeta_{i} + \xi_{ij} + \epsilon_{ijk}$$

Even though the sampling was not originally clustered among regions, could I treat it as if it were? Does this require a test to see if schools are randomly distributed across regions? What if I have all regions in my sample?

The reason why I am thinking on adding an extra level is to get rid of the region dummies and interacting regressors in the two-level equation.

• What I always would consider is that in multilevel models, the clusters are assumed to be random draws from a infinite large population of these clusters. For schools this assumptios is defenitly met imo. For regions I would say, that it depends on how many regions you have (20 regions would be too less imo but 200 would be fine). – Qaswed Jul 6 '16 at 17:37
• But it is also common to have levels where all instances are observed. For example, if you have a time dimension and measure every students during a given period of years. Thus, is randomness of observed cluster within population of clusters really required? Or this example work only because time would be the bottom level? – luchonacho Jul 6 '16 at 17:44
• I'm not shure if I got your example with the measurement of a variable of students during a time period. But in this case, I would say, that each observation in the period is nested in one student - so each student builds his own cluster. – Qaswed Jul 6 '16 at 18:01

The short answer is: Yes, you can. The fact that you didn't originally include region in your sampling is not very important.

(Additional 2 paragraphs to address the points in comments.) If schools are clustered within with regions, then it is important to account for non-independence within clusters. This can be done by including random intercepts for region, with schools within region, thus creating an additional level of cluttering, or it can be done by including region as a fixed effect.

The only justification that is needed for including a higher level is comes down to whether you have enough regions. There are also some other issues and I will address these below. Not including it as a level or a fixed effect could result in biased estimates and invalid inference. One of your comments asked why region is not included as a level in other studies. There could be several reasons: a) the region data was not available, b) there was too much data and the model could not be run (since adding further levels increases the computational burden), c) region was added but it turned out to be not needed since there was very little variation in the response at the region level.

Although you may have all regions in your sample, is this really the whole population of regions ? What about regions in other countries ? There is no universal agreement over treating a factor as fixed or random. In fact there is considerable disagreement, which is nicely summed up by Andrew Gelman Here:

(1) Fixed effects are constant across individuals, and random effects vary. For example, in a growth study, a model with random intercepts a_i and fixed slope b corresponds to parallel lines for different individuals i, or the model y_it = a_i + b t. Kreft and De Leeuw (1998) thus distinguish between fixed and random coefficients.

(2) Effects are fixed if they are interesting in themselves or random if there is interest in the underlying population. Searle, Casella, and McCulloch (1992, Section 1.4) explore this distinction in depth.

(3) “When a sample exhausts the population, the corresponding variable is fixed; when the sample is a small (i.e., negligible) part of the population the corresponding variable is random.” (Green and Tukey, 1960)

(4) “If an effect is assumed to be a realized value of a random variable, it is called a random effect.” (LaMotte, 1983)

(5) Fixed effects are estimated using least squares (or, more generally, maximum likelihood) and random effects are estimated with shrinkage (“linear unbiased prediction” in the terminology of Robinson, 1991). This definition is standard in the multilevel modeling literature (see, for example, Snijders and Bosker, 1999, Section 4.2) and in econometrics.

Even when you truly have the entire population of regions, there are still arguments in favour of treating the factor as random:

1. The model will be more parsimonious
2. Variance partitioning across the levels is much easier.

While there is no need to test whether schools are randomly distributed within regions, depending on which software you use to fit the model there is an assumption that the random effects themselves are normally distributed.

If the number if regions is very low, this would be an argument in favour of treating region as a fixed effect:

Treating factors with small numbers of levels as random will in the best case lead to very small and/or imprecise estimates of random effects; in the worst case it will lead to various numerical difficulties such as lack of convergence, zero variance estimates, etc.. (A small simulation exercise shows that at least the estimates of the standard deviation are downwardly biased in this case; it's not clear whether/how this bias would affect the point estimates of fixed effects or their estimated confidence intervals.) In the classical method-of-moments approach these problems may not arise (because the sums of squares are always well defined as long as there are at least two units), but the underlying problems of lack of power are there nevertheless. (http://glmm.wikidot.com/faq)

The question of what is the minimum number of levels is again a matter of much debate. Frequentist approaches usually call for a minimum number ranging from 5 (for example see here) to 20 (for example see Rabe-Hesketh and Skrondal, 2012, p124). Bayesian approaches on the other hand, allow as few as 3 (see Gelman, 2006)

A good approach in your situation would be to fit a model with region as a random effect, and another with region as a fixed effect and compare them.

References:

Gelman, A "Prior distributions for variance parameters in hierarchical models", Bayesian Anal. Volume 1, Number 3 (2006), 515-534.

Green, B. F., and Tukey, J. W. (1960). Complex analyses of variance: general problems. Psychometrika 127–152.

Kreft, I., and De Leeuw, J. (1998). Introducing Multilevel Modeling. London: Sage.

LaMotte, L. R. (1983). Fixed-, random-, and mixed-effects models. In Encyclopedia of Statistical Sciences, ed. S. Kotz, N. L. Johnson, and C. B. Read, , 137–141.

Rabe-Hesketh, S and Skrondal, A, Multilevel and Longitudinal Modeling Using Stata, Third Edition. Volume I: Continuous Responses, 2012, ISBN-13: 978-1-59718-108-2

Searle, S. R., Casella, G., and McCulloch, C. E. (1992). Variance Components. New York: Wiley.

Snijders, T. A. B., and Bosker, R. J. (1999). Multilevel Analysis. London: Sage.

• At first it seems to me that most of your answer (except first two paragraphs) is about whether a factor should be treated as fixed or random. Yet, my question is about adding an extra level, even when sampling method do not cluster at that level. Can you please clarify why this would be related to how to treat a factor? I cannot see the relationship. More precisely, I cannot see how the discussion supports the yes, you can statement in the first paragraph. – luchonacho Jul 12 '16 at 12:34
• As I said, it doesn't matter whether your original sampling included the region clusters. The factor I am referring to is the region ID. Since schools are clustered within regions, you can either use it as a higher grouping variable (extra level) or a fixed effect. – Robert Long Jul 12 '16 at 12:44
• In case I wasn't clear, using region as a random effect (random intercept) means including it as a higher level, with school clustered in it. – Robert Long Jul 12 '16 at 12:51
• So, it seems you are telling me yes, you can [use region as an extra level] without explaining why. And then you go into how to deal with region as a factor (not my interest), focusing on using region as either RE or FE. Please correct me if I'm wrong. I want more evidence that adding region is not arbitrary. Is there a test that can be done? – luchonacho Jul 12 '16 at 13:14
• You can use region as a random effect (and therefore as a higher "level") so long as schools are nested in region. Either they are nested or they arent. You don't say how many levels (how many regions) you have or what software you are using (hence my discussion about the number of levels necessary). Assuming you have sufficient levels then only after you fit the model will you know if it makes sense (if there is sufficient variation at the region level to warrant it's inclusion). – Robert Long Jul 12 '16 at 15:28