# Efficient Algorithms and Implementation of a Semiparametric Joint Model for Longitudinal and Competing Risk Data: With Applications to Massive Biobank Data

Authors: Shanpeng Li; Ning Li; Hong Wang; Jin Zhou; Hua Zhou; Gang Li
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

Semiparametric joint models of longitudinal and competing risk data are computationally costly, and their current implementations do not scale well to massive biobank data. This paper identifies and addresses some key computational barriers in a semiparametric joint model for longitudinal and competing risk survival data. By developing and implementing customized linear scan algorithms, we reduce the computational complexities fromOn2 or On3 to On in various steps including numerical integration, risk set calculation, and standard error estimation, where n is the number of subjects. Using both simulated and real-world biobank data, we demonstrate that these linear scan algorithms can speed up the existing methods by a factor of up to hundreds of thousands when n>104, often reducing the runtime from days to minutes. We have developed an R package, FastJM, based on the proposed algorithms for joint modeling of longitudinal and competing risk time-to-event data and made it publicly available on the Comprehensive R Archive Network (CRAN).


---

## Body

## 1. Introduction

In clinical research and other longitudinal studies, it is common to collect both longitudinal and time-to-event data on each participant, and these two endpoints are often correlated. Joint models of longitudinal and survival data have been widely used to mitigate incorrect estimation and statistical inferences associated with separate analysis of each endpoint [1, 2]. For instance, in longitudinal data analysis, joint models are often used to handle nonignorable missing data due to a terminal event, which cannot be properly accounted for by a standard mixed effects model or generalized estimating equation (GEE) method that relies on the ignorable missing-at-random or missing-completely-at-random assumption [3–5]. Joint models are also popularly employed in survival analysis to study the effects of a time-dependent covariate that is measured intermittently or subject to measurement error [6–9] and for dynamic prediction of an event outcome from the past history of a biomarker [10–14]. Comprehensive reviews of joint models for longitudinal and time-to-event data and their applications can be found in Elashoff et al. [1], Hickey et al. [15], Rizopoulos [16], Sudell et al. [17] and the references therein.Despite the explosive growth of literature on joint models for longitudinal and time-to-event data during the past three decades, efficient implementation of joint models has lagged behind, which limits the application of joint models to only small to moderate studies. Recently, massive sample size data collected from electronic health records (EHRs) and insurance claim databases warrants great opportunities to conduct clinical studies in a real-world setting. For example, the UK Biobank is a prospective cohort study with approximately 500,000 individuals, aged 37-73 years, from the general population between 2006 and 2010 in the United Kingdom [18, 19]. Aggregated and quality controlled EHR data purchasable from Optum (https://www.optum.com/EHR) includes 80+ millions patients with longitudinal lab measures. The Million Veteran Project [20] and IBM MarketScan Database [21] are some of many other big biobank databases that contain rich yet complex longitudinal and time-to-event data on 105~108 patients. However, current implementations of many joint models are inefficient and not scalable to the size of biobank data as demonstrated later in Sections 3 and 4. There is a pressing need to develop efficient implementations of joint models to enable the analysis of these rich data sources.The purpose of this paper is to develop and implement efficient algorithms for a semiparametric joint model of longitudinal and competing risk time-to-event data. As specified in Section2.1, the joint model consists of a linear mixed effects submodel for a longitudinal outcome and a semiparametric proportional cause-specific hazard submodel for a competing risk survival outcome. These two submodels are linked together by shared random effects or features of an individual’s longitudinal biomarker trajectory. In Section 2, we identify key computational bottlenecks in the semiparametric maximum likelihood inference procedure for the joint model. Specifically, we point out that in a standard implementation, the computational complexities for numerical integration, risk set calculation, and standard error estimation are of the order On2, On2, and On3, respectively, where n is the number of subjects. Consequently, current implementation grinds to a halt as n becomes large (for example, n>104). We further develop tailored linear scan algorithms to reduce the computational complexity to On in each of the aforementioned components. We illustrate by simulation and real-world data that the linearization algorithms can result in a drastic speed-up by a factor of many thousands when n>104, reducing the runtime from days to minutes for big data. Finally, we have developed a user-friendly R package FastJM to fit the shared parameter semiparametric joint model using the proposed efficient algorithms and made it publicly available on the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=FastJM.The rest of the paper is organized as follows. Section2.1 outlines the semiparametric shared random effect joint model framework and reviews a customized expectation-maximization (EM) algorithm for semiparametric maximum likelihood estimation as well as a standard error estimation method. Section 2.2 develops various linear scan algorithms to address the key computational bottlenecks in the EM algorithm and standard error estimation for large data. Section 3 presents simulation studies to illustrate the computational efficiency of the proposed linear scan algorithms. Section 4 demonstrates the improved computational performance of our FastJM R package over some established joint model R packages on two moderate to large real-world data sets. Concluding remarks are provided in Section 5.


## 2. Efficient Algorithms for a Semiparametric Joint Model of Longitudinal and Competing Risk Data

### 2.1. Notations and Preliminaries

LetYit be the longitudinal outcome at time t for subject i, i=1,2,⋯,n, and n is the total number of subjects. Suppose the longitudinal outcome Yit is observed at time points tij, j=1,2,⋯,ni, and denote Yi=Yi1,⋯,Yini. Let Ci=Ti,Di be the competing risk data on subject i, where Ti is the observed time to either failure or censoring and Di takes value in 0,1,2,⋯,K, with Di=0 indicating a censored event and Di=k implying the kth type of failure is observed on subject i, k=1,2,⋯,K. The censoring mechanism is assumed to be independent of the failure time.


#### 2.1.1. Model

Consider the following joint model in which the longitudinal outcomeYit is characterized by a linear mixed effects model:
(1)Yit=mit+εit,=Xi1tTβ+X~i1tTbi+εit,and the competing risk survival outcome is modeled by a proportional cause-specific hazard model:
(2)λkt∣Xi2,bi,γk,νk=limh⟶0Pt≤Ti<t+h,Di=k∣Ti≥t,Xi2,bih=λ0ktexpXi2Tγk+νkTbi,k=1,⋯,K.In the longitudinal submodel (1), mit is the mean of the longitudinal outcome at time t, Xi1t and X~i1t are column vectors of possibly time-varying covariates associated with the longitudinal outcome Yit, β represents a p×1 vector of fixed effects of Xi1t, bi~Nq0,Σ denotes a q×1 vector of random effects for X~i1t, εit~N0,σ2 is the measurement error independent of bi, and εit1 is independent of εit2 for any t1≠t2. In the competing risk survival submodel (2), λkt∣Xi2,bi,γk,νk is the conditional cause-specific hazard rate for type k failure at time t, given time-invariant covariates Xi2 and the shared random effects bi, and λ0kt is a completely unspecified baseline cause-specific hazard function for type k failure. The two submodels are linked together by the shared random effects bi, and the strength of association is quantified by the association parameters ν1,⋯,νK.


#### 2.1.2. Semiparametric Maximum Likelihood Estimation via an EM Algorithm

LetΨ=β,σ2,γ,ν,Σ,λ01·,⋯,λ0K· denote all unknown parameters for the joint models (1) and (2), where γ=γ1T,⋯,γKTT and ν=ν1T,⋯,νKTT. Let IDi=k be the competing event indicator for type k failure, which takes the value 1 if the condition Di=k is satisfied and 0 otherwise. The observed-data likelihood for Ψ is then given by
(3)LΨ;Y,C∝∏i=1nfYi,Ci∣Ψ=∏i=1n∫bifYi∣Ci,bi,ΨfCi∣bi,Ψfbi∣Ψdbi∝∏i=1n∫bi∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1Kλkt∣Xi2,bi,γk,νkIDi=kexp−∫0Ti∑k=1Kλkt∣Xi2,bi,γk,νkdt×12πqΣexp−12biTΣ−1bidbi,which follows from the assumption that Yi and Ci are independent conditional on the covariates and the random effects.BecauseΨ involves K unknown hazard functions, finding its maximum likelihood estimate by maximizing the above observed-data likelihood is nontrivial. However, a customized EM algorithm can be derived to compute the maximum likelihood estimate of Ψ by regarding the latent random effects bi as missing data [3]. The complete-data likelihood based on Y,C,b is
(4)LΨ;Y,C,b∝∏i=1n∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1KΔΛ0kTiexpXi2Tγk+νkTbiIDi=k×exp−∑k=1KΛ0kTiexpXi2Tγk+νkTbi×12πqΣexp−12biTΣ−1bi,where Λ0k. is the cumulative baseline hazard function for type k failure and ΔΛ0kTi=Λ0kTi−Λ0kTi−.The EM algorithm iterates between an expectation step (E-step):(5)QΨ;Ψm≡Eb∣Y,C,ΨmmlogLΨ;Y,C,b,and a maximization step (M-step):
(6)Ψm+1=argmaxΨQΨ;Ψm,until the algorithm converges, where Ψm is the estimate of Ψ from the mth iteration. The E-step (5) involves calculating the following expected values across all subjects: Embi, EmbibiT, Emexpνkbi, Embiexpνkbi, and EmbibiTexpνkbi, where
(7)Emhbi=∫hbifbi∣Yi,Ci,Ψmdbi,for any function h·. Furthermore, it can be shown that the M-step (6) has closed-form solutions for the parameters β, σ2, Σ, and Λ0kt, and that the other parameters γ and ν can be updated using the one-step Newton-Raphson method. Details are provided in equations (7.1)–(7.8) of the supplementary materials (available here).


#### 2.1.3. Standard Error Estimation

As discussed in Elashoff et al. [22] (Section 4.1, p.72), standard errors of the parametric components of the semiparametric maximum likelihood estimate can be estimated by profiled likelihood, observed information matrix, or bootstrap. All three methods can be computationally intensive when n is large. Here, we focus on the profiled likelihood-based method and show that its computation can be linearized with respect to n.LetΩ=β,Σ,σ2,γ1,⋯,γK,ν1,⋯,νK denote the parametric component of Ψ and Ω∧ its maximum likelihood estimate. The variance-covariance matrix of Ω∧ can be estimated by inverting the following approximate empirical Fisher information [23–25]:
(8)∑i=1n∇ΩliΩ∧;Y,C∇ΩliΩ∧;Y,CT,where ∇ΩliΩ;Y,C is the observed score vector from the profiled likelihood liΩ;Y,C of Ω on the ith subject by profiling out the baseline hazards. Details of the observed score vector for each parametric component are provided in equations (7.9)–(7.13) of the supplementary materials.


### 2.2. Efficient Algorithms and Implementation of the EM Algorithm and Standard Error Estimation

With naive implementation, multiple quantities in the above E-step, M-step, and standard error estimation will involveOn2 or On3 operations, which become computationally prohibitive at large sample size n. Below we identify these bottlenecks and discuss appropriate linear scan algorithms to reduce their computational complexities to On.


#### 2.2.1. Efficient Implementation of the E-Step

At each EM iteration, the E-step (5) requires calculating integral (7) across all subjects. Below we discuss how to accelerate two commonly used numerical integration methods for evaluating these integrals.(1) Standard Gauss-Hermite Quadrature Rule for Numerical Integration. A commonly used method for numerical evaluation of integral (7) is based on the standard Gauss-Hermite quadrature rule [26]:
(9)Emhbi=∫hbifYi,Ci,bi∣ΨmdbifYi,Ci∣Ψm=∫hbifYi,Ci∣bi,Ψmfbi∣Ψmdbi∫fYi,Ci∣bi,Ψmfbi∣Ψmdbi≈∑t1,t2,⋯,tqπthb~tmfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2,where
(10)fYi,Ci∣b~tm,Ψm=fYi∣b~tm,Ψym×fCi∣b~tm,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTb~tm2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkTmb~tmIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkTmb~tm,with Λ0km. the right-continuous and nondecreasing cumulative baseline hazard function for type k failure at the mth EM iteration as defined in Section 7.1 (equation (7.4)) of the supplementary material and ΔΛ0kmTi the jump size of Λ0km. at Ti, q the dimension of the random effects vector, ∑t1,t2,⋯,tq the shorthand for ∑t1=1nq⋯∑tq=1nq, nq the number of quadrature points, ct=ct1,ct2,⋯,ctqT the abscissas with corresponding weights πt, b~tm=2Σ∧m1/2ct the rescaled alternative abscissas, and Σ∧m1/2 the square root of Σ∧m [3]. However, this method is computationally intensive due to multiple factors. First, it usually requires many quadrature points to approximate an integral with sufficient accuracy because the mode of the integrand is often located in a region different from zero. Second, the computational cost increases exponentially with q because the Cartesian product of the abscissas is used to evaluate the integrand with respect to each random effect. Lastly, the alternative abscissas b~tm need to be recalculated at every EM iteration.(2) Pseudo-Adaptive Gauss-Hermite Quadrature Rule for Numerical Integration. When the number of longitudinal measurements per subject is relatively large, Rizopoulos [27] introduced a pseudo-adaptive Gauss-Hermite quadrature rule for numerical approximation of integral (7), which achieves good approximation accuracy with only a small number (nq) of quadrature points and is thus computationally more efficient. The pseudo-adaptive Gauss-Hermite quadrature rule proceeds as follows. First, fit the linear mixed effects model (1) to extract the empirical Bayes estimates of the random effects and its covariance matrix H~i−1:
(11)b~i=Σ^iX~i1V∧i−1Yi−Xi1Tβ^,(12)H~i−1=Σ^−Σ^iX~i1V∧i−1−V∧i−1Xi1T∑i=1nXi1V∧i−1Xi1T−1Xi1V∧i−1×X~i1TΣ^,where Xi1=Xi1ti1,⋯,Xi1tini, X~i1=X~i1ti1,⋯,X~i1tini, and Vi=i1TΣi1+σ2I. Then, define the alternative abscissas r~t=b~i+2H~i−1/2ct and approximate Emhbi by
(13)Emhbi≈∑t1,t2,⋯,tqπthr~tfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2,where
(14)fYi,Ci∣r~t,Ψm=fYi∣r~t,Ψym×fCi∣r~t,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTr~t2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkmTr~tIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkmTr~t,with the notations defined similarly to Equation (10).A derivation of (13) is provided in Section 7.2 of the supplementary materials. The pseudo-adaptive Gauss-Hermite quadrature rule is computationally appealing because the alternate rescaled quadrature points r~t are computed only once before the EM algorithm and do not need to be updated in the EM algorithm. Additionally, the pseudo-adaptive Gauss-Hermite quadrature rule requires fewer quadrature points than the standard Gauss-Hermite quadrature rule to achieve the same numerical approximation accuracy [27]. For example, our simulation results in the supplementary materials (Section A.4, Table A.1) illustrate that the pseudo-adaptive Gauss-Hermite quadrature rule with nq=6 quadrature points produces almost identical results to the standard Gauss-Hermite quadrature rule with nq=20 quadrature points.Remark 1 (Linear calculation ofH~i−1’s across all subjects).
At the first sight, calculatingH~i−1 (12) across all subjects would involve On2 operations since each H~i−1 involves a summation over n subjects. However, because the same quantity ∑i=1nXi1V^i−1Xi1T appears in every H~i−1, one can precompute this quantity and then use the cached value to calculate H~i−1 across all subjects. This way, one can compute the H~i−1’s across all subjects in On operations. Our simulation study in the supplementary material (Section A.5, Figure A.1) shows that applying this simple linearization algorithm can yield a speed-up by a factor of 10 to 10,000 when n grows from 10 to 105. Our implementation is also significantly faster than a popular R package lme4 [28] with over a speed-up by a factor of 10 to 500 as n grows from 10 to 105.(3) Linear Scan Algorithm for CalculatingfYi,Ci∣r~t,Ψm across All subjects. Both the standard Gauss-Hermite quadrature rule (9) and the pseudo-adaptive Gauss-Hermite quadrature rule (13) require evaluating fYi,Ci∣bi,Ψm at their prespecified abscissas across all subjects (see Equations (10) and (14)). Hence, calculating fYi,Ci∣r~t,Ψm requires evaluation of Λ0kmTi across all subjects for each k. We observe from equation (7.4) that for each k, Λ0kmtk1,⋯,Λ0kmtkqk have already been calculated from the mth EM iteration, where tk1>⋯>tkqk are qk distinct observed type k event times. For each i, calculating Λ0kmTi would involve On operations if a global search is performed to find the interval of two adjacent type k event times containing Ti. Consequently, calculating Λ0kmT1,⋯,Λ0kmTn would require On2 operations. However, by taking advantage of the fact that Λ0kmt is a right-continuous and nondecreasing step function, one can obtain Λ0kmT1,⋯,Λ0kmTn from Λ0kmtk1,⋯,Λ0kmtkqk in On operations using the following linear scan algorithm. First, sort the observation times Ti, i=1,⋯,n, in descending order. Denote by i the ranked index of a subject. Then, define a mapping
(15)Λ0kmtk1,Λ0kmtk2,⋯Λ0kmtkqk↦Λ0kmT1,Λ0kmT2,⋯,Λ0kmTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti are scanned forward to calculate Λ0kmTi as follows:
(16)Λ0kmTi=Λ0kmtk1,ifTi≥tk1,Λ0kmtkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,Ti<tkqk.Specifically, start withtk1 and scan through the first set of observation times T1≥⋯≥Tik1 where Tik1≥tk1>Tik1+1, and the corresponding Λ0kTi’s take the value Λ0ktk1. Next, move forward to tk2 and scan through the second set Tik1+1≥⋯≥Tik2, where Tik2≥tk2>Tik2+1, and the corresponding Λ0kTi’s take the value evaluated at Λ0ktk2. Repeat the same process until Tn is scanned. Because the scanned Ti’s for different tkj’s do not overlap, the entire algorithm costs only On operations.


#### 2.2.2. Linear Risk Set Scan for the M-Step

In the M-step, multiple quantities in equations (A.4)–(A.8), such as the cumulative baseline hazard functions and the Hessian matrix and score vector forγk and νk (k=1,2,⋯,K), involve aggregating information over the risk set Rtkj=r:Tr≥tkj at each uncensored event time tkj. These quantities are further aggregated across all tkj’s. If all subjects are scanned to determine the risk set Rtkj at each tkj, then aggregating information over the risk set for all uncensored event times would obviously require On2 operations. Below we explain how to linearize the number of operations for risk set calculations across all uncensored event times by taking advantage of the fact that the risk set is decreasing over time for right censored data.First, to calculateΛ0km+1tkj, j=1,⋯,qk, one needs to compute ∑r∈RtkjexpXr2TγkmEexpνkmbr, j=1,⋯,qk. Because the distinct uncensored event times tk1>⋯>tkqk are arranged in decreasing order, the risk set Rtkj+1 can be decomposed into two disjoint sets: Rtkj+1=Rtkj∪r:Tr∈tkj+1,tkj, and consequently,
(17)∑r∈Rtkj+1ar=∑r∈Rtkjar+∑r:Tr∈tkj+1,tkjar,for any sequence of real numbers a1,⋯,an. It follows from the recursive formula (17) and the fact that the subjects in Rtkj do not need to be scanned to calculate the second term of (17), one can calculate ∑r∈Rtkjar, j=1,⋯qk, in On operations when Tr’s are scanned backward in time.Next, to calculate the Hessian matrixIγkm for γk in (A.5), we first rewrite it as
(18)Iγkm=∑j=1qkΔΛ0ktkjm+1∑r∈RtkjexpXr2TγkmEexpνkmTbrXr2Xr2T,which allows one to linearize its calculation based on (17) with ar=ΔΛ0ktkjm+1expXr2TγkmEexpνkmTbrXr2Xr2T similar to the linear scan algorithm for Λ0km+1tkj’s.Finally, the above linear risk set scan algorithm can be adapted to calculate the Hessian matrix and score vector forγk and νk in equations (A.6)-(A.8) in On operations in a similar fashion.


#### 2.2.3. Linear Scan Algorithm for Standard Error Estimation

The standard error estimation formula in (8) relies on the observed score vector from the profiled likelihood where the baseline hazard is profiled out. However, for each subject i, two components of the score vector, ∇γkliΩ∧;Y,C and ∇νkliΩ∧;Y,C as given in equations (A.12) and (A.13), involve aggregating information either over r∈RTi or over both r∈Rtkj and j:tkj≤Ti. If implemented naively, the aggregation can take either On or On2 operations, respectively. As a result, the observed information matrix can take On3 operation as it requires summing up the information across all subjects. Below we describe a sequential linear scan algorithm to reduce the computational complexity from On3 to On.Our algorithm can be easily explained by considering the calculation of the following expression in the second term of∇γkliΩ∧;Y,C in (7.12):
(19)∑j:tkj≤Tidkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2,fori=1,⋯,n.In other words, we need to computeBTi for i=1,⋯,n, where Bt≡∑j:tkj≤tbkj and
(20)bkj=dkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2.Before going further, we recall that the distinct uncensored event timestk1>⋯>tkqk are in descending order and that the subjects are sorted so that the observation times Ti’s are in descending order.First of all, because the risk set is decreasing over time for right censored data, it follows from Equation (17) that Btk1,⋯,Btkqk can be computed in On operations as one scans through tk1,⋯,tkqk backward in time. Second, analogous to (15), the following linear scan algorithm can be used to calculate BT1,BT2,⋯,BTn from Btk1,⋯,Btkqk:
(21)Btk1,⋯,Btkqk↦BT1,BT2,⋯,BTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti’s are scanned forward to calculate BTi’s as follows:
(22)BTi=Btk1,ifTi≥tk1,Btkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,otherwise.The details are essentially the same as those discussed following Equation (15) and thus omitted here.


## 2.1. Notations and Preliminaries

LetYit be the longitudinal outcome at time t for subject i, i=1,2,⋯,n, and n is the total number of subjects. Suppose the longitudinal outcome Yit is observed at time points tij, j=1,2,⋯,ni, and denote Yi=Yi1,⋯,Yini. Let Ci=Ti,Di be the competing risk data on subject i, where Ti is the observed time to either failure or censoring and Di takes value in 0,1,2,⋯,K, with Di=0 indicating a censored event and Di=k implying the kth type of failure is observed on subject i, k=1,2,⋯,K. The censoring mechanism is assumed to be independent of the failure time.


### 2.1.1. Model

Consider the following joint model in which the longitudinal outcomeYit is characterized by a linear mixed effects model:
(1)Yit=mit+εit,=Xi1tTβ+X~i1tTbi+εit,and the competing risk survival outcome is modeled by a proportional cause-specific hazard model:
(2)λkt∣Xi2,bi,γk,νk=limh⟶0Pt≤Ti<t+h,Di=k∣Ti≥t,Xi2,bih=λ0ktexpXi2Tγk+νkTbi,k=1,⋯,K.In the longitudinal submodel (1), mit is the mean of the longitudinal outcome at time t, Xi1t and X~i1t are column vectors of possibly time-varying covariates associated with the longitudinal outcome Yit, β represents a p×1 vector of fixed effects of Xi1t, bi~Nq0,Σ denotes a q×1 vector of random effects for X~i1t, εit~N0,σ2 is the measurement error independent of bi, and εit1 is independent of εit2 for any t1≠t2. In the competing risk survival submodel (2), λkt∣Xi2,bi,γk,νk is the conditional cause-specific hazard rate for type k failure at time t, given time-invariant covariates Xi2 and the shared random effects bi, and λ0kt is a completely unspecified baseline cause-specific hazard function for type k failure. The two submodels are linked together by the shared random effects bi, and the strength of association is quantified by the association parameters ν1,⋯,νK.


### 2.1.2. Semiparametric Maximum Likelihood Estimation via an EM Algorithm

LetΨ=β,σ2,γ,ν,Σ,λ01·,⋯,λ0K· denote all unknown parameters for the joint models (1) and (2), where γ=γ1T,⋯,γKTT and ν=ν1T,⋯,νKTT. Let IDi=k be the competing event indicator for type k failure, which takes the value 1 if the condition Di=k is satisfied and 0 otherwise. The observed-data likelihood for Ψ is then given by
(3)LΨ;Y,C∝∏i=1nfYi,Ci∣Ψ=∏i=1n∫bifYi∣Ci,bi,ΨfCi∣bi,Ψfbi∣Ψdbi∝∏i=1n∫bi∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1Kλkt∣Xi2,bi,γk,νkIDi=kexp−∫0Ti∑k=1Kλkt∣Xi2,bi,γk,νkdt×12πqΣexp−12biTΣ−1bidbi,which follows from the assumption that Yi and Ci are independent conditional on the covariates and the random effects.BecauseΨ involves K unknown hazard functions, finding its maximum likelihood estimate by maximizing the above observed-data likelihood is nontrivial. However, a customized EM algorithm can be derived to compute the maximum likelihood estimate of Ψ by regarding the latent random effects bi as missing data [3]. The complete-data likelihood based on Y,C,b is
(4)LΨ;Y,C,b∝∏i=1n∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1KΔΛ0kTiexpXi2Tγk+νkTbiIDi=k×exp−∑k=1KΛ0kTiexpXi2Tγk+νkTbi×12πqΣexp−12biTΣ−1bi,where Λ0k. is the cumulative baseline hazard function for type k failure and ΔΛ0kTi=Λ0kTi−Λ0kTi−.The EM algorithm iterates between an expectation step (E-step):(5)QΨ;Ψm≡Eb∣Y,C,ΨmmlogLΨ;Y,C,b,and a maximization step (M-step):
(6)Ψm+1=argmaxΨQΨ;Ψm,until the algorithm converges, where Ψm is the estimate of Ψ from the mth iteration. The E-step (5) involves calculating the following expected values across all subjects: Embi, EmbibiT, Emexpνkbi, Embiexpνkbi, and EmbibiTexpνkbi, where
(7)Emhbi=∫hbifbi∣Yi,Ci,Ψmdbi,for any function h·. Furthermore, it can be shown that the M-step (6) has closed-form solutions for the parameters β, σ2, Σ, and Λ0kt, and that the other parameters γ and ν can be updated using the one-step Newton-Raphson method. Details are provided in equations (7.1)–(7.8) of the supplementary materials (available here).


### 2.1.3. Standard Error Estimation

As discussed in Elashoff et al. [22] (Section 4.1, p.72), standard errors of the parametric components of the semiparametric maximum likelihood estimate can be estimated by profiled likelihood, observed information matrix, or bootstrap. All three methods can be computationally intensive when n is large. Here, we focus on the profiled likelihood-based method and show that its computation can be linearized with respect to n.LetΩ=β,Σ,σ2,γ1,⋯,γK,ν1,⋯,νK denote the parametric component of Ψ and Ω∧ its maximum likelihood estimate. The variance-covariance matrix of Ω∧ can be estimated by inverting the following approximate empirical Fisher information [23–25]:
(8)∑i=1n∇ΩliΩ∧;Y,C∇ΩliΩ∧;Y,CT,where ∇ΩliΩ;Y,C is the observed score vector from the profiled likelihood liΩ;Y,C of Ω on the ith subject by profiling out the baseline hazards. Details of the observed score vector for each parametric component are provided in equations (7.9)–(7.13) of the supplementary materials.


## 2.1.1. Model

Consider the following joint model in which the longitudinal outcomeYit is characterized by a linear mixed effects model:
(1)Yit=mit+εit,=Xi1tTβ+X~i1tTbi+εit,and the competing risk survival outcome is modeled by a proportional cause-specific hazard model:
(2)λkt∣Xi2,bi,γk,νk=limh⟶0Pt≤Ti<t+h,Di=k∣Ti≥t,Xi2,bih=λ0ktexpXi2Tγk+νkTbi,k=1,⋯,K.In the longitudinal submodel (1), mit is the mean of the longitudinal outcome at time t, Xi1t and X~i1t are column vectors of possibly time-varying covariates associated with the longitudinal outcome Yit, β represents a p×1 vector of fixed effects of Xi1t, bi~Nq0,Σ denotes a q×1 vector of random effects for X~i1t, εit~N0,σ2 is the measurement error independent of bi, and εit1 is independent of εit2 for any t1≠t2. In the competing risk survival submodel (2), λkt∣Xi2,bi,γk,νk is the conditional cause-specific hazard rate for type k failure at time t, given time-invariant covariates Xi2 and the shared random effects bi, and λ0kt is a completely unspecified baseline cause-specific hazard function for type k failure. The two submodels are linked together by the shared random effects bi, and the strength of association is quantified by the association parameters ν1,⋯,νK.


## 2.1.2. Semiparametric Maximum Likelihood Estimation via an EM Algorithm

LetΨ=β,σ2,γ,ν,Σ,λ01·,⋯,λ0K· denote all unknown parameters for the joint models (1) and (2), where γ=γ1T,⋯,γKTT and ν=ν1T,⋯,νKTT. Let IDi=k be the competing event indicator for type k failure, which takes the value 1 if the condition Di=k is satisfied and 0 otherwise. The observed-data likelihood for Ψ is then given by
(3)LΨ;Y,C∝∏i=1nfYi,Ci∣Ψ=∏i=1n∫bifYi∣Ci,bi,ΨfCi∣bi,Ψfbi∣Ψdbi∝∏i=1n∫bi∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1Kλkt∣Xi2,bi,γk,νkIDi=kexp−∫0Ti∑k=1Kλkt∣Xi2,bi,γk,νkdt×12πqΣexp−12biTΣ−1bidbi,which follows from the assumption that Yi and Ci are independent conditional on the covariates and the random effects.BecauseΨ involves K unknown hazard functions, finding its maximum likelihood estimate by maximizing the above observed-data likelihood is nontrivial. However, a customized EM algorithm can be derived to compute the maximum likelihood estimate of Ψ by regarding the latent random effects bi as missing data [3]. The complete-data likelihood based on Y,C,b is
(4)LΨ;Y,C,b∝∏i=1n∏j=1ni12πσ2exp−12σ2Yij−Xi1tijTβ−X~i1tijTbi2×∏k=1KΔΛ0kTiexpXi2Tγk+νkTbiIDi=k×exp−∑k=1KΛ0kTiexpXi2Tγk+νkTbi×12πqΣexp−12biTΣ−1bi,where Λ0k. is the cumulative baseline hazard function for type k failure and ΔΛ0kTi=Λ0kTi−Λ0kTi−.The EM algorithm iterates between an expectation step (E-step):(5)QΨ;Ψm≡Eb∣Y,C,ΨmmlogLΨ;Y,C,b,and a maximization step (M-step):
(6)Ψm+1=argmaxΨQΨ;Ψm,until the algorithm converges, where Ψm is the estimate of Ψ from the mth iteration. The E-step (5) involves calculating the following expected values across all subjects: Embi, EmbibiT, Emexpνkbi, Embiexpνkbi, and EmbibiTexpνkbi, where
(7)Emhbi=∫hbifbi∣Yi,Ci,Ψmdbi,for any function h·. Furthermore, it can be shown that the M-step (6) has closed-form solutions for the parameters β, σ2, Σ, and Λ0kt, and that the other parameters γ and ν can be updated using the one-step Newton-Raphson method. Details are provided in equations (7.1)–(7.8) of the supplementary materials (available here).


## 2.1.3. Standard Error Estimation

As discussed in Elashoff et al. [22] (Section 4.1, p.72), standard errors of the parametric components of the semiparametric maximum likelihood estimate can be estimated by profiled likelihood, observed information matrix, or bootstrap. All three methods can be computationally intensive when n is large. Here, we focus on the profiled likelihood-based method and show that its computation can be linearized with respect to n.LetΩ=β,Σ,σ2,γ1,⋯,γK,ν1,⋯,νK denote the parametric component of Ψ and Ω∧ its maximum likelihood estimate. The variance-covariance matrix of Ω∧ can be estimated by inverting the following approximate empirical Fisher information [23–25]:
(8)∑i=1n∇ΩliΩ∧;Y,C∇ΩliΩ∧;Y,CT,where ∇ΩliΩ;Y,C is the observed score vector from the profiled likelihood liΩ;Y,C of Ω on the ith subject by profiling out the baseline hazards. Details of the observed score vector for each parametric component are provided in equations (7.9)–(7.13) of the supplementary materials.


## 2.2. Efficient Algorithms and Implementation of the EM Algorithm and Standard Error Estimation

With naive implementation, multiple quantities in the above E-step, M-step, and standard error estimation will involveOn2 or On3 operations, which become computationally prohibitive at large sample size n. Below we identify these bottlenecks and discuss appropriate linear scan algorithms to reduce their computational complexities to On.


### 2.2.1. Efficient Implementation of the E-Step

At each EM iteration, the E-step (5) requires calculating integral (7) across all subjects. Below we discuss how to accelerate two commonly used numerical integration methods for evaluating these integrals.(1) Standard Gauss-Hermite Quadrature Rule for Numerical Integration. A commonly used method for numerical evaluation of integral (7) is based on the standard Gauss-Hermite quadrature rule [26]:
(9)Emhbi=∫hbifYi,Ci,bi∣ΨmdbifYi,Ci∣Ψm=∫hbifYi,Ci∣bi,Ψmfbi∣Ψmdbi∫fYi,Ci∣bi,Ψmfbi∣Ψmdbi≈∑t1,t2,⋯,tqπthb~tmfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2,where
(10)fYi,Ci∣b~tm,Ψm=fYi∣b~tm,Ψym×fCi∣b~tm,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTb~tm2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkTmb~tmIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkTmb~tm,with Λ0km. the right-continuous and nondecreasing cumulative baseline hazard function for type k failure at the mth EM iteration as defined in Section 7.1 (equation (7.4)) of the supplementary material and ΔΛ0kmTi the jump size of Λ0km. at Ti, q the dimension of the random effects vector, ∑t1,t2,⋯,tq the shorthand for ∑t1=1nq⋯∑tq=1nq, nq the number of quadrature points, ct=ct1,ct2,⋯,ctqT the abscissas with corresponding weights πt, b~tm=2Σ∧m1/2ct the rescaled alternative abscissas, and Σ∧m1/2 the square root of Σ∧m [3]. However, this method is computationally intensive due to multiple factors. First, it usually requires many quadrature points to approximate an integral with sufficient accuracy because the mode of the integrand is often located in a region different from zero. Second, the computational cost increases exponentially with q because the Cartesian product of the abscissas is used to evaluate the integrand with respect to each random effect. Lastly, the alternative abscissas b~tm need to be recalculated at every EM iteration.(2) Pseudo-Adaptive Gauss-Hermite Quadrature Rule for Numerical Integration. When the number of longitudinal measurements per subject is relatively large, Rizopoulos [27] introduced a pseudo-adaptive Gauss-Hermite quadrature rule for numerical approximation of integral (7), which achieves good approximation accuracy with only a small number (nq) of quadrature points and is thus computationally more efficient. The pseudo-adaptive Gauss-Hermite quadrature rule proceeds as follows. First, fit the linear mixed effects model (1) to extract the empirical Bayes estimates of the random effects and its covariance matrix H~i−1:
(11)b~i=Σ^iX~i1V∧i−1Yi−Xi1Tβ^,(12)H~i−1=Σ^−Σ^iX~i1V∧i−1−V∧i−1Xi1T∑i=1nXi1V∧i−1Xi1T−1Xi1V∧i−1×X~i1TΣ^,where Xi1=Xi1ti1,⋯,Xi1tini, X~i1=X~i1ti1,⋯,X~i1tini, and Vi=i1TΣi1+σ2I. Then, define the alternative abscissas r~t=b~i+2H~i−1/2ct and approximate Emhbi by
(13)Emhbi≈∑t1,t2,⋯,tqπthr~tfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2,where
(14)fYi,Ci∣r~t,Ψm=fYi∣r~t,Ψym×fCi∣r~t,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTr~t2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkmTr~tIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkmTr~t,with the notations defined similarly to Equation (10).A derivation of (13) is provided in Section 7.2 of the supplementary materials. The pseudo-adaptive Gauss-Hermite quadrature rule is computationally appealing because the alternate rescaled quadrature points r~t are computed only once before the EM algorithm and do not need to be updated in the EM algorithm. Additionally, the pseudo-adaptive Gauss-Hermite quadrature rule requires fewer quadrature points than the standard Gauss-Hermite quadrature rule to achieve the same numerical approximation accuracy [27]. For example, our simulation results in the supplementary materials (Section A.4, Table A.1) illustrate that the pseudo-adaptive Gauss-Hermite quadrature rule with nq=6 quadrature points produces almost identical results to the standard Gauss-Hermite quadrature rule with nq=20 quadrature points.Remark 1 (Linear calculation ofH~i−1’s across all subjects).
At the first sight, calculatingH~i−1 (12) across all subjects would involve On2 operations since each H~i−1 involves a summation over n subjects. However, because the same quantity ∑i=1nXi1V^i−1Xi1T appears in every H~i−1, one can precompute this quantity and then use the cached value to calculate H~i−1 across all subjects. This way, one can compute the H~i−1’s across all subjects in On operations. Our simulation study in the supplementary material (Section A.5, Figure A.1) shows that applying this simple linearization algorithm can yield a speed-up by a factor of 10 to 10,000 when n grows from 10 to 105. Our implementation is also significantly faster than a popular R package lme4 [28] with over a speed-up by a factor of 10 to 500 as n grows from 10 to 105.(3) Linear Scan Algorithm for CalculatingfYi,Ci∣r~t,Ψm across All subjects. Both the standard Gauss-Hermite quadrature rule (9) and the pseudo-adaptive Gauss-Hermite quadrature rule (13) require evaluating fYi,Ci∣bi,Ψm at their prespecified abscissas across all subjects (see Equations (10) and (14)). Hence, calculating fYi,Ci∣r~t,Ψm requires evaluation of Λ0kmTi across all subjects for each k. We observe from equation (7.4) that for each k, Λ0kmtk1,⋯,Λ0kmtkqk have already been calculated from the mth EM iteration, where tk1>⋯>tkqk are qk distinct observed type k event times. For each i, calculating Λ0kmTi would involve On operations if a global search is performed to find the interval of two adjacent type k event times containing Ti. Consequently, calculating Λ0kmT1,⋯,Λ0kmTn would require On2 operations. However, by taking advantage of the fact that Λ0kmt is a right-continuous and nondecreasing step function, one can obtain Λ0kmT1,⋯,Λ0kmTn from Λ0kmtk1,⋯,Λ0kmtkqk in On operations using the following linear scan algorithm. First, sort the observation times Ti, i=1,⋯,n, in descending order. Denote by i the ranked index of a subject. Then, define a mapping
(15)Λ0kmtk1,Λ0kmtk2,⋯Λ0kmtkqk↦Λ0kmT1,Λ0kmT2,⋯,Λ0kmTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti are scanned forward to calculate Λ0kmTi as follows:
(16)Λ0kmTi=Λ0kmtk1,ifTi≥tk1,Λ0kmtkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,Ti<tkqk.Specifically, start withtk1 and scan through the first set of observation times T1≥⋯≥Tik1 where Tik1≥tk1>Tik1+1, and the corresponding Λ0kTi’s take the value Λ0ktk1. Next, move forward to tk2 and scan through the second set Tik1+1≥⋯≥Tik2, where Tik2≥tk2>Tik2+1, and the corresponding Λ0kTi’s take the value evaluated at Λ0ktk2. Repeat the same process until Tn is scanned. Because the scanned Ti’s for different tkj’s do not overlap, the entire algorithm costs only On operations.


### 2.2.2. Linear Risk Set Scan for the M-Step

In the M-step, multiple quantities in equations (A.4)–(A.8), such as the cumulative baseline hazard functions and the Hessian matrix and score vector forγk and νk (k=1,2,⋯,K), involve aggregating information over the risk set Rtkj=r:Tr≥tkj at each uncensored event time tkj. These quantities are further aggregated across all tkj’s. If all subjects are scanned to determine the risk set Rtkj at each tkj, then aggregating information over the risk set for all uncensored event times would obviously require On2 operations. Below we explain how to linearize the number of operations for risk set calculations across all uncensored event times by taking advantage of the fact that the risk set is decreasing over time for right censored data.First, to calculateΛ0km+1tkj, j=1,⋯,qk, one needs to compute ∑r∈RtkjexpXr2TγkmEexpνkmbr, j=1,⋯,qk. Because the distinct uncensored event times tk1>⋯>tkqk are arranged in decreasing order, the risk set Rtkj+1 can be decomposed into two disjoint sets: Rtkj+1=Rtkj∪r:Tr∈tkj+1,tkj, and consequently,
(17)∑r∈Rtkj+1ar=∑r∈Rtkjar+∑r:Tr∈tkj+1,tkjar,for any sequence of real numbers a1,⋯,an. It follows from the recursive formula (17) and the fact that the subjects in Rtkj do not need to be scanned to calculate the second term of (17), one can calculate ∑r∈Rtkjar, j=1,⋯qk, in On operations when Tr’s are scanned backward in time.Next, to calculate the Hessian matrixIγkm for γk in (A.5), we first rewrite it as
(18)Iγkm=∑j=1qkΔΛ0ktkjm+1∑r∈RtkjexpXr2TγkmEexpνkmTbrXr2Xr2T,which allows one to linearize its calculation based on (17) with ar=ΔΛ0ktkjm+1expXr2TγkmEexpνkmTbrXr2Xr2T similar to the linear scan algorithm for Λ0km+1tkj’s.Finally, the above linear risk set scan algorithm can be adapted to calculate the Hessian matrix and score vector forγk and νk in equations (A.6)-(A.8) in On operations in a similar fashion.


### 2.2.3. Linear Scan Algorithm for Standard Error Estimation

The standard error estimation formula in (8) relies on the observed score vector from the profiled likelihood where the baseline hazard is profiled out. However, for each subject i, two components of the score vector, ∇γkliΩ∧;Y,C and ∇νkliΩ∧;Y,C as given in equations (A.12) and (A.13), involve aggregating information either over r∈RTi or over both r∈Rtkj and j:tkj≤Ti. If implemented naively, the aggregation can take either On or On2 operations, respectively. As a result, the observed information matrix can take On3 operation as it requires summing up the information across all subjects. Below we describe a sequential linear scan algorithm to reduce the computational complexity from On3 to On.Our algorithm can be easily explained by considering the calculation of the following expression in the second term of∇γkliΩ∧;Y,C in (7.12):
(19)∑j:tkj≤Tidkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2,fori=1,⋯,n.In other words, we need to computeBTi for i=1,⋯,n, where Bt≡∑j:tkj≤tbkj and
(20)bkj=dkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2.Before going further, we recall that the distinct uncensored event timestk1>⋯>tkqk are in descending order and that the subjects are sorted so that the observation times Ti’s are in descending order.First of all, because the risk set is decreasing over time for right censored data, it follows from Equation (17) that Btk1,⋯,Btkqk can be computed in On operations as one scans through tk1,⋯,tkqk backward in time. Second, analogous to (15), the following linear scan algorithm can be used to calculate BT1,BT2,⋯,BTn from Btk1,⋯,Btkqk:
(21)Btk1,⋯,Btkqk↦BT1,BT2,⋯,BTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti’s are scanned forward to calculate BTi’s as follows:
(22)BTi=Btk1,ifTi≥tk1,Btkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,otherwise.The details are essentially the same as those discussed following Equation (15) and thus omitted here.


## 2.2.1. Efficient Implementation of the E-Step

At each EM iteration, the E-step (5) requires calculating integral (7) across all subjects. Below we discuss how to accelerate two commonly used numerical integration methods for evaluating these integrals.(1) Standard Gauss-Hermite Quadrature Rule for Numerical Integration. A commonly used method for numerical evaluation of integral (7) is based on the standard Gauss-Hermite quadrature rule [26]:
(9)Emhbi=∫hbifYi,Ci,bi∣ΨmdbifYi,Ci∣Ψm=∫hbifYi,Ci∣bi,Ψmfbi∣Ψmdbi∫fYi,Ci∣bi,Ψmfbi∣Ψmdbi≈∑t1,t2,⋯,tqπthb~tmfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣b~tm,Ψmfb~tm∣Ψmexpct2,where
(10)fYi,Ci∣b~tm,Ψm=fYi∣b~tm,Ψym×fCi∣b~tm,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTb~tm2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkTmb~tmIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkTmb~tm,with Λ0km. the right-continuous and nondecreasing cumulative baseline hazard function for type k failure at the mth EM iteration as defined in Section 7.1 (equation (7.4)) of the supplementary material and ΔΛ0kmTi the jump size of Λ0km. at Ti, q the dimension of the random effects vector, ∑t1,t2,⋯,tq the shorthand for ∑t1=1nq⋯∑tq=1nq, nq the number of quadrature points, ct=ct1,ct2,⋯,ctqT the abscissas with corresponding weights πt, b~tm=2Σ∧m1/2ct the rescaled alternative abscissas, and Σ∧m1/2 the square root of Σ∧m [3]. However, this method is computationally intensive due to multiple factors. First, it usually requires many quadrature points to approximate an integral with sufficient accuracy because the mode of the integrand is often located in a region different from zero. Second, the computational cost increases exponentially with q because the Cartesian product of the abscissas is used to evaluate the integrand with respect to each random effect. Lastly, the alternative abscissas b~tm need to be recalculated at every EM iteration.(2) Pseudo-Adaptive Gauss-Hermite Quadrature Rule for Numerical Integration. When the number of longitudinal measurements per subject is relatively large, Rizopoulos [27] introduced a pseudo-adaptive Gauss-Hermite quadrature rule for numerical approximation of integral (7), which achieves good approximation accuracy with only a small number (nq) of quadrature points and is thus computationally more efficient. The pseudo-adaptive Gauss-Hermite quadrature rule proceeds as follows. First, fit the linear mixed effects model (1) to extract the empirical Bayes estimates of the random effects and its covariance matrix H~i−1:
(11)b~i=Σ^iX~i1V∧i−1Yi−Xi1Tβ^,(12)H~i−1=Σ^−Σ^iX~i1V∧i−1−V∧i−1Xi1T∑i=1nXi1V∧i−1Xi1T−1Xi1V∧i−1×X~i1TΣ^,where Xi1=Xi1ti1,⋯,Xi1tini, X~i1=X~i1ti1,⋯,X~i1tini, and Vi=i1TΣi1+σ2I. Then, define the alternative abscissas r~t=b~i+2H~i−1/2ct and approximate Emhbi by
(13)Emhbi≈∑t1,t2,⋯,tqπthr~tfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2∑t1,t2,⋯,tqπtfYi,Ci∣r~t,Ψmfr~t∣Ψmexpct2,where
(14)fYi,Ci∣r~t,Ψm=fYi∣r~t,Ψym×fCi∣r~t,Ψcm=∏j=1ni12πσ2mexp−12σ2mYij−Xi1tijTβm−X~i1tijTr~t2×∏k=1KΔΛ0kmTiexpXi2Tγkm+νkmTr~tIDi=k×exp−∑k=1KΛ0kmTiexpXi2Tγkm+νkmTr~t,with the notations defined similarly to Equation (10).A derivation of (13) is provided in Section 7.2 of the supplementary materials. The pseudo-adaptive Gauss-Hermite quadrature rule is computationally appealing because the alternate rescaled quadrature points r~t are computed only once before the EM algorithm and do not need to be updated in the EM algorithm. Additionally, the pseudo-adaptive Gauss-Hermite quadrature rule requires fewer quadrature points than the standard Gauss-Hermite quadrature rule to achieve the same numerical approximation accuracy [27]. For example, our simulation results in the supplementary materials (Section A.4, Table A.1) illustrate that the pseudo-adaptive Gauss-Hermite quadrature rule with nq=6 quadrature points produces almost identical results to the standard Gauss-Hermite quadrature rule with nq=20 quadrature points.Remark 1 (Linear calculation ofH~i−1’s across all subjects).
At the first sight, calculatingH~i−1 (12) across all subjects would involve On2 operations since each H~i−1 involves a summation over n subjects. However, because the same quantity ∑i=1nXi1V^i−1Xi1T appears in every H~i−1, one can precompute this quantity and then use the cached value to calculate H~i−1 across all subjects. This way, one can compute the H~i−1’s across all subjects in On operations. Our simulation study in the supplementary material (Section A.5, Figure A.1) shows that applying this simple linearization algorithm can yield a speed-up by a factor of 10 to 10,000 when n grows from 10 to 105. Our implementation is also significantly faster than a popular R package lme4 [28] with over a speed-up by a factor of 10 to 500 as n grows from 10 to 105.(3) Linear Scan Algorithm for CalculatingfYi,Ci∣r~t,Ψm across All subjects. Both the standard Gauss-Hermite quadrature rule (9) and the pseudo-adaptive Gauss-Hermite quadrature rule (13) require evaluating fYi,Ci∣bi,Ψm at their prespecified abscissas across all subjects (see Equations (10) and (14)). Hence, calculating fYi,Ci∣r~t,Ψm requires evaluation of Λ0kmTi across all subjects for each k. We observe from equation (7.4) that for each k, Λ0kmtk1,⋯,Λ0kmtkqk have already been calculated from the mth EM iteration, where tk1>⋯>tkqk are qk distinct observed type k event times. For each i, calculating Λ0kmTi would involve On operations if a global search is performed to find the interval of two adjacent type k event times containing Ti. Consequently, calculating Λ0kmT1,⋯,Λ0kmTn would require On2 operations. However, by taking advantage of the fact that Λ0kmt is a right-continuous and nondecreasing step function, one can obtain Λ0kmT1,⋯,Λ0kmTn from Λ0kmtk1,⋯,Λ0kmtkqk in On operations using the following linear scan algorithm. First, sort the observation times Ti, i=1,⋯,n, in descending order. Denote by i the ranked index of a subject. Then, define a mapping
(15)Λ0kmtk1,Λ0kmtk2,⋯Λ0kmtkqk↦Λ0kmT1,Λ0kmT2,⋯,Λ0kmTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti are scanned forward to calculate Λ0kmTi as follows:
(16)Λ0kmTi=Λ0kmtk1,ifTi≥tk1,Λ0kmtkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,Ti<tkqk.Specifically, start withtk1 and scan through the first set of observation times T1≥⋯≥Tik1 where Tik1≥tk1>Tik1+1, and the corresponding Λ0kTi’s take the value Λ0ktk1. Next, move forward to tk2 and scan through the second set Tik1+1≥⋯≥Tik2, where Tik2≥tk2>Tik2+1, and the corresponding Λ0kTi’s take the value evaluated at Λ0ktk2. Repeat the same process until Tn is scanned. Because the scanned Ti’s for different tkj’s do not overlap, the entire algorithm costs only On operations.


## 2.2.2. Linear Risk Set Scan for the M-Step

In the M-step, multiple quantities in equations (A.4)–(A.8), such as the cumulative baseline hazard functions and the Hessian matrix and score vector forγk and νk (k=1,2,⋯,K), involve aggregating information over the risk set Rtkj=r:Tr≥tkj at each uncensored event time tkj. These quantities are further aggregated across all tkj’s. If all subjects are scanned to determine the risk set Rtkj at each tkj, then aggregating information over the risk set for all uncensored event times would obviously require On2 operations. Below we explain how to linearize the number of operations for risk set calculations across all uncensored event times by taking advantage of the fact that the risk set is decreasing over time for right censored data.First, to calculateΛ0km+1tkj, j=1,⋯,qk, one needs to compute ∑r∈RtkjexpXr2TγkmEexpνkmbr, j=1,⋯,qk. Because the distinct uncensored event times tk1>⋯>tkqk are arranged in decreasing order, the risk set Rtkj+1 can be decomposed into two disjoint sets: Rtkj+1=Rtkj∪r:Tr∈tkj+1,tkj, and consequently,
(17)∑r∈Rtkj+1ar=∑r∈Rtkjar+∑r:Tr∈tkj+1,tkjar,for any sequence of real numbers a1,⋯,an. It follows from the recursive formula (17) and the fact that the subjects in Rtkj do not need to be scanned to calculate the second term of (17), one can calculate ∑r∈Rtkjar, j=1,⋯qk, in On operations when Tr’s are scanned backward in time.Next, to calculate the Hessian matrixIγkm for γk in (A.5), we first rewrite it as
(18)Iγkm=∑j=1qkΔΛ0ktkjm+1∑r∈RtkjexpXr2TγkmEexpνkmTbrXr2Xr2T,which allows one to linearize its calculation based on (17) with ar=ΔΛ0ktkjm+1expXr2TγkmEexpνkmTbrXr2Xr2T similar to the linear scan algorithm for Λ0km+1tkj’s.Finally, the above linear risk set scan algorithm can be adapted to calculate the Hessian matrix and score vector forγk and νk in equations (A.6)-(A.8) in On operations in a similar fashion.


## 2.2.3. Linear Scan Algorithm for Standard Error Estimation

The standard error estimation formula in (8) relies on the observed score vector from the profiled likelihood where the baseline hazard is profiled out. However, for each subject i, two components of the score vector, ∇γkliΩ∧;Y,C and ∇νkliΩ∧;Y,C as given in equations (A.12) and (A.13), involve aggregating information either over r∈RTi or over both r∈Rtkj and j:tkj≤Ti. If implemented naively, the aggregation can take either On or On2 operations, respectively. As a result, the observed information matrix can take On3 operation as it requires summing up the information across all subjects. Below we describe a sequential linear scan algorithm to reduce the computational complexity from On3 to On.Our algorithm can be easily explained by considering the calculation of the following expression in the second term of∇γkliΩ∧;Y,C in (7.12):
(19)∑j:tkj≤Tidkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2,fori=1,⋯,n.In other words, we need to computeBTi for i=1,⋯,n, where Bt≡∑j:tkj≤tbkj and
(20)bkj=dkj∑r∈RtkjexpγkTXr2EexpνkTbrXr2∑r∈RtkjexpγkTXr2EexpνkTbr2.Before going further, we recall that the distinct uncensored event timestk1>⋯>tkqk are in descending order and that the subjects are sorted so that the observation times Ti’s are in descending order.First of all, because the risk set is decreasing over time for right censored data, it follows from Equation (17) that Btk1,⋯,Btkqk can be computed in On operations as one scans through tk1,⋯,tkqk backward in time. Second, analogous to (15), the following linear scan algorithm can be used to calculate BT1,BT2,⋯,BTn from Btk1,⋯,Btkqk:
(21)Btk1,⋯,Btkqk↦BT1,BT2,⋯,BTn,where tk1,⋯,tkqk are scanned forward from the largest to the smallest, and for each tkj, only a subset of the ranked observation times Ti’s are scanned forward to calculate BTi’s as follows:
(22)BTi=Btk1,ifTi≥tk1,Btkj+1,ifTi∈tkj+1,tkj,forsomej∈1,⋯,qk−1,0,otherwise.The details are essentially the same as those discussed following Equation (15) and thus omitted here.


## 3. Simulation Studies

We present a simulation study to illustrate the computational speed-up rendered by the proposed linear algorithms as the sample sizen grows from 100 to 1,000,000. All simulations were run on a MacBook Pro with 6-Core Intel Core i7 processor (2.6 GHz) and 16 GB RAM running MacOS.We generated longitudinal measurementsYij from
(23)Yij=β0+β1tij+β2X2i+b0i+b1itij+εij,which corresponds to model (1) with Xi1tijT=1,tij,X2i and X~i1tijT=1,tij, and competing risk event times from a proportional cause-specific hazard model
(24)λ1t;X1i,X2i,bi,γ1,ν1=λ01texpγ11X1i+γ12X2i+ν1Tbi,(25)λ2t;X1i,X2i,bi,γ2,ν2=λ02texpγ21X1i+γ22X2i+ν2Tbi,where the two submodels (23)–(25) are linked together through the shared random effects bi=b0i,b1iT. In the above joint model, tij=0,1,⋯ represent scheduled visit times, X1i follows N2,1.0, X2i~Bernoulli0.5 is a binary covariate, the random effects bi=b0i,b1iT follows a N20,Σ distribution with Σ11=0.5, Σ22=0.25, andΣ12=0, the measurement errors εij are iid N0,0.5 and independent of bi, and the baseline hazards λ01t and λ01t are constants 0.05 and 0.1, respectively. We simulated noninformative censoring time Vi following exp (20) and let Ti=minTi1∗,Ti2∗,Vi be the observed time (possibly censored) for subject i. The longitudinal measurements for subject i at tij are assumed missing after Ti.We first compared the runtime between three different implementations of the EM algorithm for fitting the joint models (1) and (2) as described in Section 2.2.
(1)
Method 1. This method is a standard implementation of the EM algorithm using the standard Gauss-Hermite quadrature rule in the E-step (Equation (9) with nq=20) without any linear computation.(2)
Method 2. This is a standard implementation of the EM algorithm using the pseudo-adaptive quadrature rule in the E-step (Equation (13) with nq=6) with the linear calculation of H~i−1’s described in Remark 1 and without any other linear computation.(3)
Method 3. This is a method 2+linear scan for calculating fYi,Ci∣r~t,Ψm’s+linear risk set scan for M-step as described in Section 2.2.The number of quadrature points nq for methods 1 and 2 was determined by first trying different values, {10, 20, 30} for method 1 and {6, 9, 12, 15} for method 2, and then choosing the smallest value for which the estimation results are stabilized and similar between the two implementation methods. For comparison purposes, we have also included the runtime of an established joint model R package joineR, which uses a similar EM algorithm for parameter estimation to fit a semiparametric joint model with a slightly different latent association structure in the competing risk submodel [29]. The results are depicted in Figure 1.Figure 1
Runtime (seconds) comparison between three different implementations of the EM algorithm for fitting the joint models (1) and (2) and the joineR package. The details of methods 1-3 are given in Section 3. joineR is an established R package which fits a similar semiparametric joint model with a slightly different latent association structure in the competing risk submodel [29]. Fold change is calculated as the ratio of runtime between two methods.
(a)
Runtime for four methods(b)
Runtime comparison: method 1 versus method 2(c)
Runtime comparison: method 2 versus method 3(d)
Runtime comparison: joineR versus method 3It is seen from Figure1(a) that the runtime of method 3 increases linearly with the sample size, while the runtime of the other three methods grows exponentially. For moderate sample size, method 2 is computationally more efficient than method 1 because it requires fewer quadrature points for numerical integration. However, its computational advantage diminishes as the sample size increases due to the exponentially increasing computational cost of fYi,Ci∣r~t,Ψm’s and risk set calculation in the M-step. By further linearizing the computation of these key components, method 3 has yielded more than 100-fold speed-up over method 2 when n=105, and the speed-up is expected to increase exponentially as n increases (Figure 1(b)). Furthermore, method 3 has demonstrated more than 30-fold speed-up over joineR when n=104. We also note that joineR failed to run when n=105 due to the overload of memory.We also compared the runtime of two implementations of the standard error estimation, with and without linear scan as described in Section2.2.3, and the bootstrap method employed by the joineR package [29]. The results are shown in Figure 2.Figure 2
Runtime (seconds) comparison between two implementations of standard error estimation for fitting the joint models (1) and (2), linear scan and no linear scan as described in Section 2.2.3, and the bootstrap method employed by the joineR package [29]. Fold change is calculated as the ratio of runtime between two methods.
(a)
Runtime for three methods(b)
Runtime comparison: no linear scan versus linear scan(c)
Runtime comparison: joineR versus linear scanIt is seen from Figure2(a) that the implementation with linear scan easily scales to a million subjects, taking only minutes to finish, while the naive implementation without linear scan grinds to a halt when the sample size is 10,000 or larger. Figure 2(b) shows that linear scan can generate a speed-up by a factor of greater than 100,000 when n≥10,000. Similarly, in comparison with joineR that used 100 bootstrap samples for standard error estimation, our standard error estimation method with linear scan generated a speed-up by a factor of greater than 100,000 when n≥5,000.Finally, Figures1 and 2 in Section 3 have focused on contrasting the computational efficiency of different implementations for parameter estimation and standard error estimation in terms of the runtime. We have also compared their parameter estimates and standard error in Section A.6 of the supplementary materials. As one would expect, our three different implementations (methods 1-3) yielded almost identical estimation results, whereas joineR produced similar estimation results for the longitudinal model, but slightly different results for the competing risk model due to its different latent association structure.


## 4. Real Data Examples

We have developed an R package FastJM to implement the efficient algorithms described in Section2. Below we illustrate the improved computational performance of FastJM in comparison to existing joint model R packages on a lung health study (LHS) data with n=5,887 subjects and a UK Biobank data with n=193,287 participants.


### 4.1. Lung Health Study

The lung health study (LHS) data were collected from a ten-center randomized clinical trial on 5,887 middle-aged smokers with mild to moderate chronic obstructive pulmonary disease (COPD) [30]. Patients were randomized into three arms: usual care, smoking intervention and placebo inhaler (SIP), and smoking intervention and active bronchodilator inhaler (SIA). An important objective of the study was to determine if the intervention program with the combination of intensive smoking cessation counseling and an inhaled anticholinergic bronchodilator can slow down the decline in forced expired volume in 1 s (FEV1) during a 5-year follow-up period. Patients’ FEV1 values were collected annually upon recruitment into the study. FEV1 was chosen as the primary outcome since its trajectory is an indicator of a patient’s natural loss of lung function during the progression of COPD. Since not all patients completed the whole study period, about 9.47% of longitudinal measurements were missing. One of the possible reasons for dropout is that treatment was not effective, and hence, missing longitudinal measurements after dropout are nonignorable.Joint modeling ofFEV1 together with the possible informative dropout time provides an attractive approach to deal with nonignorable missing longitudinal data due to dropout. Based on previous findings, we considered the following covariates when characterizing the trajectory of Y=FEV1: time (year), sex, age, body mass index (BMI), baseline number of cigarettes smoked per day, and the logarithm of two-point methacholine concentration-FEV1 O’Connor slope (logslope) [31]. We also included two interaction terms between treatment indicators SIP and SIA and time, so that the difference in the slope of FEV1 between SIP (or SIA) and usual care can be evaluated by testing if the interactions are zero or not. Specifically, we considered the following linear mixed effects model:
(26)Yij=β0+β1tij+β2Xagei+β3XF10CIGSi+β4Dsexi+β5Xlogslopei+β6XBMIi+β7DSIPi+β8DSIAi+β9DSIPi×tij+β10DSIAi×tij+b0i+b1itij+εij,which corresponds to model (1) with Xi1tijT=1,tij,Xagei,XF10CIGSi,Dsexi,Xlogslopei,XBMIi,DSIPi,DSIAi,DSIPi×tij,DSIAi×tij and X~i1tijT=1,tij. The random error term εij~iidN0,σ2 and the random effects bi=b0i,b1iT are assumed normally distributed with zero mean and a covariance matrix Σ. For the dropout time Ti (possibly censored at the end of the study), we assume the Cox proportional hazard submodel.
(27)λit=λ0texpγ1XBMIi+γ2DSIPi+γ3DSIAi+γ4Xlogslopei+γ5Dsexi+γ5Xagei+Wit,where λ0t denotes the baseline hazard function and Wit is a latent association structure that links the two submodels.Table1 compares the runtime of FastJM and some existing joint model packages including joineR [29], different versions of JM [2], JMbayes [32], and JSM [33] with various specifications of λ0t and Wit.Table 1
Runtime comparison between different R packages for joint modeling of a longitudinal and a single event time on the lung health study data.
PackageSemiparametric joint modelsParametric joint modelsFastJMjoineRJSMaJSMbJMbJMa1JMa2JMBayesaJMBayesb∗Baseline hazardUnspecifiedUnspecifiedUnspecifiedUnspecifiedUnspecifiedWeibullB-splineB-splineB-splineWitνTbiνX~i1tTbiνmitνX~i1tTbiνmitνmitνmitνmitνTbiRuntime0.3 min20.4 min1 h 36 min1 h 51 min∗0.9 min1 min19.8 min43 min∗Failed to produce any result due to convergence issue.Among all the semiparametric models (FastJM, joineR,JSMa, and JSMb), FastJM finished in 0.3 minutes while other methods took 20.4 minutes to 111 minutes. As a matter of fact, the runtime of FastJM was even shorter than those of some parametric joint models (JSMa and JSMb). We also observed that JMbayes based on a Bayesian MCMC framework is considerably slower than its frequentist counterpart JM. Finally, the parameter estimates and inference results for the longitudinal outcome were almost identical between all packages, but slightly different for the survival submodel because of their slightly different latent structure Wit. Detailed analysis results are summarized in Section A.7 of the supplementary materials.


### 4.2. UK Biobank Primary Care (UKB-PC) Study

The UK Biobank (UKB) is a prospective cohort study with deep genetic and phenotypic data collected on approximately 500,000 individuals, aged 37-73 years, from the general population between 2006 and 2010 in the United Kingdom [18, 19]. Participants attended assessment at their closest clinic center where they completed questionnaires, took physical measurements, and provided biological samples (blood, urine, and saliva) as a baseline assessment visit. Hospital admission records were available until February, 2018, for the full UKB cohort, whereas linkage to primary care records was available for 45% of the UKB cohort (approximately 230,000 participants) until May, 2017, for Scotland, September, 2017, for Wales, and August, 2017, for England. The detailed linkage procedures relating to primary care records are available online (https://biobank.ndph.ox.ac.uk/showcase/showcase/docs/primary_care_data.pdf).In this example, we consider a joint model of longitudinal systolic blood pressure (SBP) measurements and a competing risk event time defined as age-at-onset of type 2 diabetes (T2D) as the first risk and age-at-onset of stroke, myocardial infarction (MI), or all-cause death as the second risk, whichever occurred first. Age-at-onset of outcomes were based on participants’ primary care or hospital records, whichever occurred first. Follow-up was censored at the primary care data end date for the relevant country or the date of outcomes, if this occurred earlier. SBP measures were extracted from either baseline assessment visit or primary care data. Covariates include sex, ethnicity, and BMI measured during baseline visit. However, considering the imbalanced racial distribution in this case study, we only considered white vs. non-white ethnicity groups. Specifically, the joint model consists of a linear mixed effects model for the longitudinal outcome (SBP),(28)Yij=β0+β1tij+β2XBMIi+β3Dmalei+β4Dnon−whitei+b0i+b1itij+εij,Model−Lwhich corresponds to model (1) with Xi1tijT=1,tij,XBMIi,Dmalei,Dnon−whitei and X~i1tijT=1,tij and a proportional cause-specific hazard model for the competing risk event outcome
(29)λikt=λk0texpγk1XBMIi+γk2Dmalei+γk3Dnon−whitei+Wikt,Model−PCHfor k=1,2. In Model-L, the random error term εij~iidN0,σ2 and the random effects bi=b0i,b1iT are assumed normally distributed with zero mean and covariance matrix Σ. In Model-PCH, k=1 denotes type 2 diabetes and k=2 stroke, λk0t denotes the baseline cause-specific hazard function for cause k, and Wikt denotes the latent association structure of SBP with cause k risk.To our knowledge, besides our FastJM package, joineR and JM are two other current joint model R packages that are capable of handling competing risk event outcomes. However, because JM encountered convergence issues, we will focus on FastJM and joineR in this case study. Table2 compares the runtime of FastJM and joineR on a subset of 5,000 and 20,000 participants randomly selected from the UKB-PC data and the full UKB-PC data with 193,287 participants.Table 2
Runtime comparison between different R packages for semiparametric joint modeling of longitudinal SBP trajectory and competing risk event time on the UK Biobank primary care (UKB-PC) data.
PackageFastJMjoineRUKB-PC subset (n=5,000)1 min3.3 hUKB-PC subset (n=20,000)4.4 min33 hUKB-PC full data (n=193,287)1 h∗∗Failed to produce any result due to computational failure.Table2 shows that for the UKB-PC subset of 5,000 participants, FastJM finished within 1 minute, while joineR took 3.3 hours to finish. For the UKB-PC subset of 20,000 participants, FastJM finished within 5 minutes, while joineR took 33 hours to run. For the UKB-PC full data with 193,287 participants, FastJM finished within 1 hour, whereas joineR encountered a computational failure.Finally, the analysis results produced by FastJM and joineR are similar for the longitudinal submodel for the UKB-PC subset of 5,000 and 20,000 participants and for UKB-PC full data. For the survival submodel, the analysis results are also similar for most parameters except for the association parameters due to the different latent structureWit between two packages. Detailed analysis results are provided in Section A.8 of the supplementary materials.


## 4.1. Lung Health Study

The lung health study (LHS) data were collected from a ten-center randomized clinical trial on 5,887 middle-aged smokers with mild to moderate chronic obstructive pulmonary disease (COPD) [30]. Patients were randomized into three arms: usual care, smoking intervention and placebo inhaler (SIP), and smoking intervention and active bronchodilator inhaler (SIA). An important objective of the study was to determine if the intervention program with the combination of intensive smoking cessation counseling and an inhaled anticholinergic bronchodilator can slow down the decline in forced expired volume in 1 s (FEV1) during a 5-year follow-up period. Patients’ FEV1 values were collected annually upon recruitment into the study. FEV1 was chosen as the primary outcome since its trajectory is an indicator of a patient’s natural loss of lung function during the progression of COPD. Since not all patients completed the whole study period, about 9.47% of longitudinal measurements were missing. One of the possible reasons for dropout is that treatment was not effective, and hence, missing longitudinal measurements after dropout are nonignorable.Joint modeling ofFEV1 together with the possible informative dropout time provides an attractive approach to deal with nonignorable missing longitudinal data due to dropout. Based on previous findings, we considered the following covariates when characterizing the trajectory of Y=FEV1: time (year), sex, age, body mass index (BMI), baseline number of cigarettes smoked per day, and the logarithm of two-point methacholine concentration-FEV1 O’Connor slope (logslope) [31]. We also included two interaction terms between treatment indicators SIP and SIA and time, so that the difference in the slope of FEV1 between SIP (or SIA) and usual care can be evaluated by testing if the interactions are zero or not. Specifically, we considered the following linear mixed effects model:
(26)Yij=β0+β1tij+β2Xagei+β3XF10CIGSi+β4Dsexi+β5Xlogslopei+β6XBMIi+β7DSIPi+β8DSIAi+β9DSIPi×tij+β10DSIAi×tij+b0i+b1itij+εij,which corresponds to model (1) with Xi1tijT=1,tij,Xagei,XF10CIGSi,Dsexi,Xlogslopei,XBMIi,DSIPi,DSIAi,DSIPi×tij,DSIAi×tij and X~i1tijT=1,tij. The random error term εij~iidN0,σ2 and the random effects bi=b0i,b1iT are assumed normally distributed with zero mean and a covariance matrix Σ. For the dropout time Ti (possibly censored at the end of the study), we assume the Cox proportional hazard submodel.
(27)λit=λ0texpγ1XBMIi+γ2DSIPi+γ3DSIAi+γ4Xlogslopei+γ5Dsexi+γ5Xagei+Wit,where λ0t denotes the baseline hazard function and Wit is a latent association structure that links the two submodels.Table1 compares the runtime of FastJM and some existing joint model packages including joineR [29], different versions of JM [2], JMbayes [32], and JSM [33] with various specifications of λ0t and Wit.Table 1
Runtime comparison between different R packages for joint modeling of a longitudinal and a single event time on the lung health study data.
PackageSemiparametric joint modelsParametric joint modelsFastJMjoineRJSMaJSMbJMbJMa1JMa2JMBayesaJMBayesb∗Baseline hazardUnspecifiedUnspecifiedUnspecifiedUnspecifiedUnspecifiedWeibullB-splineB-splineB-splineWitνTbiνX~i1tTbiνmitνX~i1tTbiνmitνmitνmitνmitνTbiRuntime0.3 min20.4 min1 h 36 min1 h 51 min∗0.9 min1 min19.8 min43 min∗Failed to produce any result due to convergence issue.Among all the semiparametric models (FastJM, joineR,JSMa, and JSMb), FastJM finished in 0.3 minutes while other methods took 20.4 minutes to 111 minutes. As a matter of fact, the runtime of FastJM was even shorter than those of some parametric joint models (JSMa and JSMb). We also observed that JMbayes based on a Bayesian MCMC framework is considerably slower than its frequentist counterpart JM. Finally, the parameter estimates and inference results for the longitudinal outcome were almost identical between all packages, but slightly different for the survival submodel because of their slightly different latent structure Wit. Detailed analysis results are summarized in Section A.7 of the supplementary materials.


## 4.2. UK Biobank Primary Care (UKB-PC) Study

The UK Biobank (UKB) is a prospective cohort study with deep genetic and phenotypic data collected on approximately 500,000 individuals, aged 37-73 years, from the general population between 2006 and 2010 in the United Kingdom [18, 19]. Participants attended assessment at their closest clinic center where they completed questionnaires, took physical measurements, and provided biological samples (blood, urine, and saliva) as a baseline assessment visit. Hospital admission records were available until February, 2018, for the full UKB cohort, whereas linkage to primary care records was available for 45% of the UKB cohort (approximately 230,000 participants) until May, 2017, for Scotland, September, 2017, for Wales, and August, 2017, for England. The detailed linkage procedures relating to primary care records are available online (https://biobank.ndph.ox.ac.uk/showcase/showcase/docs/primary_care_data.pdf).In this example, we consider a joint model of longitudinal systolic blood pressure (SBP) measurements and a competing risk event time defined as age-at-onset of type 2 diabetes (T2D) as the first risk and age-at-onset of stroke, myocardial infarction (MI), or all-cause death as the second risk, whichever occurred first. Age-at-onset of outcomes were based on participants’ primary care or hospital records, whichever occurred first. Follow-up was censored at the primary care data end date for the relevant country or the date of outcomes, if this occurred earlier. SBP measures were extracted from either baseline assessment visit or primary care data. Covariates include sex, ethnicity, and BMI measured during baseline visit. However, considering the imbalanced racial distribution in this case study, we only considered white vs. non-white ethnicity groups. Specifically, the joint model consists of a linear mixed effects model for the longitudinal outcome (SBP),(28)Yij=β0+β1tij+β2XBMIi+β3Dmalei+β4Dnon−whitei+b0i+b1itij+εij,Model−Lwhich corresponds to model (1) with Xi1tijT=1,tij,XBMIi,Dmalei,Dnon−whitei and X~i1tijT=1,tij and a proportional cause-specific hazard model for the competing risk event outcome
(29)λikt=λk0texpγk1XBMIi+γk2Dmalei+γk3Dnon−whitei+Wikt,Model−PCHfor k=1,2. In Model-L, the random error term εij~iidN0,σ2 and the random effects bi=b0i,b1iT are assumed normally distributed with zero mean and covariance matrix Σ. In Model-PCH, k=1 denotes type 2 diabetes and k=2 stroke, λk0t denotes the baseline cause-specific hazard function for cause k, and Wikt denotes the latent association structure of SBP with cause k risk.To our knowledge, besides our FastJM package, joineR and JM are two other current joint model R packages that are capable of handling competing risk event outcomes. However, because JM encountered convergence issues, we will focus on FastJM and joineR in this case study. Table2 compares the runtime of FastJM and joineR on a subset of 5,000 and 20,000 participants randomly selected from the UKB-PC data and the full UKB-PC data with 193,287 participants.Table 2
Runtime comparison between different R packages for semiparametric joint modeling of longitudinal SBP trajectory and competing risk event time on the UK Biobank primary care (UKB-PC) data.
PackageFastJMjoineRUKB-PC subset (n=5,000)1 min3.3 hUKB-PC subset (n=20,000)4.4 min33 hUKB-PC full data (n=193,287)1 h∗∗Failed to produce any result due to computational failure.Table2 shows that for the UKB-PC subset of 5,000 participants, FastJM finished within 1 minute, while joineR took 3.3 hours to finish. For the UKB-PC subset of 20,000 participants, FastJM finished within 5 minutes, while joineR took 33 hours to run. For the UKB-PC full data with 193,287 participants, FastJM finished within 1 hour, whereas joineR encountered a computational failure.Finally, the analysis results produced by FastJM and joineR are similar for the longitudinal submodel for the UKB-PC subset of 5,000 and 20,000 participants and for UKB-PC full data. For the survival submodel, the analysis results are also similar for most parameters except for the association parameters due to the different latent structureWit between two packages. Detailed analysis results are provided in Section A.8 of the supplementary materials.


## 5. Discussion

We have developed customized linear scan algorithms to reduce the computational complexity fromOn2 or On3 to On within each iteration of the EM algorithm and in the standard error estimation step for a semiparametric joint model of a longitudinal biomarker and a competing risk event time outcome. Through simulation and case studies, we have demonstrated that the efficient implementation can generate a speed-up by a factor of up to hundreds of thousands and oftentimes reduce the runtime from days to minutes when the sample size is large (n>104), making it feasible to fit a joint model on a large data in real time.The ideas and techniques of this paper can potentially be adapted to improve computational efficiency for other joint models. For instance, the linear computational algorithm in Remark1 for computing the variance-covariance matrices of empirical Bayes estimates of the random effects is not specific to the joint model considered in this paper and can be used in any procedure that uses the pseudo-adaptive quadrature rule. Also, although we have focused on joint modeling of a single biomarker with a time-to-event outcome, our methodology can be easily extended to handle multiple biomarkers in a similar fashion. It is also important to note that the linear risk set scan algorithm is limited to the share random effect joint model in which the Cox submodel (2) only involves time-independent covariates. If the Cox submodel contains time-dependent covariates such as the present value of the longitudinal marker, then one may have to impose more restrictive assumptions such as assuming a parametric baseline hazard in order to linearize the computation costs with respect to the sample size.This paper has focused on linearizing the computation with respect to the sample size within the framework of classical EM algorithm that is coupled with the pseudo-adaptive quadrature rule for numerical integration in the E-step. It would be interesting to investigate if coupling our algorithms with other numerical integration methods such as quasi-Monte Carlo method [34] in the E-step or with other variations of EM algorithms such as the stochastic EM algorithm (stEM) [35] or Turbo-EM [36] could further enhance the computational efficiency, especially when there are 3 or more random effects in the model. Finally, current joint model implementations are generally not scalable as the number of longitudinal measurement grows large, rendering it infeasible to fit dense longitudinal data such as those generated from modern wearable devices for dynamic prediction of a health outcome. Future research is warranted to develop novel joint modeling procedures that are scalable to large number of subjects, random effects, and longitudinal measurements.


## 6. Software

A user-friendly R package FastJM [37] has been developed to fit the shared parameter joint model using the efficient algorithms developed in this paper and is publicly available on the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=FastJM.



---

*Source: 1362913-2022-02-08.xml*