# Joint Modelling Approaches to Survival Analysis via Likelihood-Based Boosting Techniques

Authors: Colin Griesbach; Andreas Groll; Elisabeth Bergherr
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

Joint models are a powerful class of statistical models which apply to any data where event times are recorded alongside a longitudinal outcome by connecting longitudinal and time-to-event data within a joint likelihood allowing for quantification of the association between the two outcomes without possible bias. In order to make joint models feasible for regularization and variable selection, a statistical boosting algorithm has been proposed, which fits joint models using component-wise gradient boosting techniques. However, these methods have well-known limitations, i.e., they provide no balanced updating procedure for random effects in longitudinal analysis and tend to return biased effect estimation for time-dependent covariates in survival analysis. In this manuscript, we adapt likelihood-based boosting techniques to the framework of joint models and propose a novel algorithm in order to improve inference where gradient boosting has said limitations. The algorithm represents a novel boosting approach allowing for time-dependent covariates in survival analysis and in addition offers variable selection for joint models, which is evaluated via simulations and real world application modelling CD4 cell counts of patients infected with human immunodeficiency virus (HIV). Overall, the method stands out with respect to variable selection properties and represents an accessible way to boosting for time-dependent covariates in survival analysis, which lays a foundation for all kinds of possible extensions.


---

## Body

## 1. Introduction

First suggested by [1], joint models were established as a valuable tool for analysing data where event times are measured alongside a longitudinal outcome. One naive approach of evaluating such frequently occurring data structures would be separate modelling, i.e., fitting appropriate models independently for longitudinal and time-to-event data. However, separate modelling neither corrects for event-dependent dropout in longitudinal analysis nor quantifies the relation between a time-dependent covariate and the risk for an event in survival analysis [2]. Various approaches have been proposed to address these issues, one being the Andersen-Gill model [3] for time-varying covariates in survival analysis or two-stage approaches, where longitudinal model fits are included as fixed covariates in time-to-event regression. It has been shown, however, that these methods tend to produce biased results [4, 5]. One solution therefore is combining both the survival and longitudinal models within one single joint likelihood. A wide introduction to this joint modelling framework is presented in [4] including the JM package [6]. Moreover, an evolution of joint model progression up to the year 2004 is provided in [7], and in addition, several Bayesian approaches have been carried out [8–10].Current joint modelling estimation methods, however, lack clear concepts for proper variable selection and good performance regarding prediction. Moreover they are not feasible for high-dimensional data, in particular where the number of covariates exceeds the number of observations, i.e.,p>n problems. In order to overcome these hindrances, an algorithm was initially proposed, where joint models are fitted with gradient boosting techniques, which are known for addressing exactly these issues [11]. Evolved from machine learning as an approach to classification problems originally proposed in [12], gradient boosting deals with high-dimensional data and the component-wise updating scheme offers implicit variable selection. The boosting algorithm for joint models was extended [13], but when wanting to lay more focus on the survival side of the model, gradient boosting proved to struggle with time-varying covariates in time-to-event analysis. This has also been observed for pure survival models in [14].Hence, this work focuses on likelihood-based boosting. First introduced in [15], likelihood-based boosting is designed to directly maximize a given likelihood where concepts, which gradient boosting implicitly offers, are reproduced artificially for regular optimization methods like Newton algorithms or Fisher scoring. The method was further developed for flexible semiparametric mixed models [16] and for several classes of generalised mixed models [17–20]. The R package GMMBoost [21] covers most of these approaches. Likelihood-based boosting has also been proved useful for survival analysis with time-varying effects [22], and a general overview is given in [23]. Since the random structure plays an important role as a connector between longitudinal and time-to-event data, we additionally incorporate a novel correction step within the estimation procedure for the random effects, which was first suggested in [24, 25] and reduces possible bias arising from wrongly identified random effects.The contribution of this work is the novel lbbJM boosting algorithm for joint models, which offers the first boosting-based regularization approach for time-dependent covariates in survival analysis and in addition new variable selection mechanics for joint models with focus on time-to-event analysis.The remainder of this manuscript is structured as follows: Section2 highlights the overall concepts of both joint modelling and likelihood-based boosting to give a sufficient understanding of the methods used in the following parts. Section 3 then contains a detailed description of the considered joint model together with the proposed boosting algorithm and its computational details. Sections 4 and 5 deal with applying the algorithm to different setups of simulated data as well as to the AIDS dataset [26] included in the JM package. Results and possible extensions are discussed in the final section.


## 2. Backgrounds

Before the algorithm is presented and discussed in detail, we briefly highlight the concepts of both joint modelling and likelihood-based boosting.


### 2.1. Joint Models

In general, a joint model consists of two parts, one longitudinal and one survival submodel. A popular view on joint modelling is to choose one model as themain model, whereas the other model then features the analysis occurring in the main model. With the primary outcome being longitudinal data, a survival model can be used to correct for event-dependent dropout in longitudinal analysis. For time-to-event data as outcome of interest, additional longitudinal modelling reduces measurement error on the one hand and, on the other hand, extrapolates only on single time points observed longitudinal data to continuous functions which are then included in survival analysis. We will from now on focus on joint models with time-to-event data as primary outcome.The longitudinal submodel takes the form:(1)y=ηlt,x+ε,where longitudinal outcome y is described by the longitudinal predictor function ηl depending on time t and a set of covariates x. Although t can be included in x, we will highlight it in the context of joint models, as the role of t is of greater importance. In the survival submodel, the hazard
(2)λt∣x=λ0texpηsx+αηlt,x,is modelled by a baseline hazard λ0t with multiplicative effects described by the survival predictor function ηs. In addition, the longitudinal predictor ηl is reappearing in the survival model, this time scaled by a factor α. The parameter α thus quantifies the association of the two submodels and is therefore called the association parameter. It can be interpreted as the impact a time-varying longitudinal covariate has on the hazard for an event.Parameter estimation for such joint models can be done in various ways. Two common methods are two-stage and joint likelihood approaches, respectively. In the former, the longitudinal model is estimated with the estimation method of choice leading to the model fitη^l, which is then carried forward as fixed covariate into survival analysis. In the latter, longitudinal and survival submodels are combined in a single joint likelihood. Let i=1,⋯,n denote clusters and j=1,⋯ni the repeated measurements. Assuming independent data generating processes for both submodels, the joint likelihood can then be written as
(3)Ly,T,δ=∏i=1n∏j=1niflyij∣ηlfsTi,δi∣α,ηl,ηs,λ0,with densities fl and fs for the longitudinal and survival submodels and time-to-event outcome T,δ=Ti,δii∈ℕ. Regular inference is now done by maximizing (3) using appropriate maximization methods, the most prominent one being EM algorithms.


### 2.2. Likelihood-Based Boosting

We intend to give a short description of the underlying mechanics used in the following section. The overall concept of likelihood-based boosting is to create an iterative and component-wise updating scheme, which eventually converges to a maximum likelihood estimator but is stopped early in order to prevent overfitting. Letβ model be the effect of p covariates. Likelihood-based boosting maximizes a given log-likelihood lβ by component-wise Fisher scoring in the following way:For each covariater∈1,⋯,p consider the subvector βr containing only the coefficients referring to the rth covariate. We compute the score vector and Fisher matrix as
(4)srβr=∂lβ∂βr,Frβr=−E∂2lβ∂βr∂βrT,and obtain a possible update
(5)ur≔Frβr−1srβr,for the rth component. Now, we determine the best performing covariate with respect to likelihood maximization, i.e., find the component
(6)∗=argmaxr=1,⋯,plβ~r,β~r=β1,⋯,βr+ur,⋯,βpwhere the corresponding update yields the biggest improvement of the likelihood. One receives a new model fit by weakly updating this best performing component, i.e., by scaling with a factor ν, the so called step length:
(7)βrnew=βrifr≠∗,βr+νurifr=∗,0<ν≤1,r=1,⋯,p.The step lengthν is controlling the weakness of the update to prevent overfitting and give every covariate a chance for selection. A popular choice in the literature is ν=0.1.Repeating this updating process for a sufficiently large number of iterations leads to the regular maximum likelihood estimator forβ. But instead the algorithm is stopped early to gain better prediction performance and variable selection. The optimal amount of iterations actually is a tuning parameter of the method and can be determined via cross-validation or by focusing on information criteria like AIC or BIC [27, 28].


## 2.1. Joint Models

In general, a joint model consists of two parts, one longitudinal and one survival submodel. A popular view on joint modelling is to choose one model as themain model, whereas the other model then features the analysis occurring in the main model. With the primary outcome being longitudinal data, a survival model can be used to correct for event-dependent dropout in longitudinal analysis. For time-to-event data as outcome of interest, additional longitudinal modelling reduces measurement error on the one hand and, on the other hand, extrapolates only on single time points observed longitudinal data to continuous functions which are then included in survival analysis. We will from now on focus on joint models with time-to-event data as primary outcome.The longitudinal submodel takes the form:(1)y=ηlt,x+ε,where longitudinal outcome y is described by the longitudinal predictor function ηl depending on time t and a set of covariates x. Although t can be included in x, we will highlight it in the context of joint models, as the role of t is of greater importance. In the survival submodel, the hazard
(2)λt∣x=λ0texpηsx+αηlt,x,is modelled by a baseline hazard λ0t with multiplicative effects described by the survival predictor function ηs. In addition, the longitudinal predictor ηl is reappearing in the survival model, this time scaled by a factor α. The parameter α thus quantifies the association of the two submodels and is therefore called the association parameter. It can be interpreted as the impact a time-varying longitudinal covariate has on the hazard for an event.Parameter estimation for such joint models can be done in various ways. Two common methods are two-stage and joint likelihood approaches, respectively. In the former, the longitudinal model is estimated with the estimation method of choice leading to the model fitη^l, which is then carried forward as fixed covariate into survival analysis. In the latter, longitudinal and survival submodels are combined in a single joint likelihood. Let i=1,⋯,n denote clusters and j=1,⋯ni the repeated measurements. Assuming independent data generating processes for both submodels, the joint likelihood can then be written as
(3)Ly,T,δ=∏i=1n∏j=1niflyij∣ηlfsTi,δi∣α,ηl,ηs,λ0,with densities fl and fs for the longitudinal and survival submodels and time-to-event outcome T,δ=Ti,δii∈ℕ. Regular inference is now done by maximizing (3) using appropriate maximization methods, the most prominent one being EM algorithms.


## 2.2. Likelihood-Based Boosting

We intend to give a short description of the underlying mechanics used in the following section. The overall concept of likelihood-based boosting is to create an iterative and component-wise updating scheme, which eventually converges to a maximum likelihood estimator but is stopped early in order to prevent overfitting. Letβ model be the effect of p covariates. Likelihood-based boosting maximizes a given log-likelihood lβ by component-wise Fisher scoring in the following way:For each covariater∈1,⋯,p consider the subvector βr containing only the coefficients referring to the rth covariate. We compute the score vector and Fisher matrix as
(4)srβr=∂lβ∂βr,Frβr=−E∂2lβ∂βr∂βrT,and obtain a possible update
(5)ur≔Frβr−1srβr,for the rth component. Now, we determine the best performing covariate with respect to likelihood maximization, i.e., find the component
(6)∗=argmaxr=1,⋯,plβ~r,β~r=β1,⋯,βr+ur,⋯,βpwhere the corresponding update yields the biggest improvement of the likelihood. One receives a new model fit by weakly updating this best performing component, i.e., by scaling with a factor ν, the so called step length:
(7)βrnew=βrifr≠∗,βr+νurifr=∗,0<ν≤1,r=1,⋯,p.The step lengthν is controlling the weakness of the update to prevent overfitting and give every covariate a chance for selection. A popular choice in the literature is ν=0.1.Repeating this updating process for a sufficiently large number of iterations leads to the regular maximum likelihood estimator forβ. But instead the algorithm is stopped early to gain better prediction performance and variable selection. The optimal amount of iterations actually is a tuning parameter of the method and can be determined via cross-validation or by focusing on information criteria like AIC or BIC [27, 28].


## 3. Boosting Joint Models

### 3.1. The Model

Before we introduce the boosting algorithm, we describe the specific joint model. Withi=1,⋯,n denoting the individual and j=1,⋯,ni a single specific measurement, the longitudinal submodel is given by
(8)yij=ηltij,xli+εij=β0+βttij+βlTxli+γ0i+γtitij+εij,where yij is modelled by ηl depending on specific measurement times tij and longitudinal time-independent covariates xli∈ℝpl. This represents a standard linear mixed model with intercept β0 and fixed linear effects βt and βl of time and baseline covariates as well as individual specific random effects γ0i and γti with γ0i,γti~N⊗20,Q. The error terms εij are assumed to follow a normal distribution with Eεij=0 and Varεij=σ2>0.Please note that the model can be additionally extended to interaction effects of timetij with baseline covariates xti∈ℝpt by including the term βtTxtitij in (8). This results in slightly adjusted integrals in the survival part and is omitted in the following for the sake of better readability.In the survival part, the individual hazard(9)λit=λ0texpηsxsi+αηlt,xli,is modelled with the survival predictor ηsxsi=βsTxsi containing additional linear effects βs of baseline covariates xsi∈ℝps. To execute a full likelihood approach, the baseline hazard
(10)λ0t=∑k=1Kλk1Ikt,is chosen to be piecewise-constant depending on the number of segments K and their exact locations Ik=tk−1,tk with t0=0 and tk=maxT for k=1,⋯,K. The collection of values for the baseline hazard is denoted in λ=λkk=1,⋯,K. Later, we will choose K between 7 and 10 in order to guarantee substantially more flexibility than a constant baseline hazard without becoming computationally too demanding.Given two formulas (8) and (9), we can now calculate the joint log-likelihood. Let y=yiji=1,⋯,n,j=1,⋯,ni denote the collection of all longitudinal measurements. Assuming the time-to-event process is conditionally independent from the longitudinal random structure, the joint likelihood can be decomposed into a longitudinal and a survival term. Set γT=γii=1,⋯,n with γiT=γ0i,γti and ϑl≔β0,βlT,βt,γTT. Furthermore, τ contains information on variance-covariance components σ2 and Q. The unpenalized longitudinal log-likelihood is
(11)ℓlϑl,σ2∣y=∑i=1n∑j=1nilogϕyij∣ηltij,xli,σ2,where ϕ·∣m,v denotes the density of a normal distribution with mean m and variance v. Laplace approximation follows [29] and then leads to an additional quadratic penalty term for the random effects yielding the penalized log-likelihood:
(12)ℓlpenϑl,τ∣y=∑i=1n∑j=1nilogϕyij∣ηltij,xli,σ2−12γiTQγi.Note that for the penalized log-likelihoodτ substitutes σ2 as an argument, since the penalized log-likelihood additionally contains information of the variance matrix Q.On the other hand, for given survival outcomeT,δ=Ti,δii=1,⋯,n with event times Ti and censoring indicator δi, the survival log-likelihood takes the well-known form:
(13)ℓsϑs∣T,δ=∑i=1nδilogλiTi∣ηl,ηs,α,λ−∫0Tiexpλit∣ηl,ηs,α,λdt,with ϑs≔λT,α,βsTT. For ϑ≔ϑlT,ϑsTT, we finally receive the penalized and unpenalized joint log-likelihood:
(14)ℓϑ,σ2∣y,T,δ=ℓlϑl,σ2∣y+ℓsϑs∣T,δ,ℓpenϑ,τ∣y,T,δ=ℓlpenϑl,τ∣y+ℓsϑs∣T,δ.


### 3.2. The Algorithm

The following lbbJM algorithm (likelihood-based boosting for joint models) describes a way to fit the formulated joint model by likelihood-based boosting methods discussed in Section2.2.Algorithm 1: lbbJM.
(i)Initialize starting values β^00, β^t0, β^l0, λ^0, β^s0, α^0, γ^00, and γ^t0 with variance-covariance-components σ^20,Q^0 according to Section 3.3. Choose mmax,l,mmax,s and execute the following for every possible tuple ml,ms with ml<mmax,l, ms<mmax,s:Longitudinal part(ii)form=1 to mldostep1: Update fixed effectsForr=1,⋯,pl define β~r≔β^0m−1,β^tm−1,β^rm−1T with β^rm−1 denoting the rth component of β^lm−1. Compute the score vector and Fisher matrixsrβ~r=∂ℓpen/∂β~r,Frβ~r=−E∂2ℓpen/∂β~r∂β~rT,with respect to the current interceptβ^0m−1, time effect β^tm−1, and the rth linear effect β^rm−1. Obtain pl possible updatesur=Frβ~r−1srβ~r,and find the best performing component∗∈1,⋯,pl according to Section 2.2, yielding the update u∗=u0,ut,u∗T containing the update for the effect ∗ with corresponding updates u0 for intercept and ut for the time effect. Receive β^0m, β^tm, and β^lm by updatingβ^0m=β^0m−1+νlu0,β^tm=β^tm−1+νlut,β^rm=β^rm−1ifr≠∗,β^rm−1+νlu∗ifr=∗,r=1,⋯,pl.step2: Update random effectsReceive updatesγ^0m−1⟶γ^0m,γ^tm−1⟶γ^tm,for random intercepts and slopes in an additional Fisher scoring step on the penalized log-likelihoodℓpen.step3: Update variance-covariance componentsUpdate variance-covariance componentsσ^2m−1⟶σ^2m,Q^m−1⟶Q^mfollowing the description in Section3.3.end forProceed with estimatesβ^0ml, β^tml, β^lml, γ^0ml, γ^tml as fixed values.Survival partform=1 to msdoUpdate survival effectsForr=1,⋯,ps+1 define β~r≔λ^m−1T,β^rm−1T with β^rm−1 denoting the rth component of β^sm−1 and β^ps+1m−1=α^m−1. Compute the score vector and Fisher matrixsrβ~r=∂ℓpen/∂β~r,Frβ~r=−E∂2ℓpen/∂β~r∂β~rT,with respect to the current baseline hazardλ^m−1 and the rth linear effect βrm−1 or α^m−1, respectively. Obtain ps+1 possible updatesur=Frβ~r−1srβ~r,and find the best performing effect∗∈1,⋯,ps+1 according to Section 2.2, yielding the update u∗=uλT,u∗T containing the update for the effect ∗ with corresponding baseline hazard update uλ. Receive λ^m, β^sm, and α^m by updatingλ^m=λ^m−1+νsuλ,α^m=α^m−1,if∗≠ps+1,α^m−1+νsu∗,if∗=ps+1,β^rm=β^rm−1ifr≠∗,β^rm−1+νlu∗ifr=∗,r=1,⋯,ps.end for(iii)Determine the best performing tuple m∗,l,m∗,s with respect to prediction based on the unpenalized joint log-likelihood ℓ as explained in more detail in Section 3.3.


### 3.3. Computational Details of the Algorithm

In general, the algorithm carries out a hybrid between a two-stage and a joint likelihood approach. For one single tupleml,ms, the fitting procedure goes as follows: In a first step, the longitudinal submodel is boosted up to ml iterations using the lbbLMM boosting algorithm [25]. The received estimates are carried forward into the survival model, where another boosting process up to ms iterations takes place. This fitting process is carried out for any tuples ml,ms with ml<mmax,l, ms<mmax,s, where mmax,l and mmax,s are prespecified maximum numbers of iterations per submodel. For every of these possible combinations of stopping iterations, the corresponding estimates are evaluated based on the joint likelihood using test data, which can be achieved via cross-validation or bootstrapping. Hence, the algorithm uses two-stage fitting but joint likelihood evaluation. We give a detailed description for both parts. Exact formulas for all appearing variants of score vectors and Fisher matrices can be found in the supplementary material. Please note that due to the component-wise updating scheme in both submodels, the lbbJM algorithm works with arbitrarily high numbers of candidate variables and is therefore not confined to low-dimensional data structures.For starting values, the parameters, which actually underlie the boosting process, are necessarily set to zero, thusβ^l0=β^s0=α^0=0. The baseline hazard is initialized with the intercept estimator λ^0=∑iδi/∑iTik=1,⋯,K. The remaining values are extracted from a standard linear mixed model for time and random effects
(15)yij=β0+βttij+γ0i+γtitij+εij,by using, e.g., the function lme from the R package nlme.For boosting longitudinal fixed effects, for convenience, we omit the iteration index as well as the hat indicating estimated values in the following subsections. In the first step of the longitudinal part, the effectsβl follow the classical component-wise likelihood-based boosting procedure. In each iteration, an update for every single covariate r together with intercept β0 and time effect βt is computed, leading to pl three-dimensional updates. Selection of the best performing component is then performed either by selecting the component yielding the optimal likelihood maximization or lowest information criteria like AIC or BIC, which minimizes the model complexity rather than residuals. The linear effect βt is excluded from the selection process, as it holds a very important role in a joint model and should always be included.For updating random effects, after the best performing fixed effect fromβl was updated, in every iteration, an update for the random effects is executed separately. This means that the score vector and Fisher matrix
(16)sranγ=∂ℓpen∂γ,Franγ=−E∂2ℓpen∂γ∂γT,have to be derived in order to execute the update
(17)γnew=γ+CFranγ−1sranγ.The matrixC is a correction matrix which prevents from potential correlations between the random effects estimates and any observed covariates [25], and its derivation can be traced in more detail in the supplementary material.For updating variance-covariance components, the covariance matrixQ of the random effects is updated with an approximate EM algorithm using the posterior curvatures Fii of the random effects model [30]. Receive an update by computing
(18)Q=1n∑i=1nFii−1+γiγiT.The current longitudinal model error is obtained by settingVary−ηl.For boosting the association parameter and survival effects, once the longitudinal part was updated in up toml iterations, the algorithm proceeds to boost the effects βs and α. Although being of different structures, the association parameter α is boosted alongside the effects βs, meaning the algorithm decides whether the association or some baseline survival effect is updated, based on which parameter leads to the best likelihood improvement. This means, the linear effect of the whole longitudinal trajectory is also scaled by the step length ν when being updated within the selection step, which minimizes the chance of potential overfitting also for the association parameter. An alternative method would be choosing just from the effects in βs and updating α in an additional step by optimizing the current likelihood. This approach was used in [11] and treats α as a nuisance parameter, which might not be satisfactory with regards to the importance of α. Again, only the update for α and βs is scaled by the step length νs. The baseline hazard λ receives a full update.For step lengths, apart from the stopping iterations, the step lengths are tuning parameters of the boosting algorithm. Although there is some effort in focusing on adaptive step lengths recently, we chose to set both step lengths to the constant valueνl=νs=0.1. The exact choice of the step length factor is of minor importance as long as it is sufficiently small to ensure proper performance. Setting it to 0.1 is an established choice in the boosting literature [31, 32].For stopping iterations, since the step lengths are chosen to be constant, the tupleml,ms is the main tuning parameter of the boosting algorithm. In regular boosting with only one iteration index, it is convenient to check for every single iteration and take the estimates from the estimation count leading to the best prediction. In the present two-dimensional case, this would mean finding
(19)m∗,l,m∗,s=argmaxml,ms∈Mℓϑ^ml,ms∣Xtest,with M≔1,⋯,mmax,l×1,⋯,mmax,s, ϑ^ml,ms denoting the vector of estimates received via the tuple ml,ms of total iterations and Xtest a complete set of test data for evaluation. Problem (19) is then solved via k-fold cross-validation. But since checking for every single tuple ml,ms∈M implies a very high computational effort, we suggest to coarsen the grid and check for fewer possible stopping iterations in the longitudinal part, e.g., ml∈10,20,30,⋯,mmax,l. Because of the two-stage-approach nature of the algorithm, we still can check for every single ms∈1,⋯,mmax,s without gaining additional computational effort. Furthermore, parallel computing can be executed in order to minimize computational demand.


## 3.1. The Model

Before we introduce the boosting algorithm, we describe the specific joint model. Withi=1,⋯,n denoting the individual and j=1,⋯,ni a single specific measurement, the longitudinal submodel is given by
(8)yij=ηltij,xli+εij=β0+βttij+βlTxli+γ0i+γtitij+εij,where yij is modelled by ηl depending on specific measurement times tij and longitudinal time-independent covariates xli∈ℝpl. This represents a standard linear mixed model with intercept β0 and fixed linear effects βt and βl of time and baseline covariates as well as individual specific random effects γ0i and γti with γ0i,γti~N⊗20,Q. The error terms εij are assumed to follow a normal distribution with Eεij=0 and Varεij=σ2>0.Please note that the model can be additionally extended to interaction effects of timetij with baseline covariates xti∈ℝpt by including the term βtTxtitij in (8). This results in slightly adjusted integrals in the survival part and is omitted in the following for the sake of better readability.In the survival part, the individual hazard(9)λit=λ0texpηsxsi+αηlt,xli,is modelled with the survival predictor ηsxsi=βsTxsi containing additional linear effects βs of baseline covariates xsi∈ℝps. To execute a full likelihood approach, the baseline hazard
(10)λ0t=∑k=1Kλk1Ikt,is chosen to be piecewise-constant depending on the number of segments K and their exact locations Ik=tk−1,tk with t0=0 and tk=maxT for k=1,⋯,K. The collection of values for the baseline hazard is denoted in λ=λkk=1,⋯,K. Later, we will choose K between 7 and 10 in order to guarantee substantially more flexibility than a constant baseline hazard without becoming computationally too demanding.Given two formulas (8) and (9), we can now calculate the joint log-likelihood. Let y=yiji=1,⋯,n,j=1,⋯,ni denote the collection of all longitudinal measurements. Assuming the time-to-event process is conditionally independent from the longitudinal random structure, the joint likelihood can be decomposed into a longitudinal and a survival term. Set γT=γii=1,⋯,n with γiT=γ0i,γti and ϑl≔β0,βlT,βt,γTT. Furthermore, τ contains information on variance-covariance components σ2 and Q. The unpenalized longitudinal log-likelihood is
(11)ℓlϑl,σ2∣y=∑i=1n∑j=1nilogϕyij∣ηltij,xli,σ2,where ϕ·∣m,v denotes the density of a normal distribution with mean m and variance v. Laplace approximation follows [29] and then leads to an additional quadratic penalty term for the random effects yielding the penalized log-likelihood:
(12)ℓlpenϑl,τ∣y=∑i=1n∑j=1nilogϕyij∣ηltij,xli,σ2−12γiTQγi.Note that for the penalized log-likelihoodτ substitutes σ2 as an argument, since the penalized log-likelihood additionally contains information of the variance matrix Q.On the other hand, for given survival outcomeT,δ=Ti,δii=1,⋯,n with event times Ti and censoring indicator δi, the survival log-likelihood takes the well-known form:
(13)ℓsϑs∣T,δ=∑i=1nδilogλiTi∣ηl,ηs,α,λ−∫0Tiexpλit∣ηl,ηs,α,λdt,with ϑs≔λT,α,βsTT. For ϑ≔ϑlT,ϑsTT, we finally receive the penalized and unpenalized joint log-likelihood:
(14)ℓϑ,σ2∣y,T,δ=ℓlϑl,σ2∣y+ℓsϑs∣T,δ,ℓpenϑ,τ∣y,T,δ=ℓlpenϑl,τ∣y+ℓsϑs∣T,δ.


## 3.2. The Algorithm

The following lbbJM algorithm (likelihood-based boosting for joint models) describes a way to fit the formulated joint model by likelihood-based boosting methods discussed in Section2.2.Algorithm 1: lbbJM.
(i)Initialize starting values β^00, β^t0, β^l0, λ^0, β^s0, α^0, γ^00, and γ^t0 with variance-covariance-components σ^20,Q^0 according to Section 3.3. Choose mmax,l,mmax,s and execute the following for every possible tuple ml,ms with ml<mmax,l, ms<mmax,s:Longitudinal part(ii)form=1 to mldostep1: Update fixed effectsForr=1,⋯,pl define β~r≔β^0m−1,β^tm−1,β^rm−1T with β^rm−1 denoting the rth component of β^lm−1. Compute the score vector and Fisher matrixsrβ~r=∂ℓpen/∂β~r,Frβ~r=−E∂2ℓpen/∂β~r∂β~rT,with respect to the current interceptβ^0m−1, time effect β^tm−1, and the rth linear effect β^rm−1. Obtain pl possible updatesur=Frβ~r−1srβ~r,and find the best performing component∗∈1,⋯,pl according to Section 2.2, yielding the update u∗=u0,ut,u∗T containing the update for the effect ∗ with corresponding updates u0 for intercept and ut for the time effect. Receive β^0m, β^tm, and β^lm by updatingβ^0m=β^0m−1+νlu0,β^tm=β^tm−1+νlut,β^rm=β^rm−1ifr≠∗,β^rm−1+νlu∗ifr=∗,r=1,⋯,pl.step2: Update random effectsReceive updatesγ^0m−1⟶γ^0m,γ^tm−1⟶γ^tm,for random intercepts and slopes in an additional Fisher scoring step on the penalized log-likelihoodℓpen.step3: Update variance-covariance componentsUpdate variance-covariance componentsσ^2m−1⟶σ^2m,Q^m−1⟶Q^mfollowing the description in Section3.3.end forProceed with estimatesβ^0ml, β^tml, β^lml, γ^0ml, γ^tml as fixed values.Survival partform=1 to msdoUpdate survival effectsForr=1,⋯,ps+1 define β~r≔λ^m−1T,β^rm−1T with β^rm−1 denoting the rth component of β^sm−1 and β^ps+1m−1=α^m−1. Compute the score vector and Fisher matrixsrβ~r=∂ℓpen/∂β~r,Frβ~r=−E∂2ℓpen/∂β~r∂β~rT,with respect to the current baseline hazardλ^m−1 and the rth linear effect βrm−1 or α^m−1, respectively. Obtain ps+1 possible updatesur=Frβ~r−1srβ~r,and find the best performing effect∗∈1,⋯,ps+1 according to Section 2.2, yielding the update u∗=uλT,u∗T containing the update for the effect ∗ with corresponding baseline hazard update uλ. Receive λ^m, β^sm, and α^m by updatingλ^m=λ^m−1+νsuλ,α^m=α^m−1,if∗≠ps+1,α^m−1+νsu∗,if∗=ps+1,β^rm=β^rm−1ifr≠∗,β^rm−1+νlu∗ifr=∗,r=1,⋯,ps.end for(iii)Determine the best performing tuple m∗,l,m∗,s with respect to prediction based on the unpenalized joint log-likelihood ℓ as explained in more detail in Section 3.3.


## 3.3. Computational Details of the Algorithm

In general, the algorithm carries out a hybrid between a two-stage and a joint likelihood approach. For one single tupleml,ms, the fitting procedure goes as follows: In a first step, the longitudinal submodel is boosted up to ml iterations using the lbbLMM boosting algorithm [25]. The received estimates are carried forward into the survival model, where another boosting process up to ms iterations takes place. This fitting process is carried out for any tuples ml,ms with ml<mmax,l, ms<mmax,s, where mmax,l and mmax,s are prespecified maximum numbers of iterations per submodel. For every of these possible combinations of stopping iterations, the corresponding estimates are evaluated based on the joint likelihood using test data, which can be achieved via cross-validation or bootstrapping. Hence, the algorithm uses two-stage fitting but joint likelihood evaluation. We give a detailed description for both parts. Exact formulas for all appearing variants of score vectors and Fisher matrices can be found in the supplementary material. Please note that due to the component-wise updating scheme in both submodels, the lbbJM algorithm works with arbitrarily high numbers of candidate variables and is therefore not confined to low-dimensional data structures.For starting values, the parameters, which actually underlie the boosting process, are necessarily set to zero, thusβ^l0=β^s0=α^0=0. The baseline hazard is initialized with the intercept estimator λ^0=∑iδi/∑iTik=1,⋯,K. The remaining values are extracted from a standard linear mixed model for time and random effects
(15)yij=β0+βttij+γ0i+γtitij+εij,by using, e.g., the function lme from the R package nlme.For boosting longitudinal fixed effects, for convenience, we omit the iteration index as well as the hat indicating estimated values in the following subsections. In the first step of the longitudinal part, the effectsβl follow the classical component-wise likelihood-based boosting procedure. In each iteration, an update for every single covariate r together with intercept β0 and time effect βt is computed, leading to pl three-dimensional updates. Selection of the best performing component is then performed either by selecting the component yielding the optimal likelihood maximization or lowest information criteria like AIC or BIC, which minimizes the model complexity rather than residuals. The linear effect βt is excluded from the selection process, as it holds a very important role in a joint model and should always be included.For updating random effects, after the best performing fixed effect fromβl was updated, in every iteration, an update for the random effects is executed separately. This means that the score vector and Fisher matrix
(16)sranγ=∂ℓpen∂γ,Franγ=−E∂2ℓpen∂γ∂γT,have to be derived in order to execute the update
(17)γnew=γ+CFranγ−1sranγ.The matrixC is a correction matrix which prevents from potential correlations between the random effects estimates and any observed covariates [25], and its derivation can be traced in more detail in the supplementary material.For updating variance-covariance components, the covariance matrixQ of the random effects is updated with an approximate EM algorithm using the posterior curvatures Fii of the random effects model [30]. Receive an update by computing
(18)Q=1n∑i=1nFii−1+γiγiT.The current longitudinal model error is obtained by settingVary−ηl.For boosting the association parameter and survival effects, once the longitudinal part was updated in up toml iterations, the algorithm proceeds to boost the effects βs and α. Although being of different structures, the association parameter α is boosted alongside the effects βs, meaning the algorithm decides whether the association or some baseline survival effect is updated, based on which parameter leads to the best likelihood improvement. This means, the linear effect of the whole longitudinal trajectory is also scaled by the step length ν when being updated within the selection step, which minimizes the chance of potential overfitting also for the association parameter. An alternative method would be choosing just from the effects in βs and updating α in an additional step by optimizing the current likelihood. This approach was used in [11] and treats α as a nuisance parameter, which might not be satisfactory with regards to the importance of α. Again, only the update for α and βs is scaled by the step length νs. The baseline hazard λ receives a full update.For step lengths, apart from the stopping iterations, the step lengths are tuning parameters of the boosting algorithm. Although there is some effort in focusing on adaptive step lengths recently, we chose to set both step lengths to the constant valueνl=νs=0.1. The exact choice of the step length factor is of minor importance as long as it is sufficiently small to ensure proper performance. Setting it to 0.1 is an established choice in the boosting literature [31, 32].For stopping iterations, since the step lengths are chosen to be constant, the tupleml,ms is the main tuning parameter of the boosting algorithm. In regular boosting with only one iteration index, it is convenient to check for every single iteration and take the estimates from the estimation count leading to the best prediction. In the present two-dimensional case, this would mean finding
(19)m∗,l,m∗,s=argmaxml,ms∈Mℓϑ^ml,ms∣Xtest,with M≔1,⋯,mmax,l×1,⋯,mmax,s, ϑ^ml,ms denoting the vector of estimates received via the tuple ml,ms of total iterations and Xtest a complete set of test data for evaluation. Problem (19) is then solved via k-fold cross-validation. But since checking for every single tuple ml,ms∈M implies a very high computational effort, we suggest to coarsen the grid and check for fewer possible stopping iterations in the longitudinal part, e.g., ml∈10,20,30,⋯,mmax,l. Because of the two-stage-approach nature of the algorithm, we still can check for every single ms∈1,⋯,mmax,s without gaining additional computational effort. Furthermore, parallel computing can be executed in order to minimize computational demand.


## 4. Simulations

We evaluate the lbbJM algorithm with a simulation study. The aim is to assess estimation and shrinkage characteristics in general as well as variable selection properties and performance in high dimensional, i.e.,p>n settings. The lbbJM algorithm is included in two variants. While lbbJMa executes the full approach as depicted in Section 3.2, lbbJMb performs a shortened two-stage procedure where the longitudinal submodel is fitted in advance using regular maximum likelihood inference and does not underlie any regularization. The exact lbbJMb algorithm is depicted in detail in the supplementary material. We additionally include the JM package as state of the art for convenient estimation of joint models as well as the glmnet package, which offers elastic net penalization for start-stop-data and therefore an alternative approach for regularization of time-dependent covariates in survival analysis. None of the competitors are completely suitable for a benchmark comparison and are viewed as reference points for the specific objectives addressed by the lbbJM algorithm. Regarding glmnet, as an alternative approach to regularization of time-dependent covariates, shrinkage and variable selection properties are of interest, although it focuses solely on survival analysis. JM in addition offers unregularized effect estimates with corresponding significance indicator but is neither suitable for high-dimensional setups nor offers variable selection.


### 4.1. Setup

The simulations are executed according to the model described in Section3.1 with n∈100,500 and ni=5 using inversion sampling, which is explained in detail in the supplementary material. The prespecified true parameter values are
(20)β0=1,βt=2,βlT=2,1,2,βsT=1,2,−1,α=0.5,with variance components
(21)σ=0.1,Q=20.10.10.3.The entries of the covariate vectorsxli and xsi are drawn independently from the standard normal distribution N0,1. In addition to the informative covariates with effects βl and βs, the covariate vectors xli and xsi are expanded with noninformative noise variables until the chosen numbers pl and ps of total dimensions are reached. These additional noise variables are included to evaluate variable selection properties and robustness of the approach in case of a misspecified model. The baseline hazard is chosen as λ0t=2.5t1.5 and given the censoring mechanism described in the supplementary material, the chosen parameter values result in an average censoring rate of ≈50%.Overall, we consider two scenarios. One low-dimensional setup withn=500 and pl=ps=9 mimicking a more common data structure and one high-dimensional setup, where the number of covariates included in the survival submodel exceeds the number of individuals so that conventional approaches like JM fail to return results.For the computation, JM and lbbJM useK=10 with equidistant knot placement. The grid M≔25,50,⋯,500×1,2,⋯,1000 is specified for possible tuples of stopping iterations and the optimal regularization parameter for glmnet is determined by the function cv.glmnet() based on 10-fold cross-validation.


### 4.2. Results

Since the compared estimation routines follow different approaches targeting various objectives from regular maximum likelihood estimation in joint models to regularization in pure time-to-event analysis, we focus on plain coefficient estimates averaged over 100 independent simulation runs in order to asses estimation characteristics. Variable selection properties are evaluated by considering share of true positives (TP) and false discovery rate (FDR).For low-dimensional setup (n=500, ps=9), Table 1 depicts the results for the low-dimensional setup. In the longitudinal submodel, lbbJMa has small shrinkage and therefore offers variable selection with a rather low false discovery rate of 0.23 but still selects all informative variables. The time effect βt receives comparatively high shrinkage, since time, as a cluster-varying variable, adds more information to the model. Overall, the longitudinal submodel is boosted up to 108.25 iterations on average yielding rather strongly shrunk coefficient estimates. Please note that the results for lbbJMb are simply obtained by lme() and very similar to JM. In the survival part, both boosting approaches substantially outperform glmnet in terms of variable selection while again receiving also more shrinkage. Due to the comparatively rough baseline hazard depicted in Figure 1, all full likelihood approaches, i.e., JM and lbbJM, receive additional shrinkage which is unaffected by possible regularization. As glmnet uses the partial likelihood, there are no estimates for the baseline hazard function available and the small elastic net penalty of λ∗=0.003 on average also results in weaker performance regarding the rate of false positives. Overall, the lbbJM approaches yield satisfactory results regarding both regularization and variable selection. The effect estimates clearly reflect the true values approximately obtained by JM while simultaneously receiving sufficiently large shrinkage for decent performance of identifying influential covariates in both the longitudinal and survival submodels.Table 1
Shrinkage and variable selection properties regarding longitudinal and survival outcomes averaged over 100 simulation runs of the low-dimensional scenario.
βt (sd)βl1 (sd)βl2 (sd)βl3 (sd)TPFDRml∗True2121JM1.9980.9942.0081.002———(0.03)(0.07)(0.07)(0.07)lbbJMa1.7600.9141.9220.9231.000.23108.25(0.08)(0.07)(0.07)(0.07)lbbJMb1.9920.9942.0081.002———(0.03)(0.07)(0.07)(0.07)α (sd)βs1 (sd)βs2 (sd)βs3 (sd)TPFDRms∗True0.512-2JM0.4570.9031.807-1.800———(0.04)(0.08)(0.12)(0.12)lbbJMa0.3900.7281.521-1.5161.000.27209.2(0.03)(0.06)(0.07)(0.07)lbbJMb0.3730.7131.500-1.4951.000.22196.9(0.03)(0.06)(0.07)(0.08)glmnet0.4270.9091.833-1.8231.000.51—(0.03)(0.07)(0.11)(0.10)Figure 1
Piecewise-constant baseline hazard estimates withK=10 by JM, lbbJMa and lbbJMb averaged over 100 simulation runs of the low dimensional scenario.For high-dimensional setup (n=100, ps=100), Table 2 depicts the results for the high-dimensional setup. As expected, estimates in the high-dimensional setup are regularized stronger and therefore experience more shrinkage. Again, the boosting approaches contained in lbbJM yield better variable selection properties and slightly more regularized coefficient estimates, although results seem to align with increasing dimensions. Note that JM is not capable of handling high-dimensional data structures and is therefore not included in the high-dimensional setup at all.Table 2
Shrinkage and variable selection properties regarding the longitudinal and survival outcomes averaged over 100 simulation runs of the high-dimensional scenario.
βt (sd)βl1 (sd)βl2 (sd)βl3 (sd)TPFDRml∗True2121lbbJMa1.7480.8681.8430.8751.000.36124.2(0.20)(0.14)(0.14)(0.14)lbbJMb1.9911.0081.9821.011———(0.08)(0.13)(0.16)(0.15)α (sd)βs1 (sd)βs2 (sd)βs3 (sd)TPFDRms∗True0.512-2lbbJMa0.3070.5121.242-1.2161.000.70136.7(0.06)(0.12)(0.16)(0.13)lbbJMb0.2850.4981.215-1.1911.000.67127.0(0.05)(0.13)(0.15)(0.13)glmnet0.2930.6271.449-1.4221.000.83—(0.08)(0.18)(0.30)(0.28)For computational effort, Table3 shows estimates for elapsed computation time of each routine. Times are measured in seconds and depict the computation time for one single model fit which was carried out on a 2×2.66 GHz-6-Core Intel Xeon CPU (64 GB RAM). As expected, the full boosting approach executed in lbbJMa comes with high computational costs similar to [11]. Note that the runtimes are higher in the low-dimensional scenario as the overall number of clusters is higher (n=500) leading to an increased burden in the already time-consuming longitudinal boosting procedure. The reduced approach lbbJMb, however, runs considerably faster and is therefore more desirable as long as research focus lies solely on the time-to-event analysis.Table 3
Averaged computation times for one single model fit (in seconds).
SetupJMglmnetlbbJMalbbJMbLow110.00149.1515776.1643.76High—156.444072.80248.08


## 4.1. Setup

The simulations are executed according to the model described in Section3.1 with n∈100,500 and ni=5 using inversion sampling, which is explained in detail in the supplementary material. The prespecified true parameter values are
(20)β0=1,βt=2,βlT=2,1,2,βsT=1,2,−1,α=0.5,with variance components
(21)σ=0.1,Q=20.10.10.3.The entries of the covariate vectorsxli and xsi are drawn independently from the standard normal distribution N0,1. In addition to the informative covariates with effects βl and βs, the covariate vectors xli and xsi are expanded with noninformative noise variables until the chosen numbers pl and ps of total dimensions are reached. These additional noise variables are included to evaluate variable selection properties and robustness of the approach in case of a misspecified model. The baseline hazard is chosen as λ0t=2.5t1.5 and given the censoring mechanism described in the supplementary material, the chosen parameter values result in an average censoring rate of ≈50%.Overall, we consider two scenarios. One low-dimensional setup withn=500 and pl=ps=9 mimicking a more common data structure and one high-dimensional setup, where the number of covariates included in the survival submodel exceeds the number of individuals so that conventional approaches like JM fail to return results.For the computation, JM and lbbJM useK=10 with equidistant knot placement. The grid M≔25,50,⋯,500×1,2,⋯,1000 is specified for possible tuples of stopping iterations and the optimal regularization parameter for glmnet is determined by the function cv.glmnet() based on 10-fold cross-validation.


## 4.2. Results

Since the compared estimation routines follow different approaches targeting various objectives from regular maximum likelihood estimation in joint models to regularization in pure time-to-event analysis, we focus on plain coefficient estimates averaged over 100 independent simulation runs in order to asses estimation characteristics. Variable selection properties are evaluated by considering share of true positives (TP) and false discovery rate (FDR).For low-dimensional setup (n=500, ps=9), Table 1 depicts the results for the low-dimensional setup. In the longitudinal submodel, lbbJMa has small shrinkage and therefore offers variable selection with a rather low false discovery rate of 0.23 but still selects all informative variables. The time effect βt receives comparatively high shrinkage, since time, as a cluster-varying variable, adds more information to the model. Overall, the longitudinal submodel is boosted up to 108.25 iterations on average yielding rather strongly shrunk coefficient estimates. Please note that the results for lbbJMb are simply obtained by lme() and very similar to JM. In the survival part, both boosting approaches substantially outperform glmnet in terms of variable selection while again receiving also more shrinkage. Due to the comparatively rough baseline hazard depicted in Figure 1, all full likelihood approaches, i.e., JM and lbbJM, receive additional shrinkage which is unaffected by possible regularization. As glmnet uses the partial likelihood, there are no estimates for the baseline hazard function available and the small elastic net penalty of λ∗=0.003 on average also results in weaker performance regarding the rate of false positives. Overall, the lbbJM approaches yield satisfactory results regarding both regularization and variable selection. The effect estimates clearly reflect the true values approximately obtained by JM while simultaneously receiving sufficiently large shrinkage for decent performance of identifying influential covariates in both the longitudinal and survival submodels.Table 1
Shrinkage and variable selection properties regarding longitudinal and survival outcomes averaged over 100 simulation runs of the low-dimensional scenario.
βt (sd)βl1 (sd)βl2 (sd)βl3 (sd)TPFDRml∗True2121JM1.9980.9942.0081.002———(0.03)(0.07)(0.07)(0.07)lbbJMa1.7600.9141.9220.9231.000.23108.25(0.08)(0.07)(0.07)(0.07)lbbJMb1.9920.9942.0081.002———(0.03)(0.07)(0.07)(0.07)α (sd)βs1 (sd)βs2 (sd)βs3 (sd)TPFDRms∗True0.512-2JM0.4570.9031.807-1.800———(0.04)(0.08)(0.12)(0.12)lbbJMa0.3900.7281.521-1.5161.000.27209.2(0.03)(0.06)(0.07)(0.07)lbbJMb0.3730.7131.500-1.4951.000.22196.9(0.03)(0.06)(0.07)(0.08)glmnet0.4270.9091.833-1.8231.000.51—(0.03)(0.07)(0.11)(0.10)Figure 1
Piecewise-constant baseline hazard estimates withK=10 by JM, lbbJMa and lbbJMb averaged over 100 simulation runs of the low dimensional scenario.For high-dimensional setup (n=100, ps=100), Table 2 depicts the results for the high-dimensional setup. As expected, estimates in the high-dimensional setup are regularized stronger and therefore experience more shrinkage. Again, the boosting approaches contained in lbbJM yield better variable selection properties and slightly more regularized coefficient estimates, although results seem to align with increasing dimensions. Note that JM is not capable of handling high-dimensional data structures and is therefore not included in the high-dimensional setup at all.Table 2
Shrinkage and variable selection properties regarding the longitudinal and survival outcomes averaged over 100 simulation runs of the high-dimensional scenario.
βt (sd)βl1 (sd)βl2 (sd)βl3 (sd)TPFDRml∗True2121lbbJMa1.7480.8681.8430.8751.000.36124.2(0.20)(0.14)(0.14)(0.14)lbbJMb1.9911.0081.9821.011———(0.08)(0.13)(0.16)(0.15)α (sd)βs1 (sd)βs2 (sd)βs3 (sd)TPFDRms∗True0.512-2lbbJMa0.3070.5121.242-1.2161.000.70136.7(0.06)(0.12)(0.16)(0.13)lbbJMb0.2850.4981.215-1.1911.000.67127.0(0.05)(0.13)(0.15)(0.13)glmnet0.2930.6271.449-1.4221.000.83—(0.08)(0.18)(0.30)(0.28)For computational effort, Table3 shows estimates for elapsed computation time of each routine. Times are measured in seconds and depict the computation time for one single model fit which was carried out on a 2×2.66 GHz-6-Core Intel Xeon CPU (64 GB RAM). As expected, the full boosting approach executed in lbbJMa comes with high computational costs similar to [11]. Note that the runtimes are higher in the low-dimensional scenario as the overall number of clusters is higher (n=500) leading to an increased burden in the already time-consuming longitudinal boosting procedure. The reduced approach lbbJMb, however, runs considerably faster and is therefore more desirable as long as research focus lies solely on the time-to-event analysis.Table 3
Averaged computation times for one single model fit (in seconds).
SetupJMglmnetlbbJMalbbJMbLow110.00149.1515776.1643.76High—156.444072.80248.08


## 5. Application

We showcase the lbbJM algorithm by applying it to the 1994 AIDS data [26]. The study is aimed at comparing the two antiretroviral drugs, didanosine (ddI) and zalcitabine (ddC), based on a collective of 467 patients infected with human immune deficiency virus (HIV) who were either intolerant to or failed a previous treatment with Zidovudine (AZT). Alongside several baseline covariates, the square root CD4 cell count was recorded at study entry and in multiple follow-ups after 2, 6, 12, and 18 months, respectively. The CD4 cells are attacked by the virus and thus decrease over time for infected patients; hence, they are a widely used surrogate for disease progression. In addition to the longitudinal outcome, 188 patients died during the time of the study leading to 188 observed and 279 censored events. The structure of the data is depicted in Table 4.Table 4
Structure of the dataset with primary outcomes for the joint analysis in the three columns on the left.
yTδtDrugGenderAZTprevOIID10.6716.9700ddCMaleIntoleranceAIDS18.4316.9706ddCMaleIntoleranceAIDS19.4316.97012ddCMaleIntoleranceAIDS16.3219.0000ddIMaleIntolerancenoAIDS28.1219.0006ddIMaleIntolerancenoAIDS24.5819.00012ddIMaleIntolerancenoAIDS25.0019.00018ddIMaleIntolerancenoAIDS23.4618.5300ddIFemaleIntoleranceAIDS33.6118.5302ddIFemaleIntoleranceAIDS36.1618.5316ddIFemaleIntoleranceAIDS3⋮⋮⋮⋮⋮⋮⋮⋮⋮We formulate the joint model(22)CD4it=ηlit+εit=β0+βldrugi+βt1t+βt2t·drugi+γ0i+γtit+εit,(23)λit=λ0texpβs1·genderi+βs2·AZTi+βs3·prevOIi+αηlit.The CD4 cell countCD4it for i=1,⋯,467 is described by a linear mixed model with random intercepts, random slopes for time t and linear effects of time, drug and an additional interaction between time and drug. The true, i.e. modelled by ηlt, underlying profile of the CD4 cell count is then included together with the remaining baseline covariates in the Cox full likelihood model, thus the model error εt is eliminated. Here, drug is a dummy for ddI, gender for female gender, AZT for failure of Zidovudine therapy, and prevOI for prevalence of AIDS. The number of segments for the baseline hazard was chosen to be K=7.We fit model (22) with the same methods as already used in the simulation study. The tuning parameters of the boosting algorithm were chosen to be νl=νs=0.1 for step lengths and Ma≔5,10,⋯,100×1,2,⋯,250 for the grid of possible tuples of stopping iterations for lbbJMa and Mb≔1,2,⋯,250 for lbbJMb. Again, glmnet was tuned using the cv.glmnet() function and all regularization approaches are based on 10-fold cross validation.For lbbJMa, m∗,l=10 and m∗,s=33 formed the best performing tuple of stopping iterations. The two-stage approach lbbJMb used m∗,s=40 and the optimal penalization parameter for glmnet turned out as λpen∗=0.0014. The corresponding coefficient estimates are shown in Table 5. Overall, the results reflect what was already observed in the simulation study. While glmnet shows quite conservative shrinkage properties where every variable is included in the final model, the lbbJM approaches tend to stop rather early yielding bigger shrinkage and in addition effects, which did not get selected at all. In order to give some point of reference regarding variable selection properties, the p values computed by JM are included in Table 5. The selected effects align quite nicely with being significant according to JM. While lbbJMa only includes the variable drug with only half the effect size, lbbJMb additionally sees a quite tiny impact of female gender.Table 5
Shrinkage and variable selection properties by the different packages for model (22).
β0βlβt1βt2αβs1βs2βs3JM6.970.49−0.18<0.01−0.240.310.090.66lbbJMa6.950.26−0.050−0.13000.73lbbJMb6.950.48−0.16−0.02−0.180.0300.61glmnet————−0.150.310.090.81p value (JM)<0.010.26<0.010.98<0.010.230.61<0.01The corresponding coefficient progressions for the survival submodel are visualized in Figure2. Both algorithms update the coefficients referring to the longitudinal CD4 cell profile and the variable prevOI right away and therefore see a strong connection between the risk for death and the CD4 cell count as well as whether or not a patient has AIDS. Due to early stopping, lbbJMa selects neither of the two remaining covariates into the final model. lbbJMb includes one variable more, gender, which however only has a very small effect.Figure 2
Coefficient progression in the survival part for lbbJMa ((a), with m∗,l=10 and lbbJMb (b).


## 6. Outlook and Discussion

Overall, the lbbJM algorithm introduces a novel boosting-based regularization scheme to joint models focusing on survival analysis as well as to Cox models with time-dependent covariates. The method fits in well among alternative routines and especially stands out with respect to variable selection properties. Due to its clear advantage regarding computational effort as depicted in Table3, lbbJMb is the preferred routine when research interest clearly lies on time-to-event analysis, whereas lbbJMa is capable of regularizing both submodels simultaneously. Besides the good results regarding variable selection, it is also expected that the proposed boosting methods can improve the predictive power of a joint model, since boosting algorithms are, due to their model tuning based on test errors, primarily used for prediction. A thorough investigation of improving and evaluating prediction performance of a joint model by several boosting techniques remains an interesting task.Still, the presented foundation is of a comparatively simple nature and possible extensions include more flexible modelling, e.g., based on P-splines which allow smooth effects of time-dependent as well as time-independent covariates and can additionally include time-varying effects as well as possibly time-varying association structures. As the current algorithm is confined to a time-constant association parameterα, an extension to αt would increase the flexibility of the model. While it is usually difficult to disentangle a time-dependent association parameter from the longitudinal trajectory for parameter estimation, the lbbJM algorithm could potentially avoid these identification issues due to its clear separation in a longitudinal and a survival boosting process.Similar regularization-based approaches have been proved useful for time-to-event settings [22, 33], and it can be assumed that the presented method could benefit from these concepts as well.Moreover, the presented work only lays a foundation to several extensions addressing known issues for both boosting and joint modelling. It represents an accessible way to boosting for time-dependent covariates in survival analysis. Gradient boosting has known limitations in this matter, and although efforts for a framework to overcome these issues are made [34], things rather quickly become technical and the proven robustness of likelihood-based boosting represents a flexible and far more intuitive alternative. Further developments in this direction could include multiple time-dependent covariates based on the two-stage approach of lbbJMb, where additional effort is necessary in order to provide a fair competition between time-dependent and time-independent covariates within the selection procedure.Furthermore, the component-wise updating process is capable of including allocation mechanisms into the algorithm. While the lbbJM algorithm is due to its two-stage fitting process fairly robust to estimate models, where one candidate variable is assigned to both submodels, this kind of specification is not advised in general as identification issues may arise and a proper interpretation of the resulting model might be challenging. However, it is usually tricky to decide, whether a covariate should be included in the longitudinal or the survival part of the model and these decisions are often made using prior knowledge. An allocation routine based on likelihood maximization could therefore greatly improve joint model inference and would additionally eliminate the two-stage nature of the lbbJM algorithm. This would not only decrease the computational burden but also provide a far more flexible algorithm allowing for regularization, variable selection, and allocation.



---

*Source: 4384035-2021-11-15.xml*