Highlights

2 Low back pain data

Patel et al.Reference Patel, Hee and Mistry ¹ collected data from 19 completed NSLBP trials to identify participant characteristics predicting clinical response to treatments for low back pain. Following a data-sharing agreements, we obtained data from 8 out of the 19 eligible NSLBP trials. As identifying subgroups with differential treatment effects requires data from two-arm trials, with consistent measurement of outcome variables and subgroup-defining factors across all trials, the number of trials available for the subgroup identification is further reduced. We focused on two-arm trials that included the Roland Morris Disability Questionnaire (RMDQ) measured at baseline, denoted as RMDQ_0, and at least one follow-up visit. Therefore, from the eight eligible NSLBP trials, we selected four trials for our analysis including 1780 individuals. To enable data publishing and reproducibility, we created a synthetic data set of these four trials using the R-package synthpop with the option method="parametric". The synthesis of the data results in randomization ratios that are comparable to those in the underlying real trials, as well as similar distributions for the continuous covariates. However, these ratios and distributions are not perfectly matched. In particular, the randomization ratios may yield values that appear unusual. Nevertheless, this discrepancy does not impact the main objectives of illustrating different subgroup identification approaches. The synthetic data and the R code used for the analyses in the subsequent sections can be found in the Supplementary Material.

Table 1 presents demographics and the RMDQ stratified by trial. The RMDQ score ranges from 0 to 24 with higher values indicating a poorer outcome associated with back pain. The table illustrates that the trials are of different sizes. The largest trial, trial with ID 1, includes more than $1000$ patients, while trial 3 involves only $53$ patients. RMDQ was measured at different follow-up visits. Our analyses are based on the data summarized in Table 1. The synthetic data include four studies in which RMDQ is measured at baseline (RMDQ_0) and at least at one follow-up visit. As outcome, we considered the last observation of RMDQ of each patient.

Table 1 Summary of the synthesized analysis set

Note: Categorical variables are presented as n (%), while continuous variables are reported as median [25%, 75% quantile]. Each trial is depicted in a separate column.

Table 1 illustrates that the baseline characteristics, sex, age, and RMDQ at baseline slightly differ between trials. We investigate the treatment effect in the pooled data. Using trial as fixed effect in the linear model as shown in the row “Linear Model Adjusted for Trial,” shows a significant reduction of the RMDQ in the intervention compared to the control arm (see Table 2) indicating a benefit of the intervention over the control treatment. The RMDQ scores at last follow-up differ across trials as indicated by the coefficients of the trial indicators. To account for heterogeneity in the treatment effect using random effects, we applied a linear mixed model to the entire pooled population, incorporating the treatment indicator as random effects and the trial indicator as fixed effect. The variance of the random treatment effect was estimated to be larger than zero, see Table 2. By accounting for between-trial variability in the treatment effect, this model introduces greater uncertainty in the pooled estimate of the treatment, resulting in a wider confidence interval compared to the model presented in the first row of Table 2.

Table 2 Treatment effect estimated by an unadjusted linear model, a linear model including the trial indicator as fixed effect and a linear mixed model with treatment as random and trial indicator as fixed effect

3 Subgroup identification

We illustrate the identification of subgroups in NSLBP with differential treatment effects using MOB approaches. MOB and its extensions, as metaMOB or palmtree,Reference Seibold, Hothorn and Zeileis ¹⁵ are based on generalized linear models. To describe the different variations of MOB and the underlying models, we assume that the outcome RMDQ is denoted by y. Furthermore, the treatment indicator is denoted by t, the covariates age, sex, and RMDQ at baseline are denoted by $x_{age}$ , $x_{sex}$ , and $x_{RMDQ\_0}$ . The baseline covariates $x_{age}$ , $x_{sex}$ , and $x_{RMDQ\_0}$ are considered as potential splitting variables and are therefore not involved in the underlying regression model. The analysis includes four trials $k=1,\ldots , 4$ . The corresponding trial indicator is denoted by $x_{trial}$ .

For MOB, the outcome of each subgroup j and trial k is modeled by

(1) $$ \begin{align} y_{jk}=\gamma_{j}+\beta_j t. \end{align} $$

MOB does not account for the data being pooled from different trials; therefore, the right-hand side of the equation does not depend on k. Adjusting for the different trials by fixed effects using MOB is referred to as MOB-SI. SI refers to a stratified intercept as for each trial a separate intercept is estimated. The linear model for RMDQ fitted in each subgroup j based on the method MOB-SI is

(2) $$ \begin{align} y_{jk}=\gamma_{jk}+\beta_j t, \end{align} $$

with $\gamma _{jk}$ describing the subgroup and trial-specific fixed intercept. Both MOB and MOB-SI can be fitted using the lmtree function of the partykit package. Addressing heterogeneity in the baseline with random intercepts can be achieved by applying MOB-RI. MOB-RI is based on GLMM-trees. Therefore, for analyzing the NSLBP data with RMDQ as outcome, a linear mixed model is fitted to each subgroup j and trial k:

(3) $$ \begin{align} y_{jk}= \gamma_{j}+\beta_jt_k+b_{0k} \text{ with } b_{0k}\sim \mathcal{N}(0,\tau_0^2). \end{align} $$

The random intercept $b_{0k}$ is considered to be the same for each subgroup. Accounting additionally for heterogeneity in the treatment effect is feasible using the following extended metaMOB approaches for IPD meta-analyses, namely metaMOB-RI and metaMOB-SI.

metaMOB-RI:

(4) $$ \begin{align} y_{jk}= \gamma_j+\beta_jt_k+b_{0k}+b_{1k}{t_k} \text{ with } b_{0k}\sim \mathcal{N}(0,\tau_0^2)\text{ and } b_{1k}\sim \mathcal{N}(0,\tau_1^2), \end{align} $$

and metaMOB-SI:

(5) $$ \begin{align} y_{jk}= \gamma_{jk}+\beta_jt_k+b_{1k} {t_k}\text{ with } b_{1k}\sim \mathcal{N}(0,\tau_1^2). \end{align} $$

The approaches involving random effects, MOB-RI, metaMOB-SI, and metaMOB-RI are fitted using the lmertree function of the glmertree package.

An alternative to both the MOB-RI and MOB-SI approaches, which address heterogeneity in baseline by using subgroup and trial-specific intercepts (MOB-SI) or random intercepts for the trial indicator (MOB-RI), is the Generalized Linear Model Trees with global additive effects, referred to as palmtree.Reference Seibold, Hothorn and Zeileis ¹⁵ In contrast to MOB-RI which assumes the random intercepts to be constant across subgroups, palmtree includes the treatment indicator in the model similar to MOB-SI, but assumes these intercepts to be the same across the identified subgroups:

(6) $$ \begin{align} y_{jk}=\gamma_k+\beta_j t. \end{align} $$

The regression models used for the metaMOB approach (Equations (4) and (5)) are in alignment with the models typically used in random-effects meta-analysis, here the normal–normal hierarchical model.

3.1 MOB

For the NSLBP data, approaches that account for between-trial heterogeneity are more suitable due to various differences between the trials included. Therefore, MOB using Equation (1) is not employed for data analysis. MOB-SI accounts for heterogeneity in the baseline by adjusting for the trial indicator, see Equation (2). The result of this approach is illustrated in Figure 1. The subgroups are defined by the RMDQ_0 and Age. MOB-SI partitions the group of participants with RMDQ_0 values larger than 9 into three subgroups by additional splits on RMDQ_0 and Age. MOB-SI estimates the treatment effect separately in each of the identified subgroups using a linear model with treatment and trial as factors (see Equation (2)), resulting in p-values for the treatment effect within two subgroups, denoted node 4 and node 6, that are smaller than $0.05$ . Both of these subgroups estimate a reduction of the RMDQ score of the intervention compared to the control indicating a treatment benefit. In node 6, which includes only subjects with baseline RMDQ values greater than 6 and less than or equal to 9, the estimated benefit of the intervention is the largest, with a predicted reduction of $-$ 1.66 points in RMDQ in the intervention group compared to the control group. The estimated treatment effect is based on Equation (2) that is also used for the MOB-SI approach.

Figure 1 Tree obtained by MOB-SI.

Note: Five subgroups are identified. The upper boxplots display the outcome values stratified by the treatment indicator, while the lower boxplots show the outcome values stratified by the trial indicator. Note that the RMDQ, illustrated on the y-axis, ranges from 0 to 24, with higher values indicating a poorer outcome. To the best of our knowledge, the y-axis limits are hard-coded and cannot be manually adjusted.

Figure 2 Tree obtained by metaMOB-SI.

Note: Four subgroups are defined. All splits are performed on the variable RMDQ_0. The upper boxplots display the outcome values stratified by the treatment indicator, while the lower boxplots show the outcome values stratified by the trial indicator. Note that the RMDQ, illustrated on the y-axis, ranges from 0 to 24, with higher values indicating a poorer outcome. To the best of our knowledge, the y-axis limits are hard-coded and cannot be manually adjusted.

An alternative to MOB with stratified intercept for the trial is the GLMM-tree algorithm with a random intercept to account for baseline heterogeneity. This approach is referred to by MOB-RI,Reference Huber, Benda and Friede ⁶ which imposes a restriction on the heterogeneity of the baseline by assuming a certain distribution of its random effects, see Equation (3). Nevertheless, the subgroups identified by MOB-RI are the same as for MOB-SI. However, the estimated effects of MOB-RI differ slightly from MOB-SI as the model used for the estimation of the treatment effect is based on Equation (3).