# Joint Models for Longitudinal Measurements and Competing Risks Failure Time Data

#### 2021-03-22

Joint analysis of longitudinal outcomes and survival times has gained a lot of attention in recent years. There have been extended to handle competing risks for both continuous and ordinal outcomes (Elashoff, Li, and Li 2008, @MR2758452). This vignette offers a brief introduction to the R package JMcmprsk, which implements the methods proposed to deal with such joint models, and two competing risks are assumed. The data sets for generating the initial values are provided.

## Joint modeling framework

#### Continuous outcomes

Let $$Y_i(t)$$ be the longitudinal outcome measured at time $$t$$ for subject $$i$$, $$i=1,2,\cdots,n$$, and $$n$$ is the total number of subjects in study. Let $$C_i=(T_i,D_i)$$ denote the competing risks data on subject $$i$$, where $$T_i$$ is the failure time or censoring time, and $$D_i$$ takes value in $$\{0,1,\cdots,g\}$$, with $$D_i=0$$ indicating a censored event and $$D_i=k$$ showing that subject $$i$$ fails from the $$k$$th type of failure, where $$k=1,\cdots,g$$.

The joint model is specified in terms of the following two linked components: $\begin{eqnarray*} Y_i(t)&=&X_i^{(1)}(t)^\top \beta+\tilde X_i^{(1)}(t)^\top b_i+\epsilon_i(t),\\ \lambda_k(t)&=&\lambda_{0k}(t)\exp(X_i^{(2)}(t)^\top \gamma_k+\nu_k u_i),~~\mbox{for}~~k=1,\cdots,g, \end{eqnarray*}$ where $$X_i^{(1)}(t)$$, $$X_i^{(2)}(t)$$ denote the covaraites for the fixed-effects $$\beta$$ and $$\gamma_k$$, $$\tilde X_i^{(1)}(t)$$ denotes the covaraites for the random-effects $$b_i$$, and $$\epsilon_i(t)\sim N(0,\sigma^2)$$ for all $$t\geq 0$$. The parameter $$\nu_1$$ is set to 1 to ensure identifiability. We assume that $$b_i$$ is independent of $$\epsilon_i(t)$$ and that $$\epsilon_i(t_1)$$ is independent of $$\epsilon_i(t_2)$$ for any $$t_1\neq t_2$$. We further assume the random effects $$b_i$$ and $$u_i$$ jointly have a multivariate normal distribution, denoted by $$\theta_i\sim N(0,\Sigma)$$, where $$\Sigma=(\Sigma_{b},\Sigma_{bu}^\top;\Sigma_{bu},\sigma_u)$$.

Denote $$\Psi$$ as the unknown parameters from the joint models. We propose to obtain the maximum likelihood estimate of $$\Psi$$ through an EM algorithm. The complete data likelihood is $\begin{eqnarray*} &&L(\Psi;Y,C,\theta)\\ &&\propto \Pi_{i=1}^n\Big[\Pi_{j=1}^{n_i}\frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{1}{2\sigma^2}(Y_{ij}-X_i^{(1)}(t_{ij})^\top\beta-\tilde X_i^{(1)}(t_{ij})^\top b_i)^2)\Big]\\ &&\times \Pi_{k=1}^g\lambda_k(T_i)^{I(D_i=k)}\exp\Big\{-\int_0^{T_i}\sum_{k=1}^g\lambda_k(t)dt\Big\}\\ &&\times \frac{1}{\sqrt{(2\pi)^d|\Sigma|}}\exp(-\frac{1}{2}\theta_i^\top\Sigma^{-1}\theta_i). \end{eqnarray*}$ In the E-step, we need to calculate the expected value of all the functions of $$\theta$$, The integral does not have a closed-form solution, and thus a numerical method must be employed for its evaluation. Standard option is the Gaussian quadratic rules. In the M-step, we can update $$\Psi$$ through maximizing the function obtained from the E-step.

#### Ordinal outcomes

Let $$Y_{ij}$$ denote the $$j$$th response measured on subject $$i$$, where $$i=1,\cdots,n$$, $$j=1,\cdots,n_i$$, and $$Y_{ij}$$ takes values in $$\{1,\cdots,K\}$$. The competing risks failure times on subject $$i$$ is $$(T_i,D_i)$$, and the meaning is the same as in subsection of “Continuous outcomes”.

we propose the following partial proportional odds model for $$Y_{ij}$$ $\begin{eqnarray*} P(Y_{ij}\leq k|X_{ij},\tilde X_{ij},W_{ij},b_i)=\frac{1}{1+\exp(-\theta_{k}-X_{ij}\beta-\tilde X_{ij}\alpha_{k}-W_{ij}^\top b_i)}. \end{eqnarray*}$ where $$k=1,\cdots,K-1$$, $$X_{ij}$$ and $$\tilde X_{ij}$$ are $$p\times 1$$ and $$s\times 1$$ vectors of covaraites for the fixed-effect $$\beta$$ and $$\alpha_{k}$$, with $$\alpha_{1}=0$$, and $$\tilde X_{ij}$$ is a subset of $$X_{ij}$$ for which the proportional odds assumption may not be satisfied. The $$q\times 1$$ vector $$b_i$$ represents random effects of the associated covaraites $$W_{ij}$$.

The distribution of the competing risks failure times $$(T_i,D_i)$$ are assumed to take the form of the following cause-specific hazards frailty model: $\begin{eqnarray*} \lambda_d(t|Z_i(t),u_i)&=&\lambda_{0d}(t)\exp(Z_i(t)^\top \gamma_d+\nu_d u_i),~~\mbox{for}~~d=1,\cdots,g, \end{eqnarray*}$ where the $$l\times 1$$ vector $$\gamma_d$$ and $$\nu_d$$ are cause-specific coefficients for the covariates $$Z_i(t)$$ and the random effects $$u_i$$, respectively.

The parameter $$\nu_1$$ is set to 1 to ensure identifiability. We assume the random effects $$b_i$$ and $$u_i$$ jointly have a multivariate normal distribution, denoted by $$a_i\sim N(0,\Sigma)$$.

Denote $$\Psi$$ as the unknown parameters from the joint models. We propose to obtain the maximum likelihood estimate of $$\Psi$$ through an EM algorithm. The complete data likelihood is $\begin{eqnarray*} &&L(\Psi;Y,C,a)\\ &&\propto \Pi_{i=1}^n\Big[\Pi_{j=1}^{n_i}\Pi_{k=1}^{K}\{\pi_{ij}(k)-\pi_{ij}(k-1)\}^{I(Y_{ij}=k)}\Big]\\ &&\times \Pi_{d=1}^g\lambda_d(T_i)^{I(D_i=d)}\exp\Big\{-\int_0^{T_i}\sum_{k=1}^d\lambda_d(t)dt\Big\}\\ &&\times \frac{1}{\sqrt{(2\pi)^{q+1}|\Sigma|}}\exp(-\frac{1}{2}a_i^\top\Sigma^{-1}a_i). \end{eqnarray*}$ where $$\pi_{ij}(k)$$ stands for the probability that $$Y_{ij}\leq k$$ given the covariates and the random effects. The implementation of EM algorithm is the same as that of subection of “Continuous outcomes”.

## The R package JMcmprsk

The function that fits continuous outcomes in JMcmprsk is called jmc(). As an illustrative example of jmc(), we consider Scleroderma Lung Study (Tashkin et al. 2006), a double-blinded, randomized clinical trial to evaluate effectiveness of oral cyclophosphamide (CYC) versus placebo in the treatment of lung disease due to scleroderma. This study includes 158 patients, the primary outcome is forced vital capacity (FVC, as % predicted) determined at 3-month intervals from the baseline. The event of interest is the time-to-treatment failure or death. We consider two factors baseline %FVC ($$FVC_0$$) and baseline lung fibrosis ($$FIB_0$$) and two risks, informative and noninformative. The model setups are as follows, $\begin{eqnarray*} \%FVC_{ij}&=&\beta_1+\beta_2t_{ij}+\beta_3FVC_{0i}+\beta_4FIB_{0i}+\beta_5CYC_i\\ &&+\beta_6FVC_{0i}\times CYC_i+\beta_7FIB_{0i}\times CYC_i+\beta_8 t_{ij}\times CYC_i+b_it_{ij}+\epsilon, \end{eqnarray*}$ and the cause-specific hazards frailty models are $\begin{eqnarray*} \lambda_1(t)=\lambda_{01}(t)\exp(\gamma_{11}FVC_{0i}+\gamma_{12}FIB_{0i}+\gamma_{13}CYC_i+\gamma_{14}FVC_{0i}\times CYC_i+\gamma_{15}FIB_{0i}\times CYC_i+u_i)\\ \lambda_2(t)=\lambda_{02}(t)\exp(\gamma_{21}FVC_{0i}+\gamma_{22}FIB_{0i}+\gamma_{23}CYC_i+\gamma_{24}FVC_{0i}\times CYC_i+\gamma_{25}FIB_{0i}\times CYC_i+\nu_2u_i), \end{eqnarray*}$

## Examples

We first load the package and the data.

require(JMcmprsk)
set.seed(123)
data(lung)
cread <- unique(lung[, c(1, 12, 13, 6:10)])

The model can be specified by the function jmc():

FE = c("time", "FVC0", "FIB0", "CYC", "FVC0.CYC",
"FIB0.CYC", "time.CYC"),
RE = "linear", ID = "ID",cate = NULL, intcpt = 0,
quad.points = 20, quiet = FALSE)
jmcfit

We can also call the internal function jmc_0() to get the same model fit. Below is the illustrative example of using the same dataset.

require(JMcmprsk)
data(lung)
lungY <- lung[, c(2:11)]
lungC <- unique(lung[, c(1, 12, 13, 6:10)])
lungC <- lungC[, -1]
lungM <- data.frame(table(lung$ID)) lungM <- as.data.frame(lungM[, 2]) The number of rows in lungY is the total number of measurements for all subjects in the study. The columns in yfile should start with the longitudinal outcome (column 1), the covariates for the random effects, and then the covariates for the fixed effects. For lungC, the survival / censoring time is included in the first column, and the failure type coded as 0 (censored events), 1 (risk 1), or 2 (risk 2) is given in the second column. Two competing risks are assumed. The covariates are included in the third column and on. lungM is to indicate the number of longitudinal measurements per subject. The number of rows equals to the number of subjects. Hence, the model can be specified by the function jmc_0(): res1=jmc_0(p=8, lungY, lungC, lungM, point=20, do.trace = FALSE, type_file = FALSE) res1 with $$p$$ the dimension of fixed effects (include intercept) in yfile, the option “point” is the number of points used to approximate the integral in the E-step, default is 20, and “do.trace” is used to control whether the programm prints the iteration details. “type_file” is to control whether the datasets are called from file paths / data frames, default is TRUE, i.e., file paths. If we would like to see a concise summary of the result we can simply type #because of the long running time, we save the jmo and jmc results within the package fitfile=system.file("extdata", "runfit.RData", package = "JMcmprsk") load(fitfile) jmcfit ## ## Call: ## jmc(long_data = yread, surv_data = cread, out = "FVC", FE = c("time", "FVC0", "FIB0", "CYC", "FVC0.CYC", "FIB0.CYC", "time.CYC"), RE = "linear", ID = "ID", cate = NULL, intcpt = 0, quad.points = 20, quiet = FALSE) ## ## Data Summary: ## Number of observations: 715 ## Number of groups: 140 ## ## Proportion of competing risks: ## Risk 1 : 10 % ## Risk 2 : 22.86 % ## ## Numerical intergration: ## Method: standard Guass-Hermite quadrature ## Number of quadrature points: 20 ## ## Model Type: joint modeling of longitudinal continuous and competing risks data ## ## Model summary: ## Longitudinal process: linear mixed effects model ## Event process: cause-specific Cox proportional hazard model with unspecified baseline hazard ## ## Loglikelihood: -3799.044 ## ## Longitudinal sub-model fixed effects: FVC ~ time + FVC0 + FIB0 + CYC + FVC0.CYC + FIB0.CYC + time.CYC ## ## Estimate Std. Error 95% CI Pr(>|Z|) ## Longitudinal: ## Fixed effects: ## intercept 66.0415 0.7541 ( 64.5634, 67.5196) 0.0000 ## time -0.0616 0.0790 (-0.2165, 0.0932) 0.4353 ## FVC0 0.9017 0.0365 ( 0.8302, 0.9732) 0.0000 ## FIB0 -1.7780 0.5605 (-2.8767,-0.6793) 0.0015 ## CYC 0.0150 0.9678 (-1.8819, 1.9119) 0.9876 ## FVC0.CYC 0.1380 0.0650 ( 0.0106, 0.2654) 0.0338 ## FIB0.CYC 1.7088 0.7643 ( 0.2109, 3.2068) 0.0254 ## time.CYC 0.1278 0.1102 (-0.0883, 0.3438) 0.2464 ## Random effects: ## sigma^2 22.7366 0.6575 ( 21.4478, 24.0253) 0.0000 ## ## Survival sub-model fixed effects: Surv(surv, failure_type) ~ FVC0 + FIB0 + CYC + FVC0.CYC + FIB0.CYC ## ## Estimate Std. Error 95% CI Pr(>|Z|) ## Survival: ## Fixed effects: ## FVC0_1 0.0187 0.0326 (-0.0452, 0.0826) 0.5660 ## FIB0_1 0.1803 0.3521 (-0.5098, 0.8705) 0.6086 ## CYC_1 -0.6872 0.7653 (-2.1872, 0.8128) 0.3692 ## FVC0.CYC_1 -0.0517 0.0746 (-0.1979, 0.0945) 0.4880 ## FIB0.CYC_1 -0.4665 1.1099 (-2.6419, 1.7089) 0.6743 ## FVC0_2 -0.0677 0.0271 (-0.1208,-0.0147) 0.0123 ## FIB0_2 0.1965 0.3290 (-0.4484, 0.8414) 0.5503 ## CYC_2 0.3137 0.4665 (-0.6007, 1.2280) 0.5013 ## FVC0.CYC_2 0.1051 0.0410 ( 0.0248, 0.1854) 0.0103 ## FIB0.CYC_2 0.1239 0.4120 (-0.6836, 0.9314) 0.7636 ## ## Association parameter: ## v2 1.9949 2.3093 (-2.5314, 6.5212) 0.3877 ## ## Random effects: ## sigma_b11 0.2215 0.0294 ( 0.1638, 0.2792) 0.0000 ## sigma_u 0.0501 0.0898 (-0.1259, 0.2260) 0.5772 ## Covariance: ## sigma_b1u -0.0997 0.0797 (-0.2560, 0.0565) 0.2109 The resulting table contains three parts, the fixed effects in longitudinal model, survival model and random effects. It gives the estimated parameters in the first column, standard error in the second column, 95% confidence interval and $$p$$-value for these parameters in the third and forth columns. In our example, there is only one random effect, if there are more than one random effect, the outputs will include $$sigma_b11, sigma_b12, sigma_b22, sigma_b1u, sigma_b2u$$ and so on. By the way, the supporting function coef() can be used to extract the coefficients of the longitudinal/survival process: coef(jmcfit, coeff = "beta") ## intercept time FVC0 FIB0 CYC FVC0.CYC ## 66.04146267 -0.06164756 0.90166283 -1.77799172 0.01503104 0.13798885 ## FIB0.CYC time.CYC ## 1.70883750 0.12776670 coef(jmcfit, coeff = "gamma") ## FVC0 FIB0 CYC FVC0.CYC FIB0.CYC ## [1,] 0.01871359 0.1803249 -0.6872099 -0.05172157 -0.4664724 ## [2,] -0.06772664 0.1965190 0.3136709 0.10509986 0.1239203 We proceed by checking the fit of the model using linearTest(). ## Linear hypothesis of testing all coefficients of beta's / gamma's equal 0 linearTest(jmcfit,coeff="beta") ## Chisq df Pr(>|Chi|) ## L*beta=Cb 1072.307 7 0.0000 linearTest(jmcfit,coeff="gamma") ## Chisq df Pr(>|Chi|) ## L*gamma=Cg 11.06558 10 0.3524 We can see that the hypothesis that $$\beta_1=\beta_2=\cdots=\beta_7=0$$ will be rejected, and the hypothesis $$\gamma_{11}=\gamma_{12}=\cdots= \gamma_{15}=0$$ and $$\gamma_{21}=\gamma_{22}=\cdots=\gamma_{25}=0$$ will not be rejected at the nominal level of 0.05. We extract the standard error and 95 percent confidence interval using summary(). ## Extract the standard errors for the longitudinal portion summary(jmcfit, coeff = "longitudinal") ## Longitudinal Estimate SE 95%Lower 95%Upper p-values ## 1 intercept 66.0415 0.7541 64.5634 67.5196 0.0000 ## 2 time -0.0616 0.0790 -0.2165 0.0932 0.4353 ## 3 FVC0 0.9017 0.0365 0.8302 0.9732 0.0000 ## 4 FIB0 -1.7780 0.5605 -2.8767 -0.6793 0.0015 ## 5 CYC 0.0150 0.9678 -1.8819 1.9119 0.9876 ## 6 FVC0.CYC 0.1380 0.0650 0.0106 0.2654 0.0338 ## 7 FIB0.CYC 1.7088 0.7643 0.2109 3.2068 0.0254 ## 8 time.CYC 0.1278 0.1102 -0.0883 0.3438 0.2464 ## Extract the standard errors for the survival portion summary(jmcfit, coeff = "survival") ## Survival coef exp(coef) SE(coef) 95%Lower 95%Upper p-values ## 1 FVC0_1 0.0187 1.0189 0.0326 -0.0452 0.0826 0.5660 ## 2 FIB0_1 0.1803 1.1976 0.3521 -0.5098 0.8705 0.6086 ## 3 CYC_1 -0.6872 0.5030 0.7653 -2.1872 0.8128 0.3692 ## 4 FVC0.CYC_1 -0.0517 0.9496 0.0746 -0.1979 0.0945 0.4880 ## 5 FIB0.CYC_1 -0.4665 0.6272 1.1099 -2.6419 1.7089 0.6743 ## 6 FVC0_2 -0.0677 0.9345 0.0271 -0.1208 -0.0147 0.0123 ## 7 FIB0_2 0.1965 1.2172 0.3290 -0.4484 0.8414 0.5503 ## 8 CYC_2 0.3137 1.3684 0.4665 -0.6007 1.2280 0.5013 ## 9 FVC0.CYC_2 0.1051 1.1108 0.0410 0.0248 0.1854 0.0103 ## 10 FIB0.CYC_2 0.1239 1.1319 0.4120 -0.6836 0.9314 0.7636 The implement of jmo() is the same as that of jmc(). As an illustrative example, we use the data from NINDS rt-PA stroke trial (Group 1995). In this study, 624 patients are included, and the patients treated with rt-PA were compared with those given placebo to look for an improvement from baseline in the score on the modified Rankin scale, an ordinal measure of degree of disability with categories ranging from no symptoms, no significant disability to severe disability or death, which means in this example, $$Y_{ij}$$ takes $$K=4$$ ordinal values. The following covaraites are considered: treatment group (rt-PA or placebo), modified Rankin scale prior stroke onset, time since randomization (dummy variables for 3, 6 and 12 months), and the three subtypes of acute stroke (small vessel occlusive disease, large vessel atherosclerosis or cardioembolic stroke, and unknown reasons). Similarly, we also consider the informative and noninformative risks. The model setups are as follows, $\begin{eqnarray*} P(Y_{ij}\leq k)&=&[1+\exp(-\theta_{k}-(\beta_1Group+\beta_2\mbox{Modified Rankin scale prior onset}+\beta_3time3\\ &&+\beta_4time6+\beta_5time12+\beta_6\mbox{Small vessel}+\beta_7\mbox{Large vessel or cardioembolic stroke}\\ &&+\beta_8 \mbox{Small vessel*group}+\beta_9\mbox{Large vessel or cardioembolic stroke*group})\\ &&-(\alpha_{k1}\mbox{Small vessel}+\alpha_{k2}\mbox{Large vessel or cardioembolic stroke})-b_i)]^{-1}. \end{eqnarray*}$ where $$k=1,\cdots,K-1$$. $\begin{eqnarray*} \lambda_1(t)&=&\lambda_{01}(t)\exp(\gamma_{11}Group+\gamma_{12}\mbox{Modified Rankin scale prior onset}\\ &&+\gamma_{13}\mbox{Small vessel}+\gamma_{14}\mbox{Large vessel or cardioembolic stroke}\\ &&+\gamma_{15} \mbox{Small vessel*group}+\gamma_{16}\mbox{Large vessel or cardioembolic stroke*group}+u_i)\\ \lambda_2(t)&=&\lambda_{02}(t)\exp(\gamma_{21}Group+\gamma_{22}\mbox{Modified Rankin scale prior onset}\\ &&+\gamma_{23}\mbox{Small vessel}+\gamma_{24}\mbox{Large vessel or cardioembolic stroke}\\ &&+\gamma_{25} \mbox{Small vessel*group}+\gamma_{26}\mbox{Large vessel or cardioembolic stroke*group}+\nu_2u_i), \end{eqnarray*}$ We first load the package and the data. set.seed(123) require(JMcmprsk) data(ninds) yread <- ninds[, c(1, 2:14)] cread <- ninds[, c(1, 15, 16, 6, 10:14)] cread <- unique(cread) Next, the model can be specified by the function jmo(): jmofit <- jmo(yread, cread, out = "Y", FE = c("group", "time3", "time6", "time12", "mrkprior", "smlves", "lvORcs", "smlves.group", "lvORcs.group"), cate = NULL,RE = "intercept", NP = c("smlves", "lvORcs"), ID = "ID",intcpt = 1, quad.points = 20, max.iter = 1000, quiet = FALSE, do.trace = FALSE) jmofit We can also call the internal function jmo_0() to get the same model fit. Below is the illustrative example of using the same dataset. require(JMcmprsk) data(ninds) yread <- ninds[, c(2:14)] mread <- as.data.frame(table(ninds$ID))
cread <- ninds[, c(1, 15, 16, 6, 10:14)]
jmofit

with $$p$$ the dimension of proportional odds covariates (not including intercept) in yread and $$s$$ the dimension of non-proportional odds covariates in yread. If we would like to see a concise summary of the result we can simply type

#because of the long running time, we save the jmo and jmc results within the package
fitfile=system.file("extdata", "runfit.RData", package = "JMcmprsk")
jmofit
##
## Call:
##  jmo(long_data = yread, surv_data = cread, out = "Y", FE = c("group", "time3", "time6", "time12", "mrkprior", "smlves", "lvORcs", "smlves.group", "lvORcs.group"), RE = "intercept", NP = c("smlves", "lvORcs"), ID = "ID", cate = NULL, intcpt = 1, quad.points = 20, max.iter = 1000, quiet = FALSE, do.trace = FALSE)
##
## Data Summary:
## Number of observations: 1906
## Number of groups: 587
##
## Proportion of competing risks:
## Risk 1 : 32.88 %
## Risk 2 : 4.26 %
##
## Numerical intergration:
## Number of quadrature points:  20
##
## Model Type: joint modeling of longitudinal ordinal and competing risks data
##
## Model summary:
## Longitudinal process: partial proportional odds model
## Event process: cause-specific Cox proportional hazard model with unspecified baseline hazard
##
## Loglikelihood:  -2292.271
##
## Longitudinal sub-model proportional odds:  Y ~ group + time3 + time6 + time12 + mrkprior + smlves + lvORcs + smlves.group + lvORcs.group
## Longitudinal sub-model non-proportional odds: smlves_NP + lvORcs_NP
##
##                   Estimate   Std. Error       95% CI            Pr(>|Z|)
## Longitudinal:
##  Fixed effects:
##   proportional odds:
##   group          1.6053       0.1905      ( 1.2319, 1.9786)      0.0000
##   time3          2.5132       0.1934      ( 2.1341, 2.8923)      0.0000
##   time6          2.6980       0.1962      ( 2.3134, 3.0825)      0.0000
##   time12         2.9415       0.2004      ( 2.5486, 3.3344)      0.0000
##   mrkprior       -2.1815      0.2167      (-2.6063,-1.7567)      0.0000
##   smlves         6.4358       0.4228      ( 5.6072, 7.2644)      0.0000
##   lvORcs         -1.2907      0.2861      (-1.8515,-0.7300)      0.0000
##   smlves.group   0.4903       0.7498      (-0.9793, 1.9598)      0.5132
##   lvORcs.group   -3.2277      0.4210      (-4.0528,-2.4026)      0.0000
##   Non-proportional odds:
##   smlves_NP_2    0.2725       0.4485      (-0.6066, 1.1515)      0.5435
##   lvORcs_NP_2    -0.4528      0.2466      (-0.9362, 0.0305)      0.0663
##   smlves_NP_3    1.7844       1.0613      (-0.2958, 3.8645)      0.0927
##   lvORcs_NP_3    -0.1364      0.4309      (-0.9809, 0.7081)      0.7516
##   Logit-specific intercepts:
##   theta1         -6.2336      0.1722      (-6.5712,-5.8960)      0.0000
##   theta2         -4.1911      0.1561      (-4.4971,-3.8851)      0.0000
##   theta3         3.9806       0.1896      ( 3.6091, 4.3522)      0.0000
##
## Survival sub-model fixed effects:  Surv(surv, comprisk) ~ group + mrkprior + smlves + lvORcs + smlves.group + lvORcs.group
##
##                   Estimate   Std. Error       95% CI            Pr(>|Z|)
## Survival:
##  Fixed effects:
##   group_1        -0.4630      0.2434      (-0.9400, 0.0140)      0.0571
##   mrkprior_1     0.5874       0.1371      ( 0.3187, 0.8560)      0.0000
##   smlves_1       -2.5570      0.7223      (-3.9728,-1.1413)      0.0004
##   lvORcs_1       0.5992       0.2485      ( 0.1120, 1.0863)      0.0159
##   smlves.group_1 -0.4990      1.4257      (-3.2934, 2.2955)      0.7264
##   lvORcs.group_1 1.1675       0.4692      ( 0.2479, 2.0871)      0.0128
##   group_2        0.2087       0.4834      (-0.7388, 1.1562)      0.6659
##   mrkprior_2     0.0616       0.4277      (-0.7766, 0.8998)      0.8854
##   smlves_2       0.7758       0.6217      (-0.4428, 1.9943)      0.2121
##   lvORcs_2       -0.3256      0.5120      (-1.3291, 0.6778)      0.5247
##   smlves.group_2 -0.0437      1.1573      (-2.3120, 2.2245)      0.9699
##   lvORcs.group_2 0.0991       1.0718      (-2.0015, 2.1998)      0.9263
##
## Association prameter:
##   v2             0.0101       0.1595      (-0.3025, 0.3227)      0.9496
##
## Random effects:
##   sigma_b11      55.6404       5.6560      ( 44.5547, 66.7261)      0.0000
##   sigma_u        6.6598       1.7196      ( 3.2894, 10.0303)      0.0001
##  Covariance:
##   sigma_b1u      -19.2452      0.7730      (-20.7602,-17.7302)      0.0000
##extract the parameter estimates of longitudinal proportional odds fixed effects
beta <- coef(jmofit, coeff = "beta")
beta
##        group        time3        time6       time12     mrkprior       smlves
##    1.6052525    2.5132318    2.6979598    2.9415109   -2.1815040    6.4357680
##       lvORcs smlves.group lvORcs.group
##   -1.2907438    0.4902711   -3.2276972
##extract the parameter estimates of longitudinal non-proportional odds fixed effects
alpha <- coef(jmofit, coeff = "alpha")
alpha
##      smlves_NP  lvORcs_NP
## [1,] 0.2724605 -0.4528214
## [2,] 1.7843743 -0.1363731
##extract the parameter estimates of longitudinal logit-specific intercept
theta <- coef(jmofit, coeff = "theta")
theta
## [1] -6.233618 -4.191114  3.980638
##extract the parameter estimates of survival fixed effects
gamma <- coef(jmofit, coeff = "gamma")
gamma
##           group   mrkprior     smlves     lvORcs smlves.group lvORcs.group
## [1,] -0.4629623 0.58738262 -2.5570430  0.5991513  -0.49896915   1.16751407
## [2,]  0.2087199 0.06161447  0.7757617 -0.3256335  -0.04372659   0.09914707
## Linear hypothesis of testing all coefficients of beta's / gamma's / alpha's equal 0
linearTest(jmofit,coeff="beta")
##              Chisq df Pr(>|Chi|)
## L*beta=Cb 1096.991  9 0.0000
linearTest(jmofit,coeff="gamma")
##               Chisq df Pr(>|Chi|)
## L*gamma=Cg 47.15038 12 0.0000
linearTest(jmofit,coeff="alpha")
##               Chisq df Pr(>|Chi|)
## L*alpha=Ca 8.776262  4 0.0669
## Extract the standard errors of both longitudinal and survival portions
summary(jmofit, coeff = "longitudinal")
##    Longitudinal    coef     SE 95%Lower 95%Upper p-values
## 1         group  1.6053 0.1905   1.2319   1.9786   0.0000
## 2         time3  2.5132 0.1934   2.1341   2.8923   0.0000
## 3         time6  2.6980 0.1962   2.3134   3.0825   0.0000
## 4        time12  2.9415 0.2004   2.5486   3.3344   0.0000
## 5      mrkprior -2.1815 0.2167  -2.6063  -1.7567   0.0000
## 6        smlves  6.4358 0.4228   5.6072   7.2644   0.0000
## 7        lvORcs -1.2907 0.2861  -1.8515  -0.7300   0.0000
## 8  smlves.group  0.4903 0.7498  -0.9793   1.9598   0.5132
## 9  lvORcs.group -3.2277 0.4210  -4.0528  -2.4026   0.0000
## 10  smlves_NP_2  0.2725 0.4485  -0.6066   1.1515   0.5435
## 11  lvORcs_NP_2 -0.4528 0.2466  -0.9362   0.0305   0.0663
## 12  smlves_NP_3  1.7844 1.0613  -0.2958   3.8645   0.0927
## 13  lvORcs_NP_3 -0.1364 0.4309  -0.9809   0.7081   0.7516
## 14       theta1 -6.2336 0.1722  -6.5712  -5.8960   0.0000
## 15       theta2 -4.1911 0.1561  -4.4971  -3.8851   0.0000
## 16       theta3  3.9806 0.1896   3.6091   4.3522   0.0000
summary(jmofit, coeff = "survival")
##          Survival    coef exp(coef) SE(coef) 95%Lower 95%Upper p-values
## 1         group_1 -0.4630    0.6294   0.2434  -0.9400   0.0140   0.0571
## 2      mrkprior_1  0.5874    1.7993   0.1371   0.3187   0.8560   0.0000
## 3        smlves_1 -2.5570    0.0775   0.7223  -3.9728  -1.1413   0.0004
## 4        lvORcs_1  0.5992    1.8206   0.2485   0.1120   1.0863   0.0159
## 5  smlves.group_1 -0.4990    0.6072   1.4257  -3.2934   2.2955   0.7264
## 6  lvORcs.group_1  1.1675    3.2140   0.4692   0.2479   2.0871   0.0128
## 7         group_2  0.2087    1.2321   0.4834  -0.7388   1.1562   0.6659
## 8      mrkprior_2  0.0616    1.0636   0.4277  -0.7766   0.8998   0.8854
## 9        smlves_2  0.7758    2.1722   0.6217  -0.4428   1.9943   0.2121
## 10       lvORcs_2 -0.3256    0.7221   0.5120  -1.3291   0.6778   0.5247
## 11 smlves.group_2 -0.0437    0.9572   1.1573  -2.3120   2.2245   0.9699
## 12 lvORcs.group_2  0.0991    1.1042   1.0718  -2.0015   2.1998   0.9263

# References

Elashoff, Robert M., Gang Li, and Ning Li. 2008. “A Joint Model for Longitudinal Measurements and Survival Data in the Presence of Multiple Failure Types.” Biometrics 64 (3): 762–71. https://doi.org/10.1111/j.1541-0420.2007.00952.x.

Group, Stroke Study. 1995. “Tissue Plasminogen Activator for Acute Ischemic Stroke.” N Eng J Med. 333: 1581–7.

Li, Ning, Robert M. Elashoff, Gang Li, and Jeffrey Saver. 2010. “Joint Modeling of Longitudinal Ordinal Data and Competing Risks Survival Times and Analysis of the NINDS Rt-PA Stroke Trial.” Stat. Med. 29 (5): 546–57. https://doi.org/10.1002/sim.3798.

Tashkin, Donald P, Robert Elashoff, Philip J Clements, Jonathan Goldin, Michael D Roth, Daniel E Furst, Edgar Arriola, et al. 2006. “Cyclophosphamide Versus Placebo in Scleroderma Lung Disease.” New England Journal of Medicine 354 (25): 2655–66.