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

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

Presented at Erasmus MC, Department of Biostatistics, Rotterdam

3691d1dba94a59d161a84382029b09c0?s=128

Graeme Hickey

March 14, 2018
Tweet

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 14th March 2018 GL. Hickey Joint modelling of multivariate data 1 / 55
  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 / 55
  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 / 55
  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 / 55
  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 / 55
  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 / 55
  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 / 55
  8. Introduction Model Estimation Software & Example revisited Prediction References Clinical

    example Primary biliary cirrhosis (PBC) is a chronic liver disease characterized by inflammatory destruction of the small bile ducts, which eventually leads to cirrhosis of the liver and death 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 5 / 55
  9. 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 The natural disease progression1 of untreated patients with PBC enrolled into the trial 1Biomarker transformations chosen according to univariate Box-Cox power transformations of linear model Y = β0 + β1 t + βsubject + GL. Hickey Joint modelling of multivariate data 6 / 55
  10. 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 7 / 55
  11. Introduction Model Estimation Software & Example revisited Prediction References Outline

    of the rest of the talk Model development Extend existing models for multivariate data and develop practical fitting strategies GL. Hickey Joint modelling of multivariate data 8 / 55
  12. Introduction Model Estimation Software & Example revisited Prediction References Outline

    of the rest of the talk Model development Extend existing models for multivariate data and develop practical fitting strategies Dynamic prediction Extend modelling framework for dynamic prediction purposes and explore utility of multiple longitudinal outcomes GL. Hickey Joint modelling of multivariate data 8 / 55
  13. Introduction Model Estimation Software & Example revisited Prediction References Outline

    of the rest of the talk Model development Extend existing models for multivariate data and develop practical fitting strategies Dynamic prediction Extend modelling framework for dynamic prediction purposes and explore utility of multiple longitudinal outcomes Software development Present an –package that implements these methods GL. Hickey Joint modelling of multivariate data 8 / 55
  14. 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 9 / 55
  15. 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 10 / 55
  16. 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 10 / 55
  17. 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 11 / 55
  18. 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 12 / 55
  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 Sources of correlation: GL. Hickey Joint modelling of multivariate data 12 / 55
  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 Sources of correlation: 1 Within-subject correlation between longitudinal measurements: bik ∼ N(0, Dkk) GL. Hickey Joint modelling of multivariate data 12 / 55
  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 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 12 / 55
  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 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-models2: W2i (t) = K k=1 γykW (k) 1i (t) 2Extends model proposed Henderson et al. (2000) GL. Hickey Joint modelling of multivariate data 12 / 55
  23. 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) = γy ∂ ∂t {µi (t) + W1i (t)} Linear combination of different functions3 Use of time-varying association parameters 4 . . . 3e.g. Andrinopoulou and Rizopoulos (2016) 4e.g. Andrinopoulou, Eilers, et al. (2017) GL. Hickey Joint modelling of multivariate data 13 / 55
  24. 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 14 / 55
  25. 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 15 / 55
  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 ), 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 15 / 55
  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 (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 15 / 55
  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 (bi | θ) = (2π)− r 2 |D|−1 2 exp − 1 2 bi D−1bi , with r = dim(bi ) GL. Hickey Joint modelling of multivariate data 15 / 55
  29. 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 16 / 55
  30. 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 17 / 55
  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 M-step. We maximise Q(θ | ˆ θ(m)) with respect to θ. namely, ˆ θ(m+1) = arg max θ Q(θ | ˆ θ(m)) GL. Hickey Joint modelling of multivariate data 17 / 55
  32. 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 18 / 55
  33. 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 19 / 55
  34. Introduction Model Estimation Software & Example revisited Prediction References M-step:

    non-closed form estimators I(γ) is computationally expensive to calculate I(γ) = n i=1 J j=1 ˆ λ0(tj )I(Ti ≥ tj )E ˜ v⊗2 i (tj ) exp{˜ vi (tj )γ} − J j=1 ˆ λ0(tj )2Γ(tj ) n i=1 δi I(Ti = tj ) , where Γ(tj ) = n i=1 E ˜ vi (tj ) exp{˜ vi (tj )γ} I(Ti ≥ tj ) ⊗2 Alternative option: use a Gauss-Newton-like update by approximating I(γ) using only score (i.e. first derivative) components (more later. . . ) GL. Hickey Joint modelling of multivariate data 20 / 55
  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 (Recall: EM algorithm has first-order convergence) Instead, we use the Monte Carlo Expectation-Maximization (MCEM)5 M-step updates remain the same 5Wei and Tanner (1990) GL. Hickey Joint modelling of multivariate data 21 / 55
  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 / 55
  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 / 55
  38. Introduction Model Estimation Software & Example revisited Prediction References Flavours

    of MC integration We consider 3 different flavours of Monte Carlo (MC) integration: 1 Ordinary MC 2 Variance reduction MC Antithetic variates 3 Quasi-MC Sobol sequence Halton sequence GL. Hickey Joint modelling of multivariate data 23 / 55
  39. Introduction Model Estimation Software & Example revisited Prediction References Variance

    reduction Direct sampling from the MVN distribution of bi | yi ; θ offers convergence at a rate of O(N−1/2), which is indep. of K and r = dim(bi ) Alternative: Antithetic simulation (default) 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 indep. simulations GL. Hickey Joint modelling of multivariate data 24 / 55
  40. Introduction Model Estimation Software & Example revisited Prediction References Quasi-Monte

    Carlo Replaces the (pseudo-)random sequence by a deterministic one Sobol sequence Halton sequence . . . Quasi-random sequences yield smaller errors than standard MC integration methods Convergence is O (logN)r N Early results promising + on-going research via intensive simulation . . . GL. Hickey Joint modelling of multivariate data 25 / 55
  41. 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 26 / 55
  42. 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 27 / 55
  43. Introduction Model Estimation Software & Example revisited Prediction References Convergence

    In MCEM framework, there are 3 complications to account for GL. Hickey Joint modelling of multivariate data 28 / 55
  44. Introduction Model Estimation Software & Example revisited Prediction References Convergence

    In MCEM framework, there are 3 complications to account for 1. Convergence might occur at a local maximum6 “Solution”: use sensible initial values θ0 6Applies to the conventional EM algorithm also. GL. Hickey Joint modelling of multivariate data 28 / 55
  45. Introduction Model Estimation Software & Example revisited Prediction References Convergence

    In MCEM framework, there are 3 complications to account for 1. Convergence might occur at a local maximum6 “Solution”: use sensible initial values θ0 2. Spurious convergence declared due to random chance “Solution”: require convergence for 3 iterations in succession 6Applies to the conventional EM algorithm also. GL. Hickey Joint modelling of multivariate data 28 / 55
  46. Introduction Model Estimation Software & Example revisited Prediction References Convergence

    In MCEM framework, there are 3 complications to account for 1. Convergence might occur at a local maximum6 “Solution”: use sensible initial values θ0 2. Spurious convergence declared due to random chance “Solution”: require convergence for 3 iterations in succession 3. Estimators swamped by MC error, thus precluding convergence “Solution”: increase MC size N as algorithm moves closer towards maximizer 6Applies to the conventional EM algorithm also. GL. Hickey Joint modelling of multivariate data 28 / 55
  47. Introduction Model Estimation Software & Example revisited Prediction References Initial

    values Multivariate longitudinal sub-model Ordinary EM-algorithm: E-step is closed-form, so can be made quick Time-to-event sub-model If balanced (i.e. tijk = tij ∀k), then fit a two-stage model using the BLUPs for the MV-LMM (above) using (Anderson-Gill counting process) Cox model If unbalanced, γy = 0 and γv estimated using a standard Cox PH model Alternatively, specify them manually from expert knowledge GL. Hickey Joint modelling of multivariate data 29 / 55
  48. 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 30 / 55
  49. 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 31 / 55
  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! 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 are profiled out of the likelihood) GL. Hickey Joint modelling of multivariate data 31 / 55
  51. 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 32 / 55
  52. Introduction Model Estimation Software & Example revisited Prediction References Alternative

    package options JMbayes: extends the existing R package (by Dimitris Rizopoulos, Erasmus MC) rstanarm: absorbs the now defunct R package rstanjm (by Sam Brilleman, Monash University) stjm: extends the existing Stata package (by Michael Crowther, University of Leicester) megenreg: similar to stjm, but can handle other models (also by Michael Crowther, University of Leicester) GL. Hickey Joint modelling of multivariate data 33 / 55
  53. 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 34 / 55
  54. 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, burnin = 400)) GL. Hickey Joint modelling of multivariate data 35 / 55
  55. 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 36 / 55
  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 36 / 55
  57. 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 37 / 55
  58. Introduction Model Estimation Software & Example revisited Prediction References Dynamic

    prediction So far we have only discussed inference from joint models Dynamic prediction Can we predict: 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 Utility of multivariate data Do we gain from using multiple longitudinal outcomes? GL. Hickey Joint modelling of multivariate data 38 / 55
  59. Introduction Model Estimation Software & Example revisited Prediction References Dynamic

    prediction: survival For a new subject i = n + 1 with available longitudinal measurements yn+1{t} up to time t, we want to calculate P[T∗ n+1 ≥ u | T∗ n+1 > t, yn+1{t}; θ] = E Sn+1 (u | bn+1; θ) Sn+1 (t | bn+1; θ) , where the expectation is taken with respect to the distribution p(bn+1 | T∗ n+1 > t, yn+1{t}; θ), and Sn+1 (u | bn+1; θ) = exp − u 0 λ0(s) exp{vn+1 (s)γv + Wn+1(s)}ds GL. Hickey Joint modelling of multivariate data 39 / 55
  60. Introduction Model Estimation Software & Example revisited Prediction References Dynamic

    prediction: survival Rizopoulos (2011) proposed two estimators for this: 1. First-order approximation P[T∗ n+1 ≥ u | T∗ n+1 > t, yn+1{t}; θ] ≈ Sn+1 u | ˆ bn+1 ; ˆ θ Sn+1 t | ˆ bn+1 ; ˆ θ , where ˆ bn+1 is the mode of p(bn+1 | T∗ n+1 > t, yn+1{t}; ˆ θ) 2. Simulated scheme 1 Draw θ(l) ∼ N(ˆ θ, V (ˆ θ)) 2 Draw b(l) n+1 ∼ p(bn+1 | T∗ n+1 > t, yn+1{t}; θ(l)) [Metropolis-Hastings] 3 Calculate Sn+1 u | b(l) n+1 ;θ(l) Sn+1 t | b(l) n+1 ;θ(l) 4 Repeat Steps 1–3 for l = 2, . . . , L GL. Hickey Joint modelling of multivariate data 40 / 55
  61. 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 / 55
  62. Introduction Model Estimation Software & Example revisited Prediction References Example

    code First order approximation 0.2 0.3 0.4 0.5 0.6 0.7 0.8 data.t$tk[[k]] log, serBilir 3.7 3.8 3.9 4.0 4.1 data.t$tk[[k]] albumin 0.45 0.50 0.55 0.60 0.65 data.t$tk[[k]] ^, (0.1 * prothrombin), −4 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 GL. Hickey Joint modelling of multivariate data 42 / 55
  63. Introduction Model Estimation Software & Example revisited Prediction References Example

    code Simulated scheme 0.2 0.3 0.4 0.5 0.6 0.7 0.8 data.t$tk[[k]] log, serBilir 3.7 3.8 3.9 4.0 4.1 data.t$tk[[k]] albumin 0.45 0.50 0.55 0.60 0.65 data.t$tk[[k]] ^, (0.1 * prothrombin), −4 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 GL. Hickey Joint modelling of multivariate data 42 / 55
  64. Introduction Model Estimation Software & Example revisited Prediction References Example

    code Simulated scheme + constrained B-spline smoother 0.2 0.3 0.4 0.5 0.6 0.7 0.8 data.t$tk[[k]] log, serBilir 3.7 3.8 3.9 4.0 4.1 data.t$tk[[k]] albumin 0.45 0.50 0.55 0.60 0.65 data.t$tk[[k]] ^, (0.1 * prothrombin), −4 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 GL. Hickey Joint modelling of multivariate data 42 / 55
  65. Introduction Model Estimation Software & Example revisited Prediction References Truly

    dynamic prediction GL. Hickey Joint modelling of multivariate data 43 / 55 0.20 0.25 0.30 0.35 0.40 0.45 data.t$tk[[k]] log, serBilir 2.5 3.0 3.5 4.0 4.5 5.0 5.5 data.t$tk[[k]] albumin 0.3 0.4 0.5 0.6 data.t$tk[[k]] ^, (0.1 * prothrombin), −4 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
  66. Introduction Model Estimation Software & Example revisited Prediction References Comparison

    with JM Separate univariate joint models for each biomarker (B-spline model for λ0(t)) 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Univariate joint model {B} Time Pr(Ti ≥ u | Ti > t, y ~ i (t)) R package JM joineRML 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Univariate joint model {A} Time Pr(Ti ≥ u | Ti > t, y ~ i (t)) R package JM joineRML 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Univariate joint model {P} Time Pr(Ti ≥ u | Ti > t, y ~ i (t)) R package JM joineRML GL. Hickey Joint modelling of multivariate data 44 / 55
  67. Introduction Model Estimation Software & Example revisited Prediction References Value

    of multiple biomarkers Build all 32 − 1 = 7 joint models: i.e. from trivariate ({B, A, P}) to univariate ({B}, {A}, {P}) For all patients still in the risk set at landmark time (t), predict their failure probability at horizon time (s), i.e. P[t < T∗ i ≤ s + t | T∗ i > t, yi {t}; θ] Calculate concordance statistic7 (Therneau and Watson 2015) between prediction and whether subject had event in (t, s + t] 7Using the survConcordance() function in the R survival package GL. Hickey Joint modelling of multivariate data 45 / 55
  68. Introduction Model Estimation Software & Example revisited Prediction References Ensemble

    of joint models (patient #11) GL. Hickey Joint modelling of multivariate data 46 / 55 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Visit 1 (t = 0) Time (years) Predicted survival {B, A, P} {B, A} {B, P} {A, P} {B} {A} {P}
  69. Introduction Model Estimation Software & Example revisited Prediction References Discrimination

    (concordance) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2 4 1 2 4 8 Landmark time (t; years) Horizon time (s; years) 0.6 0.7 0.8 0.9 1.0 0.6 0.7 0.8 0.9 1.0 0.6 0.7 0.8 0.9 1.0 0.6 0.7 0.8 0.9 1.0 Predicted survival Model q q q q q q q {B, A, P} {B, A} {B, P} {A, P} {B} {A} {P} GL. Hickey Joint modelling of multivariate data 47 / 55
  70. 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{t}; θ = Xn+1 (u)β + Zn+1 (u)E[bn+1], GL. Hickey Joint modelling of multivariate data 48 / 55
  71. 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{t}; θ] ≈ 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(ˆ θ, V (ˆ θ)) 2 Draw b(l) n+1 ∼ p(bn+1 | T∗ n+1 > t, yn+1{t}; θ) [Metropolis-Hastings] 3 Calculate Xn+1 (u)β(l) + Zn+1 (u)b(l) n+1 4 Repeat Steps 1–3 for l = 2, . . . , L GL. Hickey Joint modelling of multivariate data 49 / 55
  72. Introduction Model Estimation Software & Example revisited Prediction References Example

    code # NB: prothrombin time is now dropped from model and refitted # 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 50 / 55
  73. Introduction Model Estimation Software & Example revisited Prediction References Dynamic

    prediction: longitudinal GL. Hickey Joint modelling of multivariate data 51 / 55 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
  74. 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) or alternative E-steps (e.g. importance sampling) GL. Hickey Joint modelling of multivariate data 52 / 55
  75. Introduction Model Estimation Software & Example revisited Prediction References Bedankt

    voor het luisteren! Questions? Slides available from www.glhickey.com GL. Hickey Joint modelling of multivariate data 53 / 55
  76. Introduction Model Estimation Software & Example revisited Prediction References References

    I 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. Andrinopoulou, Eleni-Rosalina and Rizopoulos, Dimitris (2016). Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures. Statistics in Medicine 35(26), pp. 4813–4823. Andrinopoulou, Eleni-Rosalina, Eilers, Paul H C, et al. (2017). Improved dynamic predictions from joint models of longitudinal and survival data with time-varying effects using P-splines. Biometrics. 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. GL. Hickey Joint modelling of multivariate data 54 / 55
  77. Introduction Model Estimation Software & Example revisited Prediction References References

    II 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. Therneau, Terry M and Watson, David A (2015). The concordance statistic and the Cox model. Tech. rep. Rochester, Minnesota: Department of Health Science Research, Mayo Clinic, Technical Report Series No. 85. GL. Hickey Joint modelling of multivariate data 55 / 55