Joint modelling of multivariate longitudinal and time-to-event data

Presented at Durham University, Department of Mathematical Sciences

Transcript

1. Introduction Model Estimation Software & Example revisited Prediction References Joint modelling of multivariate longitudinal and time-to-event data Graeme L. Hickey1, Pete Philipson2, Andrea Jorgensen1, Ruwanthi-Kolamunnage-Dona1 1Department of Biostatistics, University of Liverpool, UK 2Mathematics, Physics and Electrical Engineering, University of Northumbria, UK graeme.hickey@liverpool.ac.uk Funded by Grant MR/M013227/1 13th November 2017 GL. Hickey Joint modelling of multivariate data 1 / 48
2. Introduction Model Estimation Software & Example revisited Prediction References Joint modelling Covariate data Longitudinal outcomes data Random effects GL. Hickey Joint modelling of multivariate data 2 / 48
3. Introduction Model Estimation Software & Example revisited Prediction References Joint modelling Covariate data Time-to-event outcomes data Frailties GL. Hickey Joint modelling of multivariate data 2 / 48
4. Introduction Model Estimation Software & Example revisited Prediction References Joint modelling Covariate data Time-to-event outcomes data Frailties Longitudinal outcomes data Random effects GL. Hickey Joint modelling of multivariate data 2 / 48
5. Introduction Model Estimation Software & Example revisited Prediction References Joint modelling Covariate data Time-to-event outcomes data Frailties Longitudinal outcomes data Random effects GL. Hickey Joint modelling of multivariate data 2 / 48
6. Introduction Model Estimation Software & Example revisited Prediction References Why use a joint model? Interest lies with adjustment of inferences about longitudinal measurements for possibly outcome-dependent drop-out adjustment of inferences about the time-to-event distribution conditional on intermediate and/or error prone longitudinal measurements the joint evolution of the measurement and event time processes biomarker surrogacy dynamic prediction GL. Hickey Joint modelling of multivariate data 3 / 48
7. Introduction Model Estimation Software & Example revisited Prediction References Motivation for multivariate joint models Clinical studies often repeatedly measure multiple biomarkers or other measurements and an event time Research has predominantly focused on a single event time and single measurement outcome Ignoring correlation leads to bias and reduced efficiency in estimation Harnessing all available information in a single model is advantageous and should lead to improved model predictions GL. Hickey Joint modelling of multivariate data 4 / 48
8. Introduction Model Estimation Software & Example revisited Prediction References Clinical example Figure source: https://www.medgadget.com Primary biliary cirrhosis (PBC) is a chronic liver disease char- acterized by inflammatory de- struction of the small bile ducts, which eventually leads to cirrho- sis of the liver and death GL. Hickey Joint modelling of multivariate data 5 / 48
9. Introduction Model Estimation Software & Example revisited Prediction References Clinical example Consider a subset of 154 patients randomized to placebo treatment from Mayo Clinic trial (Murtaugh et al. 1994) Multiple biomarkers repeatedly measured at intermittent times, of which we consider 3 clinically relevant ones: 1 serum bilirunbin (mg/dl) 2 serum albumin (mg/dl) 3 prothrombin time (seconds) GL. Hickey Joint modelling of multivariate data 6 / 48
10. Introduction Model Estimation Software & Example revisited Prediction References Objective 1 1 Determine if longitudinal biomarker trajectories are associated with death GL. Hickey Joint modelling of multivariate data 7 / 48
11. Introduction Model Estimation Software & Example revisited Prediction References Objective 1 1 Determine if longitudinal biomarker trajectories are associated with death Objective 2 1 Dynamically predict the biomarker trajectories and time to death for a new patient GL. Hickey Joint modelling of multivariate data 7 / 48
12. Introduction Model Estimation Software & Example revisited Prediction References Objective 1 1 Determine if longitudinal biomarker trajectories are associated with death Objective 2 1 Dynamically predict the biomarker trajectories and time to death for a new patient Objective 3 1 Wrap it all up into a freely available software package GL. Hickey Joint modelling of multivariate data 7 / 48
13. Introduction Model Estimation Software & Example revisited Prediction References Prothrombin time (0.1 x seconds)−4 Albumin (mg dL) Serum bilirubin (loge mg dL) 0 5 10 0 5 10 0 5 10 0.0 0.5 1.0 1.5 2 4 6 8 −2 0 2 4 Time from registration (years) Alive Dead GL. Hickey Joint modelling of multivariate data 8 / 48
14. Introduction Model Estimation Software & Example revisited Prediction References 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Time from registration (years) Survival probability nevents = 69, 44.8% GL. Hickey Joint modelling of multivariate data 9 / 48
15. Introduction Model Estimation Software & Example revisited Prediction References Data For each subject i = 1, . . . , n, we observe yi = (yi1 , . . . , yiK ) is a K-variate continuous outcome vector, where each yik denotes an (nik × 1)-vector of observed longitudinal measurements for the k-th outcome type: yik = (yi1k, . . . , yinik k) Observation times tijk for j = 1, . . . , nik, which can differ between subjects and outcomes (Ti , δi ), where Ti = min(T∗ i , Ci ), where T∗ i is the true event time, Ci corresponds to a potential right-censoring time, and δi is the failure indicator equal to 1 if the failure is observed (T∗ i ≤ Ci ) and 0 otherwise GL. Hickey Joint modelling of multivariate data 10 / 48
16. Introduction Model Estimation Software & Example revisited Prediction References Longitudinal sub-model Following Henderson et al. (2000) for the univariate case yi (t) = µi (t) + W1i (t) + εi (t), where εi (t) is the model error term, which is i.i.d. N(0, σ2) and independent of W1i (t) µi (t) = xi (t)β is the mean response xi (t) is a p-vector of (possibly) time-varying covariates with corresponding fixed effect terms β W1i (t) is a zero-mean latent Gaussian process GL. Hickey Joint modelling of multivariate data 11 / 48
17. Introduction Model Estimation Software & Example revisited Prediction References Longitudinal sub-model We can extend it to K-separate sub-models (with k = 1, . . . , K) yik(t) = µik(t) + W (k) 1i (t) + εik(t), where εik(t) is the model error term, which is i.i.d. N(0, σ2 k ) and independent of W (k) 1i (t) µik(t) = xik (t)βk is the mean response xik(t) is a pk-vector of (possibly) time-varying covariates with corresponding fixed effect terms βk W (k) 1i (t) is a zero-mean latent Gaussian process GL. Hickey Joint modelling of multivariate data 11 / 48
18. Introduction Model Estimation Software & Example revisited Prediction References Time-to-event sub-model λi (t) = lim dt→0 P(t ≤ Ti < t + dt | Ti ≥ t) dt = λ0(t) exp vi (t)γv + W2i (t) , where λ0(·) is an unspecified baseline hazard function vi (t) is a q-vector of (possibly) time-varying covariates with corresponding fixed effect terms γv W2i (t) is a zero-mean latent Gaussian process, independent of the censoring process GL. Hickey Joint modelling of multivariate data 12 / 48
19. Introduction Model Estimation Software & Example revisited Prediction References Correlation Following Laird and Ware (1982): W (k) 1i (t) = zik (t)bik for k = 1, . . . , K GL. Hickey Joint modelling of multivariate data 13 / 48
20. Introduction Model Estimation Software & Example revisited Prediction References Correlation Following Laird and Ware (1982): W (k) 1i (t) = zik (t)bik for k = 1, . . . , K Three sources of correlation: GL. Hickey Joint modelling of multivariate data 13 / 48
21. Introduction Model Estimation Software & Example revisited Prediction References Correlation Following Laird and Ware (1982): W (k) 1i (t) = zik (t)bik for k = 1, . . . , K Three sources of correlation: 1 Within-subject correlation between longitudinal measurements: bik ∼ N(0, Dkk) GL. Hickey Joint modelling of multivariate data 13 / 48
22. Introduction Model Estimation Software & Example revisited Prediction References Correlation Following Laird and Ware (1982): W (k) 1i (t) = zik (t)bik for k = 1, . . . , K Three sources of correlation: 1 Within-subject correlation between longitudinal measurements: bik ∼ N(0, Dkk) 2 Between longitudinal outcomes correlation: cov(bik, bil ) = Dkl for k = l GL. Hickey Joint modelling of multivariate data 13 / 48
23. Introduction Model Estimation Software & Example revisited Prediction References Correlation Following Laird and Ware (1982): W (k) 1i (t) = zik (t)bik for k = 1, . . . , K Three sources of correlation: 1 Within-subject correlation between longitudinal measurements: bik ∼ N(0, Dkk) 2 Between longitudinal outcomes correlation: cov(bik, bil ) = Dkl for k = l 3 Correlation between sub-models1: W2i (t) = K k=1 γykW (k) 1i (t) 1Extends model proposed Henderson et al. (2000) GL. Hickey Joint modelling of multivariate data 13 / 48
24. Introduction Model Estimation Software & Example revisited Prediction References Association structure: alternatives Many other proposals for association structures in the literature: Current value parameterisation: W2i (t) = γy {µi (t) + W1i (t)} Random effects parameterisation: W2i (t) = γy1 bi Bivariate distribution: (W1i , W2i ) ∼ N(0, Ω) Random-slopes parameterisation: W2i (t) = γy1 {µi (t) + W1i (t)} + γy2 ∂ ∂t {µi (t) + W1i (t)} . . . GL. Hickey Joint modelling of multivariate data 14 / 48
25. Introduction Model Estimation Software & Example revisited Prediction References Likelihood We can re-write the longitudinal sub-model as yi | bi , β, Σi ∼ N(Xi β + Zi bi , Σi ), with bi | D ∼ N(0, D), where β = (β1 , . . . , βK ), bi = (bi1 , . . . , biK ) , and Xi =      Xi1 · · · 0 . . . ... . . . 0 · · · XiK      , Zi =      Zi1 · · · 0 . . . ... . . . 0 · · · ZiK      , D =      D11 · · · D1K . . . ... . . . D1K · · · DKK      Σi =      σ2 1 Ini1 · · · 0 . . . ... . . . 0 · · · σ2 K IniK      GL. Hickey Joint modelling of multivariate data 15 / 48
26. Introduction Model Estimation Software & Example revisited Prediction References Likelihood The observed data likelihood is given by n i=1 ∞ −∞ f (yi | bi , θ)f (Ti , δi | bi , θ)f (bi | θ)dbi where θ = (β , vech(D), σ2 1 , . . . , σ2 K , λ0(t), γv , γy ) GL. Hickey Joint modelling of multivariate data 16 / 48
27. Introduction Model Estimation Software & Example revisited Prediction References Likelihood The observed data likelihood is given by n i=1 ∞ −∞ f (yi | bi , θ)f (Ti , δi | bi , θ)f (bi | θ)dbi where θ = (β , vech(D), σ2 1 , . . . , σ2 K , λ0(t), γv , γy ), and f (yi | bi , θ) = K k=1 (2π)−nik 2 |Σi |−1 2 exp − 1 2 (yi − Xi β − Zi bi ) Σ−1 i (yi − Xi β − Zi bi ) GL. Hickey Joint modelling of multivariate data 16 / 48
28. Introduction Model Estimation Software & Example revisited Prediction References Likelihood The observed data likelihood is given by n i=1 ∞ −∞ f (yi | bi , θ)f (Ti , δi | bi , θ)f (bi | θ)dbi where θ = (β , vech(D), σ2 1 , . . . , σ2 K , λ0(t), γv , γy ), and f (Ti , δi | bi ; θ) = λ0(Ti ) exp vi γv + W2i (Ti , bi ) δi exp − Ti 0 λ0(u) exp vi γv + W2i (u, bi ) du GL. Hickey Joint modelling of multivariate data 16 / 48
29. Introduction Model Estimation Software & Example revisited Prediction References Likelihood The observed data likelihood is given by n i=1 ∞ −∞ f (yi | bi , θ)f (Ti , δi | bi , θ)f (bi | θ)dbi where θ = (β , vech(D), σ2 1 , . . . , σ2 K , λ0(t), γv , γy ), and f (bi | θ) = (2π)− r 2 |D|−1 2 exp − 1 2 bi D−1bi , with r = dim(bi ) GL. Hickey Joint modelling of multivariate data 16 / 48
30. Introduction Model Estimation Software & Example revisited Prediction References Estimation Multiple approaches have been considered over the years: Markov chain Monte Carlo (MCMC) Direct likelihood maximisation (e.g. Newton-methods) Generalised estimating equations EM algorithm (treating the random effects as missing data) . . . GL. Hickey Joint modelling of multivariate data 17 / 48
31. Introduction Model Estimation Software & Example revisited Prediction References EM algorithm (Dempster et al. 1977) E-step. At the m-th iteration, we compute the expected log-likelihood of the complete data conditional on the observed data and the current estimate of the parameters. Q(θ | ˆ θ(m)) = n i=1 E log f (yi , Ti , δi , bi | θ) , = n i=1 ∞ −∞ log f (yi , Ti , δi , bi | θ) f (bi | Ti , δi , yi ; ˆ θ(m))dbi GL. Hickey Joint modelling of multivariate data 18 / 48
32. Introduction Model Estimation Software & Example revisited Prediction References EM algorithm (Dempster et al. 1977) E-step. At the m-th iteration, we compute the expected log-likelihood of the complete data conditional on the observed data and the current estimate of the parameters. Q(θ | ˆ θ(m)) = n i=1 E log f (yi , Ti , δi , bi | θ) , = n i=1 ∞ −∞ log f (yi , Ti , δi , bi | θ) f (bi | Ti , δi , yi ; ˆ θ(m))dbi M-step. We maximise Q(θ | ˆ θ(m)) with respect to θ. namely, ˆ θ(m+1) = arg max θ Q(θ | ˆ θ(m)) GL. Hickey Joint modelling of multivariate data 18 / 48
33. Introduction Model Estimation Software & Example revisited Prediction References M-step: closed form estimators ˆ λ0(t) = n i=1 δi I(Ti = t) n i=1 E exp vi γv + W2i (t, bi ) I(Ti ≥ t) ˆ β = n i=1 Xi Xi −1 n i=1 Xi (yi − Zi E[bi ]) ˆ σ2 k = 1 n i=1 nik n i=1 (yik − Xikβk) (yik − Xikβk − 2ZikE[bik]) +trace Zik ZikE[bikbik ] ˆ D = 1 n n i=1 E bi bi GL. Hickey Joint modelling of multivariate data 19 / 48
34. Introduction Model Estimation Software & Example revisited Prediction References M-step: non-closed form estimators There is no closed form update for γ = (γv , γy ), so use a one-step Newton-Raphson iteration ˆ γ(m+1) = ˆ γ(m) + I ˆ γ(m) −1 S ˆ γ(m) , where S(γ) = n i=1 δi E [˜ vi (Ti )] − Ti 0 λ0(u)E ˜ vi (u) exp{˜ vi (u)γ} du I(γ) = − ∂ ∂γ S(γ) with ˜ vi (t) = vi , zi1 (t)bi1, . . . , ziK (t)biK a (q + K)–vector GL. Hickey Joint modelling of multivariate data 20 / 48
35. Introduction Model Estimation Software & Example revisited Prediction References MCEM algorithm E-step requires calculating several multidimensional integrals of form E h(bi ) | Ti , δi , yi ; ˆ θ Gauss-quadrature can be slow if dim(bi ) is large ⇒ might not scale well as K increases Instead, we use the Monte Carlo Expectation-Maximization (MCEM; Wei and Tanner 1990) M-step updates remain the same GL. Hickey Joint modelling of multivariate data 21 / 48
36. Introduction Model Estimation Software & Example revisited Prediction References Monte Carlo E-step Conventional EM algorithm: use quadrature to compute E h(bi ) | Ti , δi , yi ; ˆ θ = ∞ −∞ h(bi )f (bi | yi ; ˆ θ)f (Ti , δi | bi ; ˆ θ)dbi ∞ −∞ f (bi | yi ; ˆ θ)f (Ti , δi | bi ; ˆ θ)dbi , where h(·) = any known fuction, bi | yi , θ ∼ N Ai Zi Σ−1 i (yi − Xi β) , Ai , and Ai = Zi Σ−1 i Zi + D−1 −1 GL. Hickey Joint modelling of multivariate data 22 / 48
37. Introduction Model Estimation Software & Example revisited Prediction References Monte Carlo E-step MCEM algorithm E-step: use Monte Carlo integration to compute E h(bi ) | Ti , δi , yi ; ˆ θ ≈ 1 N N d=1 h b(d) i f Ti , δi | b(d) i ; ˆ θ 1 N N d=1 f Ti , δi | b(d) i ; ˆ θ where h(·) = any known fuction, bi | yi , θ ∼ N Ai Zi Σ−1 i (yi − Xi β) , Ai , and Ai = Zi Σ−1 i Zi + D−1 −1 b(1) i , b(2) i , . . . , b(N) i ∼ bi | yi , θ a Monte Carlo draw GL. Hickey Joint modelling of multivariate data 22 / 48
38. Introduction Model Estimation Software & Example revisited Prediction References Speeding up convergence Monte Carlo integration converges at a rate of O(N−1/2), which is independent of K and r = dim(bi ) EM algorithm convergences linearly Can we speed this up? GL. Hickey Joint modelling of multivariate data 23 / 48
39. Introduction Model Estimation Software & Example revisited Prediction References Speeding up convergence Monte Carlo integration converges at a rate of O(N−1/2), which is independent of K and r = dim(bi ) EM algorithm convergences linearly Can we speed this up? 1 Antithetic variates 2 Quasi-Monte Carlo GL. Hickey Joint modelling of multivariate data 23 / 48
40. Introduction Model Estimation Software & Example revisited Prediction References Variance reduction Instead of directly sampling from the MVN distribution for bi | yi ; θ, we apply a variance reduction technique Antithetic simulation Sample Ω ∼ N(0, Ir ) and obtain the pairs Ai Zi Σ−1 i (yi − Xi β) ± Ci Ω, where Ci is the Cholesky decomposition of Ai such that Ci Ci = Ai Negative correlation between the N/2 pairs ⇒ smaller variance in the sample means than would be obtained from N independent simulations GL. Hickey Joint modelling of multivariate data 24 / 48
41. Introduction Model Estimation Software & Example revisited Prediction References Convergence In standard EM, convergence usually declared at (m + 1)-th iteration if one of the following criteria satisfied Relative change: ∆(m+1) rel = max |ˆ θ(m+1)−ˆ θ(m)| |ˆ θ(m)|+ 1 < 0 Absolute change: ∆(m+1) abs = max |ˆ θ(m+1) − ˆ θ(m)| < 2 for some choice of 0, 1, and 2 GL. Hickey Joint modelling of multivariate data 25 / 48
42. Introduction Model Estimation Software & Example revisited Prediction References Convergence In MCEM framework, there are 2 complications to account for GL. Hickey Joint modelling of multivariate data 26 / 48
43. Introduction Model Estimation Software & Example revisited Prediction References Convergence In MCEM framework, there are 2 complications to account for 1 spurious convergence declared due to random chance GL. Hickey Joint modelling of multivariate data 26 / 48
44. Introduction Model Estimation Software & Example revisited Prediction References Convergence In MCEM framework, there are 2 complications to account for 1 spurious convergence declared due to random chance ⇒ Solution: require convergence for 3 iterations in succession GL. Hickey Joint modelling of multivariate data 26 / 48
45. Introduction Model Estimation Software & Example revisited Prediction References Convergence In MCEM framework, there are 2 complications to account for 1 spurious convergence declared due to random chance ⇒ Solution: require convergence for 3 iterations in succession 2 estimators swamped by Monte Carlo error, thus precluding convergence GL. Hickey Joint modelling of multivariate data 26 / 48
46. Introduction Model Estimation Software & Example revisited Prediction References Convergence In MCEM framework, there are 2 complications to account for 1 spurious convergence declared due to random chance ⇒ Solution: require convergence for 3 iterations in succession 2 estimators swamped by Monte Carlo error, thus precluding convergence ⇒ Solution: increase Monte Carlo size N as algorithm moves closer towards maximizer GL. Hickey Joint modelling of multivariate data 26 / 48
47. Introduction Model Estimation Software & Example revisited Prediction References Dynamic MC size Using large N when far from maximizer = computationally inefficient Using small N when close to maximizer = unlikely to detect convergence Solution (proposed by Ripatti et al. 2002): after a ‘burn-in’ phase, calculate the coefficient of variation statistic cv(∆(m+1) rel ) = sd(∆(m−1) rel , ∆(m) rel , ∆(m+1) rel ) mean(∆(m−1) rel , ∆(m) rel , ∆(m+1) rel ) , and increase N to N + N/δ if cv(∆(m+1) rel ) > cv(∆(m) rel ) for some small positive integer δ GL. Hickey Joint modelling of multivariate data 27 / 48
48. Introduction Model Estimation Software & Example revisited Prediction References Quasi-Monte Carlo Replaces the (pseudo-)random sequence by a deterministic one Quasi-random sequences yield smaller errors than standard Monte Carlo integration methods Convergence is O (logN)r N Research on-going. . . GL. Hickey Joint modelling of multivariate data 28 / 48
49. Introduction Model Estimation Software & Example revisited Prediction References Quasi-Monte Carlo q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q OMC AMC QMC Uniform 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 U1 U2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q OMC AMC QMC Normal −4 −2 0 2 −4 −2 0 2 −4 −2 0 2 −2 0 2 4 Z1 Z2 Key: OMC = ordinary Monte Carlo; AMC = antithetic Monte Carlo; QMC = quasi-Monte Carlo GL. Hickey Joint modelling of multivariate data 29 / 48
50. Introduction Model Estimation Software & Example revisited Prediction References Standard error estimation Method 1: Bootstrap Conceptually simple + theoretically superior (Hsieh et al. 2006). . . but computationally slow! GL. Hickey Joint modelling of multivariate data 30 / 48
51. Introduction Model Estimation Software & Example revisited Prediction References Standard error estimation Method 1: Bootstrap Conceptually simple + theoretically superior (Hsieh et al. 2006). . . but computationally slow! Method 2: Empirical information matrix approximation Following McLachlan and Krishnan (2008), SE(θ) ≈ I−1/2 e (ˆ θ), where Ie(θ) = n i=1 si (θ)si (θ) − 1 n S(θ)S (θ), S(θ) = n i=1 si (θ) is the score vector for θ−λ0(t) (baseline hazards a profiled out of the likelihood) GL. Hickey Joint modelling of multivariate data 30 / 48
52. Introduction Model Estimation Software & Example revisited Prediction References joineRML + getVarCov() vcov() fixef() ranef()* AIC() BIC() confint() formula() sampleData() dynSurv()* dynLong()* print() summary() plot() sigma() coef() update() baseHaz() residuals() fitted() logLik() bootSE() Rich collection of associated methods * associated with additional plot methods mjoint() Version 0.3.0 available on CRAN https://cran.r-project.org/web/packages/joineRML/ Developmental version available on GitHub https://github.com/graemeleehickey/joineRML Parallel Computing + GL. Hickey Joint modelling of multivariate data 31 / 48
53. Introduction Model Estimation Software & Example revisited Prediction References Alternative options Pre-2017: none! 2017-onwards: joineRML: discussed today stjm: a new extension to the Stata package2 written by Michael Crowther megenreg: similar to stjm, but can handle other models rstanarm: development branch that absorbs package written by Sam Brilleman3 JMbayes: a new extension4 to the R package written by Dimitris Rizopoulos 2Crowther MJ. Joint Statistical Meeting. Seattle; 2015. 3 github.com/sambrilleman/rstanjm 4 github.com/drizopoulos/JMbayes GL. Hickey Joint modelling of multivariate data 32 / 48
54. Introduction Model Estimation Software & Example revisited Prediction References Proposed model for PBC data Longitudinal sub-model log(serBilir) = (β0,1 + b0i,1 ) + (β1,1 + b1i,1 )year + εij1, albumin = (β0,2 + b0i,2 ) + (β1,2 + b1i,2 )year + εij2, (0.1 × prothrombin)−4 = (β0,3 + b0i,3 ) + (β1,3 + b1i,3 )year + εij3, bi ∼ N6 (0, D), and εijk ∼ N(0, σ2 k ) for k = 1, 2, 3; Time-to-event sub-model λi (t) = λ0 (t) exp {γv age + W2i (t)} , W2i (t) = γbil (b0i,1 + b1i,1 t) + γalb (b0i,2 + b1i,2 t) + γpro (b0i,3 + b1i,3 t). GL. Hickey Joint modelling of multivariate data 33 / 48
55. Introduction Model Estimation Software & Example revisited Prediction References Example code data(pbc2) placebo <- subset(pbc2, drug == "placebo") fit.pbc <- mjoint( formLongFixed = list( "bil" = log(serBilir) ˜ year, "alb" = albumin ˜ year, "pro" = (0.1 * prothrombin)ˆ-4 ˜ year), formLongRandom = list( "bil" = ˜ year | id, "alb" = ˜ year | id, "pro" = ˜ year | id), formSurv = Surv(years, status2) ˜ age, data = placebo, timeVar = "year", control = list(tol0 = 0.001, burin = 400)) GL. Hickey Joint modelling of multivariate data 34 / 48
56. Introduction Model Estimation Software & Example revisited Prediction References Results Parameter Estimate SE 95% CI β0,1 0.5541 0.0858 (0.3859, 0.7223) β1,1 0.2009 0.0201 (0.1616, 0.2402) β0,2 3.5549 0.0356 (3.4850, 3.6248) β1,2 -0.1245 0.0101 (-0.1444, -0.1047) β0,3 0.8304 0.0212 (0.7888, 0.8719) β1,3 -0.0577 0.0062 (-0.0699, -0.0456) γv 0.0462 0.0151 (0.0166, 0.0759) γbil 0.8181 0.2046 (0.4171, 1.2191) γalb -1.7060 0.6181 (-2.9173, -0.4946) γpro -2.2085 1.6070 (-5.3582, 0.9412) GL. Hickey Joint modelling of multivariate data 35 / 48
57. Introduction Model Estimation Software & Example revisited Prediction References Results Parameter Estimate SE 95% CI β0,1 0.5541 0.0858 (0.3859, 0.7223) β1,1 0.2009 0.0201 (0.1616, 0.2402) β0,2 3.5549 0.0356 (3.4850, 3.6248) β1,2 -0.1245 0.0101 (-0.1444, -0.1047) β0,3 0.8304 0.0212 (0.7888, 0.8719) β1,3 -0.0577 0.0062 (-0.0699, -0.0456) γv 0.0462 0.0151 (0.0166, 0.0759) γbil 0.8181 0.2046 (0.4171, 1.2191) γalb -1.7060 0.6181 (-2.9173, -0.4946) γpro -2.2085 1.6070 (-5.3582, 0.9412) GL. Hickey Joint modelling of multivariate data 35 / 48
58. Introduction Model Estimation Software & Example revisited Prediction References Results Effect of multivariate inference over univariate joint model: Parameter Model Estimate 95% CI γbil UV 1.2182 (0.9789, 1.6130) γbil MV 0.8181 (0.4171, 1.2191) γalb UV -3.0770 (-4.4865, -2.3466) γalb MV -1.7060 (-2.9173, -0.4946) γpro UV -7.2078 (-10.5410, -5.3917) γpro MV -2.2085 (-5.3582, 0.9412) UV = univariate joint model (fitted with joineR package); MV = multivariate joint model GL. Hickey Joint modelling of multivariate data 36 / 48
59. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction So far we have only discussed inference from joint models How we can use them for prediction? Predict what? 1 Failure probability at time u > t given longitudinal data observed up until time t 2 Longitudinal trajectories at time u > t given longitudinal data observed up until time t GL. Hickey Joint modelling of multivariate data 37 / 48
60. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction: example Bivariate joint model We will consider the PBC data again (as above) with K = 2 biomarkers only: serurm bilirubin (log-transformed) and albumin (untransformed), since prothrombin time was non-significant in the trivariate model GL. Hickey Joint modelling of multivariate data 38 / 48
61. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction: survival For a new subject i = n + 1, we want to calculate P[T∗ n+1 ≥ u | T∗ n+1 > t, yn+1; θ] = E Sn+1 (u | W2,n+1(u, bn+1; θ); θ) Sn+1 (t | W2,n+1(t, bn+1; θ); θ) , where W2i (t, bi ; θ) = {W2i (s, vi ; θ); 0 ≤ s < t} and the expectation is taken with respect to the distribution p(bn+1 | T∗ n+1 > t, yn+1; θ) GL. Hickey Joint modelling of multivariate data 39 / 48
62. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction: survival Rizopoulos (2011) proposed two estimators for this: 1 A first-order approximation P[T∗ n+1 ≥ u | T∗ n+1 > t, yn+1 ; θ] ≈ Sn+1 u | W2,n+1 (u, ˆ bn+1 ; ˆ θmle ); ˆ θmle Sn+1 t | W2,n+1 (t, ˆ bn+1 ; ˆ θmle ); ˆ θmle , where ˆ bn+1 is the mode of p(bn+1 | T∗ n+1 > t, yn+1 ; θ) 2 A simulated scheme 1 Draw θ(l) ∼ N(ˆ θmle , V (ˆ θmle )) 2 Draw b(l) n+1 ∼ p(bn+1 | T∗ n+1 > t, yn+1 ; θ) [Metropolis-Hastings] 3 Calculate Sn+1 u | W2,n+1(u,b(l) n+1 ;θ(l));θ(l) Sn+1 t | W2,n+1(t,b(l) n+1 ;θ(l));θ(l) 4 Repeat Steps 1–3 l = 2, . . . , L times GL. Hickey Joint modelling of multivariate data 40 / 48
63. Introduction Model Estimation Software & Example revisited Prediction References Example code # New patient nd <- subset(placebo, id == "11") # patient 11 # First-order prediction (default) pred1 <- dynSurv(fit.pbc, nd[1:5, ]) pred1 plot(pred1) # Simulated prediction pred2 <- dynSurv(fit.pbc, nd[1:5, ], type = "simulated", scale = 2) pred2 plot(pred2) GL. Hickey Joint modelling of multivariate data 41 / 48
64. Introduction Model Estimation Software & Example revisited Prediction References Dynamic predicton: survival GL. Hickey Joint modelling of multivariate data 42 / 48 0.20 0.25 0.30 0.35 0.40 0.45 data.t$tk[[k]] log(bilirubin) 2.5 3.0 3.5 4.0 4.5 5.0 5.5 data.t$tk[[k]] albumin 0 0 2 4 6 8 10 12 14 xpts ypts 0.0 0.2 0.4 0.6 0.8 1.0 Time Event−free probability
65. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction: longitudinal For a new subject i = n + 1, we want to calculate E yn+1(u) | T∗ n+1 > t, yn+1; θ = Xn+1 (u)β + Zn+1 (u)E[bn+1], GL. Hickey Joint modelling of multivariate data 43 / 48
66. Introduction Model Estimation Software & Example revisited Prediction References Dynamic prediction: longitudinal Again, we can use the same estimation proposals: 1 A first-order approximation E [yn+1 (u) | T∗ n+1 > t, yn+1 ; θ] ≈ Xn+1 (u)ˆ β + Zn+1 (u)ˆ bn+1, where ˆ bn+1 is the mode of p(bn+1 | T∗ n+1 > t, yn+1 ; θ) 2 A simulated scheme 1 Draw θ(l) ∼ N(ˆ θmle , V (ˆ θmle )) 2 Draw b(l) n+1 ∼ p(bn+1 | T∗ n+1 > t, yn+1 ; θ) [Metropolis-Hastings] 3 Calculate Xn+1 (u)β(l) + Zn+1 (u)b(l) n+1 4 Repeat Steps 1–3 l = 2, . . . , L times GL. Hickey Joint modelling of multivariate data 44 / 48
67. Introduction Model Estimation Software & Example revisited Prediction References Example code # First-order prediction (default) pred1 <- dynLong(fit.pbc, nd[1:5, ]) pred1 plot(pred1) # Simulated prediction pred2 <- dynLong(fit.pbc, nd[1:5, ], type = "simulated", scale = 2) pred2 plot(pred2) GL. Hickey Joint modelling of multivariate data 45 / 48
68. Introduction Model Estimation Software & Example revisited Prediction References Dynamic predicton: longitudinal GL. Hickey Joint modelling of multivariate data 46 / 48 0 1 2 3 4 5 xpts log(bilirubin) 1.5 2.0 2.5 3.0 3.5 4.0 xpts albumin 0 2 4 6 8 10 12 14 Time
69. Introduction Model Estimation Software & Example revisited Prediction References Open challenges How can we incorporate high-dimensional K? E.g. K = 10? Data reduction techniques: can we project high-dimensional K onto a lower order plane? Speed-up calculations using approximations (e.g. Laplace approximations) GL. Hickey Joint modelling of multivariate data 47 / 48
70. Introduction Model Estimation Software & Example revisited Prediction References References Murtaugh, Paul A et al. (1994). Primary biliary cirrhosis: prediction of short-term survival based on repeated patient visits. Hepatology 20(1), pp. 126–134. Henderson, R et al. (2000). Joint modelling of longitudinal measurements and event time data. Biostatistics 1(4), pp. 465–480. Laird, NM and Ware, JH (1982). Random-effects models for longitudinal data. Biometrics 38(4), pp. 963–74. Dempster, AP et al. (1977). Maximum likelihood from incomplete data via the EM algorithm. J Roy Stat Soc B 39(1), pp. 1–38. Wei, GC and Tanner, MA (1990). A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. J Am Stat Assoc 85(411), pp. 699–704. Ripatti, S et al. (2002). Maximum likelihood inference for multivariate frailty models using an automated Monte Carlo EM algorithm. Life Dat Anal 8(2002), pp. 349–360. Hsieh, F et al. (2006). Joint modeling of survival and longitudinal data: Likelihood approach revisited. Biometrics 62(4), pp. 1037–1043. McLachlan, GJ and Krishnan, T (2008). The EM Algorithm and Extensions. Second. Wiley-Interscience. Rizopoulos, Dimitris (2011). Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. Biometrics 67(3), pp. 819–829. GL. Hickey Joint modelling of multivariate data 48 / 48