Fitting multilevel models with ordinal outcomes: performance of alternative specifications and methods of estimation

Abstract

Previous research has compared methods of estimation for fitting multilevel models to binary data, but there are reasons to believe that the results will not always generalize to the ordinal case. This article thus evaluates (a) whether and when fitting multilevel linear models to ordinal outcome data is justified and (b) which estimator to employ when instead fitting multilevel cumulative logit models to ordinal data, maximum likelihood (ML), or penalized quasi-likelihood (PQL). ML and PQL are compared across variations in sample size, magnitude of variance components, number of outcome categories, and distribution shape. Fitting a multilevel linear model to ordinal outcomes is shown to be inferior in virtually all circumstances. PQL performance improves markedly with the number of ordinal categories, regardless of distribution shape. In contrast to binary data, PQL often performs as well as ML when used with ordinal data. Further, the performance of PQL is typically superior to ML when the data include a small to moderate number of clusters (i.e., ≤ 50 clusters).

Figure 1

Marginal category distributions used in the simulation study (averaged over predictors and random effects). Notes. Within 3-, 5-, and 7-category outcome conditions, marginal frequencies are held constant but permuted across categories to manipulate the distribution shape (bell-shaped, skewed, or polarized) without changing sparseness. Within 2-category outcome conditions, it is impossible to hold marginal frequencies constant while manipulating shape (balanced or unbalanced).

Figure 2

Average bias for the three fixed effect estimates (excluding thresholds) across estimator, number of outcome categories, and distribution shape. Notes. The normal-theory REML (Restricted Maximum Likelihood) estimator was used when fitting the linear multilevel model. The estimators of PQL (Penalized Quasi-Likelihood) or ML (with adaptive quadrature) were used when fitting the multilevel cumulative logit (logistic) model. Points for two-category conditions are not connected to points for 3-7 category conditions because their distribution shapes do not correspond. Results show that bias is large and sensitive to distribution shape when using the linear model but not when using the cumulative logit model (either estimator). Results are collapsed over the number of clusters, cluster size, and the magnitude of the random effects.

Figure 3

Average bias for the three fixed effect estimates (excluding thresholds) across logistic estimators, number of outcome categories and cluster size. Notes. Logistic estimators were either PQL (Penalized Quasi-Likelihood) or ML (Maximum Likelihood) with adaptive quadrature. Results show that PQL produces somewhat negatively biased fixed effect estimates, particularly when random effects have large variances, whereas the estimates obtained from logistic ML show small, positive bias. In both cases, bias decreases with the number of categories of the outcome. Results are collapsed over number of clusters and distribution shape and do not include linear multilevel model conditions.

Figure 4

Mean-Squared Error (MSE) for the fixed effects (excluding thresholds) across number of outcome categories, number of clusters, and cluster size. Notes. MSE is indicated by the height of the vertical lines, and it is broken into components representing squared bias (portion of the line below the symbol) and sampling variance (portion of the line above the symbol). The scale differs across panels and is discontinuous in the upper right panel. MSE is averaged across the three fixed effects.Results are plotted for multilevel cumulative logit models; PQL denotes Penalized Quasi-Likelihood and ML denotes Maximum Likelihood with adaptive quadrature. This plot does include linear multilevel model conditions.

Figure 5

Average bias for the standard errors (SE) of the three fixed effect estimates (excluding thresholds) across number of outcome categories and cluster size. Notes. Results are plotted for multilevel cumulative logit models; PQL denotes Penalized Quasi-Likelihood and ML denotes Maximum Likelihood with adaptive quadrature. This plot does include linear multilevel model conditions.

Figure 6

Mean-Squared Error (MSE) for the standard deviation of the random intercept, across number of outcome categories, number of clusters, and cluster size. Notes. The scale differs across panels and is discontinuous in the upper right panel. See Figure 4 notes for definition of quantities in this plot.

Figure 7

Mean-Squared Error (MSE) for the standard deviation of the random slope, across number of outcome categories, number of clusters, and cluster size. Notes. The scale differs across panels and is discontinuous in the upper right panel. See Figure 4 notes for definition of quantities in this plot.