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:
- The model will be more parsimonious
- 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.
