Small Sample Methods for Multilevel Modeling: A Colloquial Elucidation of REML and the Kenward-Roger Correction

ABSTRACT

Studies on small sample properties of multilevel models have become increasingly prominent in the methodological literature in response to the frequency with which small sample data appear in empirical studies. Simulation results generally recommend that empirical researchers employ restricted maximum likelihood estimation (REML) with a Kenward-Roger correction with small samples in frequentist contexts to minimize small sample bias in estimation and to prevent inflation of Type-I error rates. However, simulation studies focus on recommendations for best practice, and there is little to no explanation of why traditional maximum likelihood (ML) breaks down with smaller samples, what differentiates REML from ML, or how the Kenward-Roger correction remedies lingering small sample issues. Due to the complexity of these methods, most extant descriptions are highly mathematical and are intended to prove that the methods improve small sample performance as intended. Thus, empirical researchers have documentation that these methods are advantageous but still lack resources to help understand what the methods actually do and why they are needed. This tutorial explains why ML falters with small samples, how REML circumvents some issues, and how Kenward-Roger works. We do so without equations or derivations to support more widespread understanding and use of these valuable methods.

KEYWORDS:

• Restricted maximum likelihood
• Kenward-Roger
• tutorial
• explanation
• mixed model

Article information

Conflict of Interest Disclosures: The author signed a form for disclosure of potential conflicts of interest. The author did not report any financial or other conflicts of interest in relation to the work described.

Ethical Principles: The author affirms having followed professional ethical guidelines in preparing this work. These guidelines include obtaining informed consent from human participants, maintaining ethical treatment and respect for the rights of human or animal participants, and ensuring the privacy of participants and their data, such as ensuring that individual participants cannot be identified in reported results or from publicly available original or archival data.

Funding: This work was not supported by any external funding.

Role of the Funders/Sponsors: None of the funders or sponsors of this research had any role in the design and conduct of the study; collection, management, analysis, and interpretation of data; preparation, review, or approval of the manuscript; or decision to submit the manuscript for publication.

Acknowledgments: The author would like to thank Dan Bauer, Patrick Curran, and Dave Thissen for their comments on prior versions of this manuscript. The ideas and opinions expressed herein are those of the authors alone, and endorsement by the author's institution is not intended and should not be inferred.

Notes

1 There are some minor differences between these types of models but, for the purpose of this manuscript, we consider the models to be more or less interchangeable.

2 Technically, when implemented in software, this step is conducted with an error contrast of the outcome variable and the projection matrix rather than actually fitting a single-level model with OLS and saving the residuals. The result of the process has the same effect, so we describe the process using OLS to increase the intuition of the process and minimize the reliance on mathematical terminology.

3 In MLMs, the total variance V is equal to V = ZτZ^T + R where τ is the covariance matrix of the Level-2 random effects, Z is the random effect design matrix, and R is the covariance matrix of the Level-1 residuals. Using the REML estimates of τ and R to compute V, the GLS fixed effects are computed by (X^TV^{− 1}X)^{− 1}X^TV^{− 1}y where X is the fixed effect design matrix and y is a vector of the outcome variable values. The addition of the V term accounts for the dependence of observations within clusters which is ignored by OLS. The OLS fixed effect estimates are calculated by (X^TX)^{− 1}X^Ty.

4 The GLS and ML estimates are equal on average but are not perfectly identical. That is, in the long run, GLS and ML are expected to yield the same values but there may be small differences when estimates are compared for a single data set.

5 The full mathematical expression in Equation 4(4) is much more complex than the function y = e^x used in Figure 3, so we do not wish to imply that the a Taylor series expansion of Equation 4(4) will necessarily result in the remarkably close approximation seen in Figure 3. The example in Figure 3 was chosen because it is a textbook example of how a Taylor series expansion operates. As a result, the Taylor series expansion happens to be quite good.

6 We refer to univariate tests here, such as testing a single coefficient. If multiparameter tests are desired, the statement should be revised to say that researchers should use F tests rather than χ² tests because the F distribution is a finite sample version of the χ² distribution (i.e., an F distribution with infinite denominator degrees of freedom is equal to a χ² distribution if the numerator F degrees of freedom match the χ² degrees of freedom).

7 Some software programs such as Stata and Mplus report Z-tests for fixed effects which assume infinitely large samples and therefore do not have degrees of freedom. With smaller samples, Z-tests are not appropriate, generally speaking.

8 The Prasad-Rao-Jeske-Kackar-Harville is not available as a formal option in SAS though the relevant information can be pieced together. Using the DDFM=KR(FIRSTORDER) option in the MODEL statement uses the Prasad-Rao-Jeske-Kackar-Harville method for the standard error correction rather than the Kenward-Roger correction. The output still contains the Kenward-Roger degrees of freedom, so we manually calculated the p-values using the containment degrees of freedom.

9 Containment degrees of freedom are calculated by N − rank(X, Z) where N is the overall sample size, X is a design matrix for the fixed effects and Z is a design matrix for the random effects. In linear algebra, the rank of a matrix is the dimension of the vector space spanned by the columns. The computation can be simplified with a random intercepts model for balanced data (as is the case in the Stapleton et al. data). For a fixed effect that does not have a random effect, the containment degrees of freedom are equal to (total sample size) minus (number of predictors) minus (number of clusters minus Level-2 predictors). This model has 84 total people, 3 predictors, 14 clusters, and 2 Level-2 predictors, so the containment degrees of freedom for a nonvarying fixed effect is 84 − 3 − (14 − 2) = 69. For effects that do have random effects (the intercept in this model), the containment degrees of freedom is equal to (number of clusters) minus (number of Level-2 random effects) minus (number of Level-2 predictors). In this data, there are 14 clusters, 1 Level-2 random effect, and 2 Level-2 predictors, so the degrees of freedom for the intercept is 14 − 1 − 2 = 11. These formulas are not as straightforward for unbalanced data or models with random slopes because the output from the rank operator will be dependent on aspects like the degree of unbalancedness and the random effect covariance structure.