Multilevel Imputation For Ordinal And Nominal Variables

Why categorical imputation must reflect the hierarchy

A categorical variable observed for people within clinics, schools, or other groups has two features to represent: its category structure and the dependence among observations in the same cluster. An ordinal response, such as a three-level severity rating, has a meaningful order. A nominal response, such as treatment preference, does not. Choosing a model that ignores either the order or the clustering can distort the category probabilities used to fill in missing values.

The imputation model should also distinguish predictors measured for individuals from predictors measured for clusters. Individual-level predictors can vary within a cluster; a cluster-level predictor has one value shared by all individuals in that cluster. If there is a reason to separate within-cluster and between-cluster associations for an individual-level predictor, its cluster mean can be included alongside the individual’s deviation from that mean, as in the earlier material on nonlinear effects within and between clusters.

Model assumptions and category probabilities

Let \(Y_{ij}\) be a categorical variable for individual \(i\) in cluster \(j\). For an ordinal variable, suppose its \(K\) categories are ordered from 1 to \(K\). A cumulative-logit model with a cluster random intercept is

Here \(X_{ij}\) contains individual-level predictors, \(W_j\) contains cluster-level predictors, \(\tau_1<\cdots<\tau_{K-1}\) are ordered cutpoints, and \(b_j\) is a shared cluster effect. A common assumption is \(b_j\sim N(0,\sigma_b^2)\), independently across clusters. Conditional on predictors and \(b_j\), individuals’ outcomes are usually treated as independent. The shared \(b_j\) induces dependence among people in the same cluster.

The model assumes proportional odds: each predictor has the same coefficient at every cumulative split \(Y\leq k\) versus \(Y>k\). With the displayed sign convention, a positive predictor coefficient shifts probability toward higher categories. To obtain the probability of each category, define \(F(t)=1/(1+\exp(-t))\), \(\tau_0=-\infty\), and \(\tau_K=+\infty\). Set \(\eta_{ij}=X_{ij}^{T}\beta+W_j^{T}\gamma+b_j\). Then

For a nominal variable, categories have no natural order. Choose category \(K\) as a reference and use a baseline-category multinomial-logit model:

The vector \(\mathbf{u}_j=(u_{j1},\ldots,u_{j,K-1})^T\) represents shared cluster effects and is commonly modeled as multivariate normal with mean zero and covariance matrix \(\Sigma\). The category probabilities are the softmax probabilities: if \(q_{ijk}=\exp(\alpha_k+X_{ij}^{T}\beta_k+W_j^{T}\gamma_k+u_{jk})\) for \(k<K\), and \(q_{ijK}=1\), then \(P(Y_{ij}=k)=q_{ijk}/\sum_{\ell=1}^K q_{ij\ell}\). Category-specific coefficients allow predictor effects to differ across contrasts with the reference category. Changing the reference changes coefficient descriptions, not the fitted probabilities.

These models require the category order to be substantively meaningful for the ordinal model, a suitable link and predictor specification, and enough observed data to estimate the cutpoints or category-specific coefficients. The random-effects distribution and conditional-independence assumptions should be reasonable. For imputation under a missing-at-random assumption, the conditional model should include observed variables that predict the missing category or its missingness, including variables used in the substantive analysis. If a category is absent or extremely rare in many clusters, cluster-specific probabilities may be poorly determined.

Why the category formula is coherent

For the ordinal model, write \(C_k=P(Y_{ij}\leq k\mid X_{ij},W_j,b_j)\), with \(C_0=0\) and \(C_K=1\). The category events are disjoint, and \(\{Y=k\}=\{Y\leq k\}\setminus\{Y\leq k-1\}\). Therefore \(P(Y=k)=C_k-C_{k-1}\). The ordered cutpoints make \(F(\tau_k-\eta)\) nondecreasing in \(k\), so these differences are nonnegative. They sum to

Thus the cumulative model yields a valid probability distribution over the original categories without fitting a separate binary model for every split. In a multilevel model, these probabilities are conditional on the cluster effect. A shared effect matters: people in the same cluster use the same \(b_j\) or \(\mathbf{u}_j\), which can make their category outcomes associated even when their individual-level predictors differ.

Using predictors from both levels

The predictor specification should reflect how the data were generated and how the substantive analysis will use them. For example, including only an individual-level measurement \(X_{ij}\) can blend its within-cluster and between-cluster associations. A useful specification may include the cluster mean \(\bar{X}_j\) and the deviation \(X_{ij}-\bar{X}_j\) separately, together with a cluster-level predictor \(W_j\). The example below shows how these terms change the ordinal imputation probabilities.

Imputation, alternatives, and diagnostics

In a chained-equations procedure, the relevant conditional model is fitted using observed values and the current values of other variables in the imputation system. For each missing categorical value, the model supplies probabilities for all allowed categories, and a category is drawn from that distribution. The draw should reflect uncertainty in the fitted model and cluster effects; substituting a single fitted most-likely category would suppress variation and understate uncertainty. A cluster-level categorical variable is different: it has one value per cluster, so it should be imputed once at the cluster level and then represented consistently for its members, not independently for each row.

A pooled single-level ordinal or multinomial model is simpler, but treats observations as conditionally independent across clusters and may misrepresent both uncertainty and cluster-specific category patterns. A cumulative-link model uses the ordinal information efficiently when proportional odds is reasonable. A multinomial model avoids that restriction but spends more parameters and ignores ordering. Fitting separate binary models for ordinal thresholds is another possible strategy, but their predicted cumulative probabilities can violate the required ordering unless constraints are imposed. The model choice should be based on the variable’s meaning and diagnostics, not on which model produces the most convenient result.

Useful checks include examining category counts overall and by cluster, inspecting whether imputations reproduce important observed associations, and checking convergence and stability of model estimates. For an ordinal model, assess whether the common-slope proportional-odds assumption is plausible; systematic differences across cumulative splits warn against it. For a nominal model, check whether some category contrasts are nearly unsupported by the data or whether estimates are unstable. Compare observed and imputed category frequencies within relevant predictor and cluster groups. Implausible cluster-specific predictions, persistent computational warnings, or chains that do not stabilize are reasons to revisit the specification, category coding, or available information.

No model can recover a category that is unsupported by the observed data without strong assumptions. Sparse clusters, very rare categories, few clusters, and a misspecified random-effects distribution can make estimates fragile. Do not collapse categories solely to make a model run: collapsing changes the scientific meaning and may discard a meaningful order or distinction. Instead, assess whether the variable can be modeled more parsimoniously, whether predictors can be simplified without omitting essential analysis terms, or whether the limits of the available information must be reported.