Individual participant data meta-analysis with mixed-effects transformation models

Summary One-stage meta-analysis of individual participant data (IPD) poses several statistical and computational challenges. For time-to-event outcomes, the approach requires the estimation of complicated nonlinear mixed-effects models that are flexible enough to realistically capture the most important characteristics of the IPD. We present a model class that incorporates general normally distributed random effects into linear transformation models. We discuss extensions to model between-study heterogeneity in baseline risks and covariate effects and also relax the assumption of proportional hazards. Within the proposed framework, data with arbitrary random censoring patterns can be handled. The accompanying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\textsf{R}$\end{document} package tramME utilizes the Laplace approximation and automatic differentiation to perform efficient maximum likelihood estimation and inference in mixed-effects transformation models. We compare several variants of our model to predict the survival of patients with chronic obstructive pulmonary disease using a large data set of prognostic studies. Finally, a simulation study is presented that verifies the correctness of the implementation and highlights its efficiency compared to an alternative approach.

. In the application presented in this text, the results are very similar from models with polynomial order above four, so we used this order for all six model variants.
The integral in the log-likelihood function (see Equation 3.6 of the article) is evaluated using the Laplace approximation. The Laplace approximation is based on the second order Taylor expansion of the joint log-likelihood function at the conditional mode of the random effects, b i (ϑ, Σ) = arg max bi L i (ϑ, Σ, b i ). In the context of IPD meta-analysis, these random effects predictions are the deviations of the individual studies from the mean effect sizes. The approximation of the log-likelihood function is .
In order to improve numerical properties of the models estimated in example application, we rescaled the covariates using their ranges. This way, convergence for the otherwise complicated models was more easily attained, and the results became more stable. Wherever it was necessary, the results were cast back to their original scale. The corresponding standard errors were recalculated using the delta method.
The maximum likelihood estimation of mixed-effects transformation models is implemented in the R add-on package tramME (Tamási and Hothorn, 2021), and the fixed-effects transformation models are estimated using the tram and mlt packages (Hothorn, 2020). For the implementation details and several worked example applications, see the package vignettes.
The model specifications described in Table 1 and used in our example application in Section 4 can be estimated with the following commands.
R> library("tram") R> library("tramME") R> sup <-c(1, 150) ## setting the support for basis approximation R> opt <-optim_control(iter.max = 1e4, eval.max = 1e4, rel.tol = 1e-8) The estimation on a subset of the original 3CIA sample can be performed by running the demo "IPD-MA" of the tramME package: R> demo("IPD-MA", package = "tramME") A.2 Comparison of tramME and merlin To further validate the proposed models and implementation, we have also replicated the core results in Stata using the mixed-effects model developed by Crowther and others (2014), and implemented in the merlin package (Crowther, 2020). merlin fits a wide range of multivariate multilevel regressions, of which the mixed-effects Royston-Parmar model is a special case. Similarly to the transformation models proposed in the main text, the model by Royston and Parmar (2002) parameterizes the log-cumulative hazards.
Although the two models are formally very similar, there are several important computational differences between the two implementations. Table S-1 gives a summary of these. merlin uses the original Royston-Parmar formulation and utilizes natural cubic splines to approximate the baseline log-cumulative hazards without introducing explicit constraints to ensure the monotonicity of the function (see Royston andParmar, 2002, pp. 2177). merlin estimates the fixed effects only, univariate models with analytic gradients, and hence fitting Models 5 and 6 are computationally highly efficient. However, models with random effects firstly use adaptive Gauss-Hermite quadrature to calculate the likelihood, and secondly use finite differences for all gradient calculations. In contrast, tramME uses polynomials in Bernstein form with corresponding explicit constraints on their parameters to capture the baseline log-cumulative hazard and the time-varying effects.
Building on the automatic differentiation and Laplace approximation provided by the TMB package (Kristensen and others, 2016), the estimation of the presented models with tramME is very efficient. The slower computations render the fitting of the mixed-effects models with merlin to the 3CIA dataset practically infeasible.
The Stata code for estimating the six model variants with the merlin package is available from Github at https://github.com/btamasi/tramme-ipd-ma.
We compare the numerical results from the fitting Model 5 (homogeneous, time-invariant effects of the predictors, stratified baseline risk) with both tramME and merlin to verify the

B.1 Censoring
The outcome variable, the survival time of the participants, is recorded with different precision in different cohorts. In most of the cases, the follow-up time is rounded to months, while in other cohorts the same information is available in fractions of months. It is clear, that having only rounded follow-up times represents partial information on the actual times of events, and a correct statistical approach must take this issue into account. For this reason, we treat the observations recorded with follow-up times rounded to zero and one decimal places as interval-censored observations, while survival times with two decimal places as exact observations. The likelihood contributions of exact continuous and interval-censored observations are given by Equation (3.4) and by the third case in Equation (3.5) of the article, respectively. Our proposed estimation procedure can easily accommodate both types of observations.
The distribution of censoring times for each cohort is given by Figure     stratified baseline risks, Table 1). In this case, the mean and standard deviation of the linear predictor describe the location and the variability of a weighted sum of the prognostic factors in each cohort. As Table S-3 shows, there are differences in both the location of the linear predictor and its variability across studies.
To gauge the level of heterogeneity in the prognostic factor effects as well as the baseline risks, we first consider the "fully stratified model", i.e., fitting a Cox proportional hazards model with a fully parametric baseline cumulative hazard to each cohort separately. This approach corresponds to the first step of a two-stage IPD meta-analysis using a specification analogous to Model 2 (see Table 1). Fitting separate models to each cohort is considered suboptimal, because the estimates for the different studies do not borrow information from others, but it is a natural first step in assessing the heterogeneity in the prognostic factor effects, as it does not introduce too many additional assumptions. For simplicity, we only estimate the model with time-constant effects.  which are interpreted as log-hazard ratios in our model specifications.  Table 1) is violated, as some of the baseline transformation functions cross.

C. Additional results from model estimation
This section presents some additional results from the IPD meta-analysis of the 3CIA dataset using transformation models.   cohort in the 3CIA dataset from the mixed-effects model specifications (Models 1 to 4, in Table 1).
Figure S-7 presents the baseline transformation functions, (log-cumulative baseline hazards) estimated using Models 1 to 6 (see Table 1) for each study in the 3CIA dataset. Note that Models 1 and 3 assume proportional baseline hazards, which is reflected in the fitted curves in this plot.

D. Additional simulation results
In this section, we discuss further results of the simulation study that compares the mixedeffects transformation model approach for IPD-MA with the results obtained from the procedure outlined in Garcia and others (2019).
The estimation of the variance components parameters can be evaluated based on Figure S-8.
We generated a four-dimensional random-effects vector for each study in the dataset (random frailty term and three random slopes). Because the procedure by Garcia and others (2019) consists of re-estimating binary mixed-effects models at each time point, at which we want to evaluate the distribution of the outcome, we obtain five sets of random-effects covariance parameters, which we present separately. Similarly to our previous results, both of the approaches fare well in the specific estimation task, with the mixed-effects transformation model having slightly smaller variances. Slight downward bias can be observed in the variance estimates of the random effects with both approaches, which is in line with previous findings about maximum likelihood estimates of variance parameters in a mixed-effects IPD meta-analysis setting (Crowther and others, 2012).
We also compare the models by their ability to predict the true values of the random effects. As To evaluate the effect of misspecifying the inverse link function in a mixed-effects transformation model, we run a separate simulation experiment, which uses the data generating mechanism that corresponds to Model 1, with parameter values estimated from the 3CIA dataset. Instead of fitting the correct model, we fit a specification that uses expit = logit −1 , the CDF of the standard logistic regression, as the inverse link function: This specification corresponds to a proportional odds (PO) models. To compare the results with the correctly specified model, as well as the true effects, we transform the ef- fect estimates from both models to the probability scale and calculate probabilistic indices, Thas and others, 2012). The probabilistic index (PI) gives the probability that an individual with a vector of covariate values (x) has shorter survival time than an other individual in the same study that has the same predictor values except for one, which is larger by one unit (x ). By integrating over the distribution of the corresponding random effects, we can calculate the expected PI that will serve as a common scale to compare our