Upgrade to Pro — share decks privately, control downloads, hide ads and more …

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

Graeme Hickey
November 13, 2017

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

Presented at Durham University, Department of Mathematical Sciences

Graeme Hickey

November 13, 2017
Tweet

More Decks by Graeme Hickey

Other Decks in Research

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 [email protected] 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