Ordinal Mixed-Effects Random Forest

• Ordinal Mixed-Effects Random Forest (OMERF) is a framework that models ordinal outcomes in grouped data by integrating random forests with mixed-effects cumulative-link models.
• It employs an iterative back-fitting algorithm to jointly estimate nonlinear fixed effects and group-specific random effects, capturing complex interactions.
• OMERF delivers accurate group-level inference and unbiased variance estimation, validated through simulation studies and a PISA 2022 case study.

The Ordinal Mixed-Effects Random Forest (OMERF) is a statistical framework designed to model ordinal outcomes in hierarchical or grouped data, combining nonparametric random forest techniques with mixed-effects modeling in a cumulative-link setting. By jointly estimating nonlinear fixed effects and group-specific random effects, OMERF enables flexible modeling of complex interactions among covariates while rigorously accounting for clustering and enabling inference at multiple data levels. The method operates via an iterative back-fitting algorithm that alternates between nonparametric regression for fixed effects and likelihood-based estimation for random effects and thresholds, thereby facilitating interpretable group-level inferences alongside high predictive accuracy in nonlinear regimes (Bergonzoli et al., 2024).

1. Model Architecture

OMERF addresses ordinal responses $y_{ij} \in \{1,2,\dots,C\}$ observed for unit $j$ in group $i$, structured within a cumulative-link or proportional-odds framework. The cumulative probability for category $c$,

$\gamma_{ijc} = \Pr(y_{ij} \leq c)$

is linked to a set of fixed and random effects via

$$g(\gamma_{ijc}) = \theta_c - \left\{ f(\mathbf{x}_{ij}) + \mathbf{z}_{ij}^T \mathbf{b}_i \right\}$$

where $g$ is a monotone link function (typically logistic), $\theta_c$ are the ordinal thresholds, $f(\mathbf{x}_{ij})$ is an unknown (potentially nonlinear) fixed-effects function, and $\mathbf{z}_{ij}^T \mathbf{b}_i$ encodes group-specific random effects sampled as $$\mathbf{b}_i \sim \mathcal{N}_Q(\mathbf{0}, \Sigma_b)$$.

A latent process representation $$y^*_{ij} = f(\mathbf{x}_{ij}) + \mathbf{z}_{ij}^T \mathbf{b}_i + \epsilon_{ij}$$ with $\epsilon_{ij} \sim \mathrm{Logistic}(0,1)$ underpins the ordinal formulation, with observed responses determined by thresholding $y^*_{ij}$ at $(\theta_{c-1},\theta_c]$.

Unlike standard models assuming $f(\mathbf{x}) = \mathbf{x}^T \beta$, OMERF uses a regression random forest to model $f$ nonparametrically, thereby capturing nonlinear and interaction effects among predictor variables.

2. Fitting Algorithm

Estimation proceeds via a multi-step back-fitting loop, alternating between regression forests and cumulative-link mixed-model (CLMM) fitting:

1. Initialization: Fit a working ordinal forest or regression to obtain initial linear predictor estimates $\hat{\eta}_{ijc}^{(0)}$; set $\hat{\mathbf{b}}_i^{(0)} = \mathbf{0}$.
2. Back-fitting Iterations: For $t = 1,2,\ldots$ until convergence:
   • Compute working targets:

   $$T_{ij}^{(t)} = \hat{\eta}_{ijc}^{(t-1)} + \mathbf{z}_{ij}^T \hat{\mathbf{b}}_i^{(t-1)}$$

- Regress these continuous targets on $\mathbf{x}_{ij}$ using a random forest, yielding an updated fixed-effects function estimate $\hat{f}^{(t)}$. - Fit a CLMM with offset $\hat{f}^{(t)}(\mathbf{x}_{ij})$ (fixed), estimating group random effects $\hat{\mathbf{b}}_i^{(t)}$ and thresholds $\hat{\theta}^{(t)}$ via Laplace or Gauss–Hermite approximation. - Check convergence (e.g., $$\max_i \| \hat{\mathbf{b}}_i^{(t)} - \hat{\mathbf{b}}_i^{(t-1)} \| < \varepsilon$$).

At convergence, OMERF provides joint estimates of the nonparametric $f$, the random-effects $\hat{\mathbf{b}}_i$, and the ordinal thresholds.

3. Prediction and Inference

Predicting for a new observation in group $i$ with features $(\mathbf{x}_{\rm new}, \mathbf{z}_{\rm new})$ uses the estimated model as: $$\hat{\eta}_{{\rm new},\,c} = \hat{\theta}_c - \left\{ \hat{f}(\mathbf{x}_{\rm new}) + \mathbf{z}_{\rm new}^T \hat{\mathbf{b}}_i \right\}$$

$$\hat{\gamma}_{{\rm new},c} = g^{-1}(\hat{\eta}_{{\rm new},c}), \quad \Pr(y_{\rm new} = c) = \hat{\gamma}_c - \hat{\gamma}_{c-1}$$

Fixed thresholds ($\hat{\theta}_c$) and random effect variance components ($\hat{\Sigma}_b$) are accompanied by approximate standard errors from the CLMM step, permitting Wald-type confidence intervals.

The intraclass correlation (ICC), quantifying group-level latent variance, is estimated as: $$\widehat{\mathrm{ICC}} = \frac{\hat{\sigma}_b^2}{\hat{\sigma}_b^2 + \pi^2 / 3}$$ for logit-link scenarios.

Variable importance from the random forest enables interpretation of influential predictors, and partial-dependence plots illustrate the effect of covariates on the underlying latent variable.

4. Hyperparameter Selection

Key tunable parameters include the number of trees $K$ (typically 200–1000), number of variables per split (mtry), minimum node size, maximum tree depth (optional), convergence criterion $\varepsilon$, maximum number of iterations, and random-effects structure (choice of $\mathbf{z}_{ij}$). Hyperparameters can be optimized via out-of-bag error for the forest, cross-validation of the OMERF loop, or information criteria such as AIC or BIC at the CLMM step. Commonly, random forest settings are fixed to defaults while random-effects structure is selected more carefully.

5. Simulation Study

A simulation study assessed OMERF against competitors on synthetic datasets with $I = 10$ groups of $n_i = 100$ units and $C = 3$ ordinal categories, generated from latent models with both polynomial and tree-structured fixed effects, and either random intercept or random intercept plus random slope group effects. The R function genOrdCat() mapped latent outcomes to ordinal responses.

Competing methods included:

• Cumulative Link Model (CLM)
• Cumulative Link Mixed Model (CLMM)
• Ordinal Random Forest (without random effects)
• OMERF

For each scenario, data were split 80/20 for train/test, with evaluation via accuracy, MSE (numeric label surrogate), Adjusted Rand Index, Cohen’s $\kappa$, and Cardoso’s and Ballante’s ordinal-specific indices.

Results showed that OMERF attained the best test-set scores in highly nonlinear generative processes; CLM/CLMM dominated in strictly linear cases. OMERF produced nearly unbiased random-effects estimates, and partial dependence plots from $\hat{f}$ tracked latent fixed-effects surfaces closely.

6. Application: PISA 2022 Case Study

OMERF was applied to data from 7,639 Italian 15-year-old students nested within 293 schools, predicting math proficiency categorized into low, mid, and high (three ordered groups). Fixed-effect covariates encompassed demographics, ICT and school climate variables, personality traits, study behaviors, parental education (HISCED), and socio-economic status (ESCS). A school-level random intercept was modeled.

Implementation employed randomForest::randomForest() for the forest step and ordinal::clmm() for the mixed model, with a convergence threshold of $0.05$ and up to 100 iterations.

Predictive performance across methods is summarized as follows:

Method	Accuracy	MSE	ARI	$\kappa$
CLM	0.594	0.457	0.074	0.248
CLMM	0.670	0.334	0.188	0.396
OrdForest	0.601	0.438	0.091	0.260
OMERF	0.644	0.373	0.150	0.355

CLMM achieved the best overall scores, reflecting the dominance of the nearly linear ESCS predictor. OMERF nevertheless enabled detection of nonlinear influences (notably in personality and affective scales), provided school-level variance decomposition (estimated $\hat{\sigma}_b^2 = 1.695$, ICC ≈ 0.34), and produced highly concordant school effects with CLMM.

Variable importance analyses placed ESCS and HISCED at the top, followed by personality attributes. Partial dependence plots, e.g., for curiosity or emotional control, revealed discrete threshold effects not captured in the linear model. Both OMERF and CLMM indicated substantial between-school heterogeneity.

7. Summary and Context

OMERF integrates the nonparametric flexibility of random forests with the mixed-effects methodology essential for analyzing ordinal, hierarchical data, via a stable and interpretable back-fitting algorithm. It delivers state-of-the-art prediction in nonlinear settings, unbiased random-effects estimation, and multilevel inference, with straightforward implementation in R. Empirical and theoretical analyses demonstrate its superiority over purely parametric or forest-only approaches when complex nonlinearities or interactions are present, while retaining interpretability via variance components and group-level decomposition (Bergonzoli et al., 2024).