A review on longitudinal data analysis with random forest in precision medicine

In plain words

In medicine you often measure the same patient many times: blood samples every few months, several stool samples at one visit, and so on. Random forests are a popular way to build prediction models from such data. They assume every row in the table is independent, though, and rows from one patient are not. If you ignore this, the forest can memorise patients instead of learning patterns, its trees end up too similar to each other, and its built-in error estimate looks better than it really is.

This paper is a survey of the ways people have adapted random forests to data with repeated measurements. The simplest fixes change how data is sampled: average each patient's rows into one, or resample whole patients instead of individual rows. More involved fixes change the model. One family summarises each predictor's recent history into features. Another combines the forest with a classical mixed-effects model, which describes patient-specific offsets explicitly and lets the forest learn the rest. A third family treats a patient's whole series of outcomes as one vector and changes how trees choose splits.

The paper does not run its own experiments. It gives the model equations, explains how each method is fitted, lists the R software available for each, and points out gaps: few methods give variable importance, few have been tested on high-dimensional data, and none have been compared with each other systematically.

TL;DR

A narrative review of CART-based random forest (RF) extensions for longitudinal and clustered data, organised by data structure: univariate response without time effect (averaging, subject-level bootstrapping / RF++), univariate response with time effect (historical RF, LMM/GLMM hybrids such as MERT/MERF, RE-EM tree/REEMforest, SMERF/SREEMforest, GMERT, GMET/GMERF, Bayesian BiMM), and multivariate responses (distance- and likelihood-based split criteria, spline-coefficient responses) §1 §2 §3. It gives the core equations and EM-style fitting loops, a software table tab. 3, and lists open problems §5. It has no experiments of its own and no code. Use it as a map of the methods and their trade-offs. Implement from the original papers, because several formulas here contain typos or ambiguities (see Pitfalls).

Problem

Standard RF assumes independent observations. In longitudinal data each subject i contributes n_i rows (layout in tab. 1), so rows are clustered and may be serially correlated §1. Observation times can be equally or irregularly spaced, and missing visits are common and may be informative §1. Applying standard RF directly causes two problems: (a) observation-level bootstrap samples tend to contain rows from nearly every subject, so trees are correlated or near-identical, which hurts prediction §2.1 fig. 1; (b) out-of-bag (OOB) error is too optimistic, because OOB rows closely resemble in-bag rows from the same subject §2.1. Ignoring clustering in tree methods has also been reported to produce spurious subgroups and inaccurate variable selection §1 §2.1. Scope: CART-based RF only, categorical and quantitative outcomes; GUIDE and conditional inference trees are out of scope §1 §5.

Method

The review sorts the methods into three groups. Notation: subject i=1..N has observations j=1..n_i, with m responses and p predictors per row §1 tab. 1.

1. Clustered data, no time effect §2.2

Assumed model: subject mean plus independent noise eq. 1:

y_{ij} = \mu_i + \varepsilon_{ij} eq. 1

with \varepsilon_{ij} mean zero and independent within and across subjects §2.2.

• Averaging §2.2.1: collapse each subject's replicates to one row, then run standard RF. This loses within-subject variation, hides unbalanced n_i, and cannot average categorical values that differ between replicates. It was reported to degrade more than standard RF as N shrinks §2.2.1.
• Subject-level bootstrap (RF++) §2.2.2 fig. 2: for each tree, bootstrap subjects and include all of their rows as in-bag. Subject-level prediction is possible by majority vote over a subject's rows, which also gives a subject-level OOB misclassification rate §2.2.2.
• Two-stage bootstrap (Adler et al.) §2.2.2: resample subjects, then keep one randomly chosen row per selected subject. The cited simulations found this best against standard RF, averaging and RF++, with the caveat that results depend on the setting §2.2.2.
• Limitation: sampling only adjusts for clustering. There are no random effects, and time order is not used §2.2.2.

# Per-tree training sample for clustered data  [§2.2.2] [fig. 2]
def tree_sample(rows_by_subject, two_stage):
    subjects = list(rows_by_subject.keys())
    # Bootstrap at the subject level, not the row level  [§2.2.2]
    drawn = sample_with_replacement(subjects, size=len(subjects))
    in_bag = []
    for subject in drawn:
        rows = rows_by_subject[subject]
        if two_stage:
            # Adler et al.: one random row per drawn subject  [§2.2.2]
            in_bag.append(random_choice(rows))
        else:
            # RF++: every row of the drawn subject is in-bag  [§2.2.2]
            in_bag.extend(rows)
    # OOB = rows of subjects never drawn  [§2.2.2]
    out_of_bag = [
        row for subject in subjects if subject not in drawn
        for row in rows_by_subject[subject]
    ]
    return in_bag, out_of_bag

The OOB rule in the last step (whole subjects out of bag) follows from subject-level resampling. The review does not spell it out.

2. Univariate response with time effect §2.3

Historical RF §2.3.1. Data \{y_{ij}, t_{ij}, X_{ij}\}. The model for y_{ij} uses the concurrent (t_{ij}, X_{ij}) plus all earlier observations of subject i before t_{ij}. Time-invariant predictors are split as in standard CART. For a time-varying variable k, the history is first reduced by a summary function. An example is a count of past values, over both responses and predictors, inside a look-back window that fall below a threshold:

s_c(\eta;\bar{z}_{ijk}) = \sum_{t_{il}\in [t_{ij}-\eta_1,\, t_{ij})} I(z_{ilk} < \eta_2) eq. 2

where \bar z_{ij} = \{z_{il}=(y_{il},x_{il}) : t_{il}<t_{ij}\} §2.3.1. Alternatives are the relative frequency above a threshold, or the mean. Each has a windowed variant with t_{il}\in[t_{ij}-\eta_1, t_{ij}-\eta_3) §2.3.1. Only the frequency-based summaries are scale-invariant §2.3.1. A node splits on the predictor with the lowest Gini impurity (classification) or sum of squared errors (regression). For time-varying predictors the search also optimises \eta, which is costly when there are many such predictors §2.3.1. To reduce the cost, the htree authors optimise the cut-off on a random sub-sample of the windowed observations, and use subject-level resampling to build the forest §2.3.1. Variable importance is computed by deleting a variable: the error with predictor k removed minus the error with it included §2.3.1.

# Candidate split on a time-varying variable k  [§2.3.1] [eq. 2]
def count_summary(history, t_now, k, window, threshold):
    # history: past (time, z) pairs of one subject, z = (y, x)
    total = 0
    for time, z in history:
        # Look back at most `window` time units, excluding now  [eq. 2]
        if t_now - window <= time < t_now:
            if z[k] < threshold:
                total += 1
    return total

def best_split_time_varying(node_rows, k, eta_grid):
    best = None
    for window, threshold in eta_grid:
        # One summary value per row for this eta  [§2.3.1]
        values = [
            count_summary(r.history, r.time, k, window, threshold)
            for r in subsample(node_rows)  # extra randomisation [§2.3.1]
        ]
        # Standard CART cut-point search on the summary  [§2.3.1]
        candidate = best_cutpoint(values, node_rows)
        best = better_of(best, candidate)
    return best

The candidate set of \eta values (eta_grid) and the sub-sample size are not given in the review.

LMM-based hybrids, Gaussian outcome §2.3.2 tab. 2. The classic LMM is y_i = X_i\beta + Z_ib_i + \epsilon_i with b_i\sim N(0,D) and \epsilon_i\sim N(0,R_i), fitted by ML or REML §2.3.2. The hybrids replace X_i\beta by a tree or forest:

y_i = f(X_i) + Z_ib_i + \epsilon_i,\quad b_i\sim N(0,D),\ \epsilon_i\sim N(0,R_i) eq. 3

with R_i=\sigma^2 I_{n_i} assumed for MERT/MERF §2.3.2.

# MERT / MERF EM-style fitting  [§2.3.2]
def fit_merf(y, X, Z, subjects, use_forest):
    # Initialise random effects and variances  [§2.3.2]
    b = {i: zeros(q) for i in subjects}
    sigma2 = 1.0
    D = identity(q)
    while not converged(gll(f_hat, b, sigma2, D)):
        # Remove current random effects from the outcome  [§2.3.2]
        y_star = {i: y[i] - Z[i] @ b[i] for i in subjects}
        # Tree (MERT) or RF (MERF) on adjusted outcome  [§2.3.2]
        f_hat = fit_tree_or_rf(y_star, X, forest=use_forest)
        # Fit LMM y = f_hat(X) + Z b + eps for b, sigma2, D  [§2.3.2]
        b, sigma2, D = fit_random_effects(y, f_hat(X), Z)
    return f_hat, b, sigma2, D

The convergence criterion is the generalized log-likelihood §2.3.2:

GLL(f,b_i\mid y)=\sum_{i=1}^N\big\{[y_i-f(X_i)-Z_ib_i]^T R_i^{-1}[y_i-f(X_i)-Z_ib_i] + b_i^T D^{-1} b_i + \log|D| + \log|R_i|\big\} §2.3.2

The review states that MERF draws bootstrap samples at the observation level and uses OOB predictions to limit overfitting. Row-level resampling is justified because the random effects are assumed to absorb all within-subject correlation §2.3.2.

RE-EM tree §2.3.2. Step 1 is the same as MERT. In step 2 the tree only defines a partition: build the n_i\times T indicator matrix \Phi^i with \Phi^i_{jt}=I(y_{ij}\in g_t) over the T leaves g_t, then fit the LMM y_i=\Phi^i\mu + Z_ib_i+\epsilon_i, where \mu holds the per-leaf fixed effects §2.3.2. It uses nlme::lme, so R_i may be non-diagonal (for example autocorrelated) §2.3.2. REEMforest averages K RE-EM trees: \hat f(X_i)=\frac1K\sum_k \Phi^{i,k}\hat\mu_k §2.3.2.

Stochastic variants (SMERT, SMERF, SREEMtree, SREEMforest) §2.3.2 add a subject-specific centred Gaussian process \omega_i with \mathrm{Cov}(\omega_i(s),\omega_i(t))=\gamma^2\Gamma(s,t): y_i=f(X_i)+Z_ib_i+\omega_i+\epsilon_i, with b_i, \epsilon_i and \omega_i mutually independent. The adjusted outcome becomes y_i^\ast=y_i-Z_i\hat b_i-\hat\omega_i §2.3.2.

Prediction §2.3.2: for a new subject, use only the fixed part (\hat f(X) or the leaf effect \hat\mu_t). For a subject seen in training, add its estimated random effect.

GLMM-based hybrids §2.3.2. The GLMM is g(\mu_i)=\eta_i=X_i\beta+Z_ib_i with b_i\sim N(0,D) and \mu_i=E(y_i\mid b_i), fitted by PQL, IRLS or Newton-Raphson §2.3.2.

• GMERT: PQL in which a weighted MERT pseudo-model replaces the weighted LMM pseudo-model. It needs initial \hat\mu_i. Its simulation used \hat\mu_{ij}=0.25 if y_{ij}=0 and 0.75 if y_{ij}=1, with no further justification given §2.3.2.
• GMET / GMERF: a GLMM version of RE-EM, with leaf indicators as fixed effects and tree outcome y_i^\ast=\eta_i-Z_ib_i. Here \eta_i is first estimated with a standard GLM on all predictors, so the method fails when p greater than the number of observations. GMERF replaces the tree with an RF §2.3.2.

BiMM tree/forest (Bayesian, binary outcome) §2.3.3: \mathrm{logit}(\mu_{it})=\mathrm{CART}(X_{it})\beta+Z_{it}b_{it}, where \mathrm{CART}(X_{it}) are leaf-membership indicators §2.3.3.

# BiMM tree / forest  [§2.3.3]
def fit_bimm(y, X, Z, tol, use_forest):
    y_work = y
    previous = -inf
    while True:
        # Grow CART (or RF) on current outcome, get probabilities  [§2.3.3]
        p_tree = fit_tree_or_rf(y_work, X, forest=use_forest).predict_proba(X)
        # Bayesian GLMM on (y, p_tree, Z)  [§2.3.3]
        glmm = fit_bayesian_glmm(y, p_tree, Z)
        q = glmm.predicted_probability()
        # New outcome: discretise y + q  [§2.3.3]
        y_work = discretize(y + q)
        # Stop when posterior log-likelihood change < tol  [§2.3.3]
        if abs(glmm.posterior_loglik - previous) < tol:
            break
        previous = glmm.posterior_loglik
    return glmm

The discretisation rule for y+q is not stated in the review.

3. Multivariate responses §3

These methods treat a subject's repeated outcomes, or several outcomes, as one vector, and change the node impurity §3.1. They apply directly only to time-invariant predictors. Time-varying predictors need the techniques from Group 2 as well §3.

Distance-based impurity §3.1.1.

• Segal (1992), generalized sum of squares with a node covariance \mathbf V(\theta,t) eq. 4:

  SS(t)=\sum_{i\in t}(y_i-\bar y_t)'\mathbf V(\theta,t)^{-1}(y_i-\bar y_t) eq. 4

  Split gain \phi(s,t)=SS(t)-SS(t_L)-SS(t_R). To keep \phi\ge 0, the covariance parameters are estimated in the parent and shared by both children. Covariance choices: independence, AR(1), compound symmetry, or sample covariance. All subjects must share the same observation times §3.1.1.

• With independence, this reduces to the multivariate sum of squares \sum_{i\in t}\sum_{j}(y_{ij}-\bar y_j)^2, used in mvpart (De'Ath) and in Segal and Xiao's multivariate RF §3.1.1.

• Larsen and Speckman: Mahalanobis distance with one covariance \mathbf S estimated once on the full data and kept fixed during growth §3.1.1.

• db-MRT / Sim et al.: impurity from a distance matrix, Imp(t)=\sum_{i,j\in t}D_{ij}^2, so responses may have different lengths or irregular times as long as a distance is defined. Prediction within a leaf then needs care, because a plain average may not exist §3.1.1.

Likelihood-based impurity §3.1.2.

• Zhang (1998), multiple binary responses: exponential-family model f(y_i;\Psi,\theta)=\exp(\Psi'y_i+\theta w_i-A(\Psi,\theta)) with w_i=\sum_{j<k}y_{ij}y_{ik}. Node impurity is the maximised log-likelihood within the node §3.1.2. Ordinal responses are first turned into binary indicators (Zhang and Ye 2008) §3.1.2.
• Abdolell et al. (2002), multivariate normal: per-observation deviance D(\mu_i;y_i)=(y_i-\mu_i)'\Sigma^{-1}(y_i-\mu_i) with \Sigma known and constant. The node deviance sums this using the REML \hat\mu, and the impurity is its negation. This equals a Mahalanobis distance §3.1.2.

Multiple longitudinal responses §3.2: a subject's response becomes a matrix. Options: db-MRT with a suitable distance; combining Group 3 with subject-level bootstrapping or mixed effects; or Yu and Lambert (1999), who fit a spline to each trajectory and regress the basis coefficients with a multivariate tree. This also handles irregular times §3.2.

Hyperparameters

The review defines no training configuration. The only concrete values it reports are below.

Results

The review runs no experiments. All findings below come from the cited papers as summarised by the authors.

Software overview tab. 3 (T = tree, F = forest, R = regression, C = classification, P = permutation importance):

The authors note that MERT/MERF code is not publicly available, and that the packages have never been benchmarked against each other §4 §5. Except for BiMM and REEMforest, the methods have not been evaluated on high-dimensional data §5.

Paper vs code

No official code.

Pitfalls

• Use the original papers for implementation details. The review describes each algorithm "briefly" and says some are "slightly modified" §2.3.2. Number of trees, mtry, stopping rules, convergence tolerances and the \eta search are not given.
• db-MRT split sign. The review writes \phi(s)=-(Imp(t)-Imp(t_L)-Imp(t_R)) and says to take the maximum §3.1.1. Taken literally, this picks the split with the smallest impurity reduction. You almost certainly want to maximise Imp(t)-Imp(t_L)-Imp(t_R). Check the original. The impurity is also shown without normalisation by node size §3.1.1. Check whether the source divides by |t|, because unnormalised pairwise sums grow quadratically with node size.
• "l_1-norm" is really squared Euclidean. The identity linking db-MRT to Segal's independence case, \sum|y_i-\bar y|^2=\frac1{2N}\sum_i\sum_j|y_i-y_j|^2, holds for squared Euclidean (l_2) distances, not l_1 as the text says §3.1.1.
• Segal independence formula. It is written with \bar y_j (a global per-component mean) §3.1.1. Inside a node, the mean must be the node mean \bar y_{t,j}, consistent with eq. 4.
• Same covariance for parent and children. In Segal's criterion, re-estimating \mathbf V in each child can make \phi negative. Estimate it in the parent and reuse it §3.1.1.
• RE-EM indicator. \Phi^i_{jt}=I(y_{ij}\in g_t) §2.3.2 means "observation j (its predictors x_{ij}) falls in leaf t", not the response value.
• Two-stage bootstrap wording. The review first says all of a chosen subject's rows are in-bag, then that one row per subject is kept §2.2.2. The intended scheme, as compared in the simulations, is one row per resampled subject.
• OOB leakage. With observation-level bootstrapping, OOB error on clustered data is optimistic §2.1. If you use subject-level resampling, compute OOB per subject. MERF deliberately keeps observation-level resampling, which is valid only if the random effects fully capture the within-subject correlation §2.3.2.
• New vs known subjects. Mixed-effects hybrids predict differently for unseen subjects (fixed part only) and for training subjects (fixed part plus \hat b_i) §2.3.2. Do not add random effects you do not have.
• High-dimensional GLMM hybrids. GMET/GMERF initialise \eta_i with a plain GLM on all predictors, which is impossible when p exceeds the number of observations §2.3.2 §5.
• Historical RF cost. Optimising \eta inside every split for every time-varying predictor is expensive with omics-scale p §2.3.1. Delete-variable importance needs refitting once per variable, and it can hide important variables that are correlated with others §2.3.1.
• Scale invariance. Only the count and frequency summaries keep RF's invariance to monotone transforms. Mean summaries do not §2.3.1.
• Typo in tab. 1. The last row lists x_{Nn_11}\dots x_{Nn_1p}; it should be x_{Nn_N1}\dots x_{Nn_Np}.
• Missing data are named as a major challenge, but only Segal's surrogate splits are mentioned as a remedy §5.

Lineage

• Builds on: Breiman's random forest and CART §2.1.
• Reviewed methods: RF++ (Karpievitch et al. 2009), Adler et al. 2011a/b, htree (Sexton and Laake), MERT/MERF (Hajjem et al. 2011, 2014), RE-EM tree (Sela and Simonoff 2012), REEMforest and stochastic variants (Capitaine et al. 2021), GMERT (Hajjem et al. 2017), GMET/GMERF (Fontana et al. 2018; Pellagatti et al. 2021), BiMM (Speiser et al. 2019, 2020), Segal 1992, De'Ath 2002, Larsen and Speckman 2004, Sim et al. 2013, Zhang 1998, Abdolell et al. 2002, Yu and Lambert 1999 §2 §3.
• Out of scope: GUIDE and conditional-inference tree forests §5.