Research

A permutation-based random forest diagnostic for linear mixed models

Independent research · Sole author · Department of Statistics, University of Washington · 2025

A linear mixed model can be wrong in several ways at once, and the standard diagnostics are not well suited to finding out which. RandomForestSpecCheck takes a different route: it asks a random forest to predict the fitted model's residuals from the covariates, and uses a permutation test that respects the clustered design to decide whether the forest's success is real. The result is a single verdict on whether the model has left systematic structure behind, without any commitment in advance to the form that structure might take.

The problem with checking a mixed model

Linear mixed models are the standard instrument for correlated data: repeated measures on the same subject, students inside classrooms, patients inside clinics. Their appeal is that they model the correlation rather than ignoring it, and the price of that appeal is that several assumptions have to hold at once. The fixed effects must have the right functional form. The random effects must have the right structure, which usually means deciding whether a slope varies across clusters and not only an intercept. The residual covariance must be specified correctly. A failure in any one of these can bias estimates or, more often and more insidiously, leave the estimates roughly intact while corrupting the standard errors.

The conventional checks are the ones inherited from ordinary regression. Residuals are plotted against fitted values and inspected for curvature. Normal quantile plots are consulted for the random effects and the within-cluster errors. Information criteria are compared across a handful of candidate specifications. Each of these is useful and each has a characteristic blind spot. A residual plot is a judgement call, and a dense cloud of several thousand points will absorb a good deal of structure before it looks anything other than round. Information criteria compare only the models an analyst thought to fit, so a misspecification that occurred to nobody is never in the comparison at all. Formal tests for mixed-model adequacy exist, but they tend to be aimed at one alternative at a time, which returns the analyst to the same difficulty: you have to guess what is wrong before you can test for it.

The statistic

Fit the mixed model, extract the conditional residuals, and regress them on a set of candidate predictors using a random forest. The predictors are the ones an analyst suspects might have been mishandled: time, the cluster identifier, a covariate that entered only linearly. The measure of the forest's success is its out-of-bag coefficient of determination,

R²_oob = 1 − MSE_oob / var(ê_ij)

where each tree is grown on a bootstrap sample and evaluated only on the observations that sample left out. Aggregating those out-of-bag predictions gives the forest an internal cross-validation, and so a nearly unbiased estimate of how well it generalises, without setting aside a test set. This matters more here than in an ordinary prediction problem, because a forest with enough depth will fit residual noise perfectly if permitted to score itself on the data it was grown from, and the diagnostic would then flag every model ever fitted.

The framing is a familiar hypothesis test wearing unfamiliar clothes. Under correct specification the conditional expectation of a residual given any predictor is zero, and a sufficiently large forest approximates conditional expectations, so the best achievable out-of-bag R² is zero up to sampling variability. Under misspecification the forest has something to find. The null and alternative are

H₀: E[e | X] = 0   versus   H₁: E[e | X] ≠ 0

tested without specifying any parametric alternative. That is the whole point of using a forest rather than an expanded regression: the analyst supplies the variables, not the functional form.

Where the null distribution comes from

The out-of-bag R² of a random forest has no tractable sampling distribution, so it has to be simulated. The obvious device is a permutation test, and the obvious hazard is that the data are clustered. Shuffling residuals freely across the whole dataset would destroy the correlation structure along with the signal, and the resulting null would describe a scenario the data were never drawn from.

The procedure therefore permutes within clusters. Under the null hypothesis of no association between residuals and predictors, residuals are exchangeable inside a cluster: any ordering is as likely as any other, so permuting them leaves the joint distribution alone while severing the link to the predictors. Each cluster keeps its size and its membership; only the internal ordering is scrambled. Refitting the forest to the permuted residuals and repeating B times gives a distribution of out-of-bag R² values under a design identical to the real one, in a world with nothing to find. The test is valid in finite samples under exchangeability, with no appeal to asymptotics and no assumption that residuals are normal.

The decision rule

Because large values of the statistic indicate leftover structure, the test is one-sided and the comparison is against the upper tail. The procedure takes the 97.5th percentile of the permutation distribution as its critical value. On its own, however, statistical significance is a weak standard for a diagnostic: with enough data, structure that explains a fraction of a percent of residual variance will clear any percentile, and an analyst who redesigns a model on that basis has been badly served.

The procedure therefore applies two criteria and flags a model only if both are met:

R²_oob > Q_0.975   and   R²_oob > τ

The second criterion is a practical threshold, set by default at τ = 0.05, which is to say that the leftover structure must account for at least five percent of residual variance before it is worth an analyst's attention. The default is a judgement rather than a derivation, and it is exposed as an argument for exactly that reason. In the simulations below it balanced sensitivity against specificity well; in the real-data application it turns out to do essentially all of the discriminating work, which is a finding worth stating plainly rather than burying.

The diagnostic, running

The panel below is a working implementation of the procedure in the browser, not an illustration of it. It simulates clustered longitudinal data from a model with a quadratic time trend, fits the random-intercepts model that omits the quadratic term, grows a forest on the resulting residuals, builds the permutation null by shuffling those residuals within clusters, and applies the rule. Nothing is precomputed.

Running.

What the method detects, and what it does not

Four simulation studies were run, each generating data under a specific violation and fitting a model that omitted it. Three concern the mean structure: an omitted quadratic fixed effect, an omitted random slope, and unmodelled AR(1) correlation among the within-cluster errors. The fourth concerns the variance: residual variance increasing with time, fitted with a model assuming homoscedasticity. Every cell of every design was replicated one hundred times, and the quantity reported is the proportion of replicates in which the diagnostic flagged the model.

The three mean-structure studies behave as a diagnostic should. False positives are held between zero and three percent throughout, which is conservative against a nominal five. Power then depends on the two things it ought to depend on, the size of the omitted effect and the signal-to-noise ratio of the design. A strong quadratic trend is caught in eighty to one hundred percent of replicates. A weak one, contributing a coefficient of 0.05, is caught rarely, and the honest reading is that the method requires either a large sample or a quiet one to detect subtle nonlinearity. Omitted random slopes show the same pattern with an added dependence on cluster size: with ten observations per cluster and low noise the diagnostic finds a slope variance of 0.4 in eighty-one to ninety-six percent of replicates, and with five observations per cluster it finds almost nothing, because the temporal information needed to reveal a slope simply is not there. Unmodelled autocorrelation is the hardest of the three, reaching only about a quarter of replicates at ρ = 0.6.

The fourth study is the one worth dwelling on. Heteroscedastic residuals were essentially never flagged, at any variance strength, in any design. This is not a failure of the implementation; it is a consequence of what the statistic measures. The forest is predicting the residual, not its magnitude, and in a design where the variance grows with time while the mean stays at zero there is no conditional mean for the forest to recover. From the point of view of the statistic the fitted model is adequate, and the statistic says so.

Application: the Framingham Heart Study

The diagnostic was applied to a subset of the Framingham Heart Study, comprising more than 2,600 cholesterol measurements taken at two-year intervals over ten years, alongside baseline age, sex, and body mass index. Five linear mixed models were fitted. Four are specifications a careful analyst might reasonably choose between: random intercepts, random slopes for time, an AR(1) residual correlation structure, and a quadratic term in time. The fifth is a deliberate act of sabotage, omitting time altogether from the fixed effects, a variable known on biological grounds to matter.

Two things are worth reading off the table. The first is the intended result: the four reasonable specifications were cleared and the model missing a known predictor was flagged, with an observed R² of 0.191 against a threshold of 0.05. Standard residual-versus-fitted plots for all five models produce dense clouds centred at zero, and the sabotaged model's cloud is not visibly different from the others; the quantile plots are likewise unremarkable in all five cases. The misspecification was in the mean structure, which is where the forest was looking and where the eye was not.

The second is less comfortable and more interesting. Every one of the five models, including the four that passed, has an observed R² above its own permutation 97.5th percentile. The permutation criterion alone separated nothing here. What separated the sabotaged model was the practical threshold, and the reason lies in the geometry of the null: because permuted residuals contain no signal at all, the out-of-bag R² under the null is reliably negative, so the upper percentile of the null sits below zero and almost any real fit clears it. In a dataset this size the first criterion is close to free. That is precisely the situation the second criterion was introduced for, and the application shows it earning its place rather than merely accompanying the first. It also sharpens a limitation: on data of this scale the procedure's operating characteristics are governed by τ, and an analyst who accepted the default without thinking about the effect size it encodes would be delegating the decision to a number they had not chosen.

Limitations

The procedure detects mean structure and only mean structure, as the fourth simulation study establishes directly. Power against subtle effects is modest and depends jointly on the number of clusters, the number of observations within them, and the noise level, so a null result in a small or noisy study is weak evidence of adequacy. The forest can only find structure in variables it is given, so a misspecification involving a covariate that was never measured is invisible by construction. And because the permutation scheme deliberately preserves clusters, a wholly omitted cluster-level mean effect is not disrupted under the null and will not be flagged; that failure mode is caught by simpler means, such as comparing average residuals across clusters, but it has to be caught by something. The method is a screening instrument for nuanced within-cluster structure, and it is at its most useful when treated as the first question asked rather than the last.

The function

The procedure is packaged as a single R function over ranger for the forests and dplyr for the grouped permutations. It takes a fitted model object, the data, a formula naming the predictors to search, the number of permutations, and the threshold, and it prints a verdict alongside the numbers behind it.

RandomForestSpecCheck(
  fitted_model = my_lmm,
  data         = my_data,
  formula_rf   = resid_model ~ time + cluster,
  B            = 500,
  threshold    = 0.05
)

RandomForestSpecCheck: Random Forest Model Specification Diagnostic for LMMs

Practical threshold:                        0.05
Observed R-squared:                         0.191
Permutation mean R-squared and 95% CI:     -0.009 (-0.011, -0.006)

Result: Model likely misspecified

Interpretation: The observed R-squared exceeds both the upper bound of the
null distribution and the practical threshold, indicating that the model leaves
systematic structure in the residuals. This suggests that the model may be
misspecified.

The printed output is deliberately verbose. A diagnostic that returns a bare logical invites an analyst to treat it as a verdict, and it is not one; it is a signal that some named variables are still doing work the model has not accounted for. The natural next steps are to read the forest's variable-importance measures to see which ones, and then to attempt the corresponding parametric extension and check whether the flag clears.

Contribution

Permutation inference and tree ensembles are both long established, and residual analysis with machine learning has been explored for ordinary and generalised linear models. What had not been assembled was the combination in a mixed-model setting: a within-cluster permutation scheme that respects the nested design, a forest-derived statistic with an interpretable scale, a decision rule that separates significance from magnitude, and a simulation programme that maps where the resulting procedure has power and where it does not, including an explicit negative result. The work was written for submission to Computational Statistics & Data Analysis.

← All statistical work