1 Introduction
Joint modeling of longitudinal and survival data has drawn a lot of attention over the past two decades. Much of the research has been focused on data with a single event time and a single type of failure, usually under the assumption of independent censoring of event times (Tsiatis and Davidian 2004). However, in some situations interest lies with competing risks data, where there is more than one possible cause of an event or where the censoring is informative (Williamson et al. 2008). Typically, a standard linear mixed model or its extensions are used for the longitudinal submodel. Cause-specific hazards model with either unspecified or spline baseline hazards are studied for the competing risk submodels. Various types of random effects are assumed to account for the association between these submodels.
Despite various theoretical and methodological developments (Hickey et al. 2018b; Papageorgiou et al. 2019), there are still limited software packages to deal with specific problems in the analysis of follow-up data in clinical studies. To our knowledge, currently, there are three related CRAN R packages, namely JM (Rizopoulos 2012), joineR (Williamson et al. 2008), and lcmm (Proust-Lima et al. 2017), which support the modeling of longitudinal and survival data with competing risks.
The JM package provides support for competing risks via the
"CompRisk" option in the jointModel() function. In JM, a linear
mixed-effects submodel is modeled for the longitudinal part and a
relative risk submodel is assumed for each competing event. In the
current version (1.4-8), only the piecewise proportional hazards model,
where the log baseline hazard is approximated using B-splines, is
supported for the survival component. The joineR package fits the
joint model (Williamson et al. 2008) for joint models of longitudinal data
and competing risks using the joint() function. In their model, the
time-to-event data is modeled using a cause-specific Cox proportional
hazards regression model with time-varying covariates. The longitudinal
outcome is modeled using a linear mixed effects model. The association
is captured by a zero-mean shared latent Gaussian process. Parameters in
the model are estimated using an Expectation Maximization(EM) algorithm.
The lcmm package implements the support for competing risks joint
modeling in the Jointlcmm() function. Radically different from the
above two R packages, the lcmm package uses a less well-known
framework called the joint latent class model (Proust-Lima et al. 2014), which
assumes that dependency between the longitudinal markers and the
survival risk can be captured by a latent class structure entirely.
However, the lcmm package is mainly designed for prediction purpose
and may not be suited to evaluate specific assumptions regarding the
characteristics of the marker trajectory that are the most influential
on the event risk (Proust-Lima et al. 2014).
In all these packages, a time-independent shared random effects vector is usually assumed in modeling the longitudinal and survival data. However, they are not capable of fitting more flexible models with separate random effects in these submodels (Elashoff et al. 2008; Li et al. 2010). In many biomedical applications, sometimes, it is necessary to have a model which takes into account longitudinal ordinal outcomes for the longitudinal part. Yet, due to the complex nature of joint modeling, most of the available software does not support longitudinal ordinal variables (Armero et al. 2016; Ferrer 2017). We thus decided to fill this gap and implemented a joint model which supports ordinal disease markers based on our previous work (Li et al. 2010).
Both JM and joineR packages depend heavily on the R nlme and
survival packages. In JM, the linear mixed-effects submodel and
the survival submodel are first fitted usinglme() and coxph() R
function in these packages before a joint modeling process. In
joineR, lme() and coxph() functions are applied to obtain
initial values for parameters in the joint model, which are further
estimated by an EM algorithm. The major advantage of using available
packages such as survival and nlme lies that joint modeling R
packages can be built quickly with adequate efficiency as most of these
base R packages have been optimized for speed. However, if required
functionality is not available in these packages, as is the case of
(Elashoff et al. 2008) and (Li et al. 2010), implementing new joint modeling
methods is a non-trivial task.
Compared with JM and joineR packages, the JMcmprsk package introduced here can be regarded as a "stand-alone" R package, which does not required initial estimates for the linear mixed effects model or survival submodel to compute parameters of the joint model in question. In particular, the JMcmprsk package is built within the Rcpp (Eddelbuettel et al. 2011) and GSL(The GNU Scientific Library)(Galassi et al. 2002) framework, which make R functions have access to a wide range of fast numerical routines such as Monte Carlo integration, numerical integration and differentiation.
2 Joint Models with Competing Risks
A joint model for competing risk data consists of two linked components: the longitudinal submodel, which takes care of repeatedly measured information and the survival submodel, which deals with multiple failure times. The combination of different longitudinal and survival components leads to a variety of joint models (Hickey et al. 2018a).
In the current version of JMcmprsk, we have implemented two joint models for competing risk data, namely joint modeling with continuous longitudinal outcomes (Elashoff et al. 2008), and joint modeling with ordinal longitudinal outcomes (Li et al. 2010). Both models have adopted a cause-specific Cox submodel with a frailty term for multiple survival endpoints. The difference between these two models lies in the longitudinal part. The former model applies a linear mixed submodel for the continuous longitudinal outcome, while the latter model includes a partial proportional odds submodel for the ordinal longitudinal outcome.
Different from previous approaches (Williamson et al. 2008; Rizopoulos 2012), we assume a flexible separate random effects structure for the longitudinal submodel and the survival submodel. Furthermore, the association between both submodels is modeled by the assumption that the random effects in two submodels jointly have a multivariate normal distribution.
Model 1: Joint modeling with continuous longitudinal 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 submodels: \[\begin{aligned} 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),~~for~~k=1,\cdots,g, \end{aligned}\] where \(X_i^{(1)}(t)\), \(X_i^{(2)}(t)\) denote the covariates for the fixed-effects \(\beta\) and \(\gamma_k\), \(\tilde X_i^{(1)}(t)\) denotes the covariates 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{aligned} &&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{aligned}\]
In the E-step, we need to calculate the expected value of all the functions of \(\theta\). Since the integral over the random effects does not have a closed-form solution, an iterative numerical method has to be employed.
In JMcmprsk, the integral over time is approximated using a Gauss-Kronrod quadrature and the computation of the integral over the individual random effects is achieved using a Gauss-Hermite quadrature. The quadrature approximates the integral using a weighted sum of function values at specified points within the domain of integration; the Gaussian quadrature is based on the use of polynomial functions. A standard option here is the Gaussian quadratic rules. In the M-step, \(\Psi\) is updated by maximizing the functions obtained from the E-step.
Model 2: Joint modeling with ordinal longitudinal 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 notations have the same meaning as in Model 1.
We propose the following partial proportional odds model for \(Y_{ij}\) \[\begin{aligned} 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{aligned}\] where \(k=1,\cdots,K-1\), \(X_{ij}\) and \(\tilde X_{ij}\) are \(p\times 1\) and \(s\times 1\) vectors of covariates 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 covariates \(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{aligned} \lambda_d(t|Z_i(t),u_i)&=&\lambda_{0d}(t)\exp(Z_i(t)^\top \gamma_d+\nu_d u_i),~~for~~d=1,\cdots,g, \end{aligned}\] where the \(l\times 1\) vector \(\gamma_d\) and \(\nu_d\) are the 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{aligned} &&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{aligned}\] 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 in this case is similar to the procedure of Model 1.
3 Package structure and functionality
The R package JMcmprsk implements the above two joint models on the basis of R package Rcpp (Eddelbuettel et al. 2011) and GSL library(Galassi et al. 2002) and is hosted at CRAN. After setting the GSL environment by following the instructions in the INSTALL file from the package, we can issue the following command in the R console to install the package:
> install.packages("JMcmprsk")
There are two major functions included in the JMcmprsk package: the
function that fits continuous outcomes jmc() and the function that
fits ordinal outcomes jmo().
jmc() function
As an illustrative example of jmc(), we consider Scleroderma Lung
Study (Tashkin et al. 2006), a double-blinded, randomized
clinical trial to evaluate the effectiveness of oral cyclophosphamide
(CYC) versus placebo in the treatment of lung disease due to
scleroderma. This study consists of 158 patients and 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 covariates, baseline
%FVC (\(FVC_0\)) and baseline lung fibrosis (\(FIB_0\)) and two risks,
informative and noninformative. The model setups are as follows:
\[\begin{aligned}
\%FVC_{ij}&=&\beta_0+\beta_1t_{ij}+\beta_2FVC_{0i}+\beta_3FIB_{0i}+\beta_4CYC_i\\
&&+\beta_5FVC_{0i}\times CYC_i+\beta_6FIB_{0i}\times CYC_i+\beta_7 t_{ij}\times CYC_i+b_it_{ij}+\epsilon,
\end{aligned}\]
and the cause-specific hazards frailty models are
\[\begin{aligned}
\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{aligned}\]
We first load the package and the data.
library(JMcmprsk)
set.seed(123)
data(lung)
yread <- lung[, c(1,2:11)]
cread <- unique(lung[, c(1, 12, 13, 6:10)])The number of rows in "yread" is the total number of measurements for all subjects in the study. For "cread", 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.
Then, "yread" and "cread" are used as the longitudinal and survival
input data for the model specified by the function jmc()as shown
below:
jmcfit <- 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 = TRUE, do.trace = FALSE)where out is the name of the outcome variable in the longitudinal
sub-model, FE the list of covariates for the fixed effects in the
longitudinal sub-model, RE the types/vector of random effects in the
longitudinal sub-model, ID the column name of subject id, cate the
list of categorical variables for the fixed effects in the longitudinal
sub-model, intcpt the indicator of random intercept coded as 1 (yes,
default) or 0(no). The option quiet is used to print the progress of
function, the default is TRUE (no printing).
A concise summary of the results can be obtained using jmcfitas shown
below:
>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.2109The 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, the standard error in the second column, and 95% confidence interval and \(p\)-value for these parameters in the third and fourth columns. In our example, there is only one random effect. If there is more than one random effect, the output will include \(sigma_b11, sigma_b12, sigma_b22, sigma_b1u, sigma_b2u\), and so on.
The supporting function coef() can be used to extract the coefficients
of the longitudinal/survival process by specifying the argument coeff,
where"beta" and "gamma" denotes the longitudinal and survival
submodel fixed effects, respectively.
beta <- coef(jmcfit, coeff = "beta")
>beta
intercept time.1 FVC0 FIB0 CYC FVC0.CYC FIB0.CYC
66.04146267 -0.06164756 0.90166283 -1.77799172 0.01503104 0.13798885 1.70883750
time.CYC
0.12776670
gamma <- coef(jmcfit, coeff = "gamma")
>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.1239203The supporting function summary() can be used to extract the point
estimate, the standard error, 95%CI, and \(p\)-values of the coefficients
of both sub-models with the option coeff to specify which submodel
fixed effects one would like to extract, and digits, the number of
digits to be printed out. We proceed below to extract the fixed effects
for both submodels:
>summary(jmcfit, coeff = "longitudinal", digits = 4)
Longitudinal coef 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>summary(jmcfit, coeff = "survival", digits = 4)
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.7636We proceed to test the global hypothesis for the longitudinal and the
survival submodels using linearTest().
>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.3524The results suggest that the hypothesis \(\beta_1=\beta_2=\cdots=\beta_7=0\) is rejected, and the hypothesis \(\gamma_{11}=\gamma_{12}=\cdots= \gamma_{15} = \gamma_{21}=\gamma_{22}=\cdots=\gamma_{25}=0\) is not rejected at the significance level of 0.05.
linearTest() can also be used to test any linear hypothesis about the
coefficients for each submodel. For example, if one wants to test
\(H_0 : \beta_1 = \beta_2\) in the longitudinal submodel, then we start
with a linear contrast Lb and pass it to linearTest().
Lb <- matrix(c(1, -1, 0, 0, 0, 0, 0), ncol = length(beta)-1, nrow = 1)
>linearTest(jmcfit, coeff="beta", Lb = Lb)
Chisq df Pr(>|Chi|)
L*beta=Cb 124.8179 1 0.0000Note that we do not include intercept for linear hypotheses testing. It is seen that the hypothesis \(\beta_1 = \beta_2\) is rejected at level 0.05 in the above example.
Similarly, a linear hypotheses testing can also be done in the survival
submodel using linearTest(). For example, if we want to test
\(H_0: \gamma_{11} = \gamma_{21}\), then we start with another linear
contrast Lg and pass it to linearTest().
Lg <- matrix(c(1, 0, 0, 0, 0, -1, 0, 0, 0, 0), ncol = length(gamma), nrow = 1)
>linearTest(jmcfit, coeff="gamma", Lg = Lg)
Chisq df Pr(>|Chi|)
L*gamma=Cg 4.301511 1 0.0381It is seen that the hypothesis \(\gamma_{11} = \gamma_{21}\) is rejected at level 0.05.
For categorical variables, jmc() function will create the appropriate
dummy variables automatically as needed within the function. The
reference group in a categorical variable is specified as the one that
comes first alphabetically. Below is another example:
First, we add two categorical variables "sex" and "race" to the longitudinal data set "yread", in which "sex" is coded as "Female" or "Male", and race is coded as "Asian", "White", "Black", or "Hispanic".
#make up two categorical variables and add them into yread
set.seed(123)
sex <- sample(c("Female", "Male"), nrow(cread), replace = TRUE)
race <- sample(c("White", "Black", "Asian", "Hispanic"),
nrow(cread), replace = TRUE)
ID <- cread$ID
cate_var <- data.frame(ID, sex, race)
if (require(dplyr)) {
yread <- dplyr::left_join(yread, cate_var, by = "ID")
}Second, we rerun the model with the two added categorical variables.
# run jmc function again for yread file with two added categorical variables
res2 <- jmc(long_data = yread, surv_data = cread,
out = "FVC", cate = c("sex", "race"),
FE = c("time", "FVC0", "FIB0", "CYC", "FVC0.CYC",
"FIB0.CYC", "time.CYC"),
RE = "time", ID = "ID", intcpt = 0,
quad.points = 20, quiet = FALSE)
res2We can obtain the estimated coefficients of the longitudinal process
using coef().
> coef(res2, coeff = "beta")
intercept time FVC0 FIB0 CYC FVC0.CYC FIB0.CYC time.CYC
67.05760799 -0.07340060 0.91105151 -1.75007966 0.02269507 0.13045588 1.58807248 0.15876200
Male Black Hispanic White
-0.77110697 -0.94635182 -0.45873814 -1.19910638jmo() function
The implementation of jmo() is very similar to that of jmc(). As an
illustrative example, we use the data from (Stroke Study 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 the 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 covariates 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{aligned}
P(Y_{ij}\leq k)&=&[1+\exp(-\theta_{k}-(\beta_1Group+\beta_2\text{ Modified Rankin scale prior onset }+\beta_3time3\\
&&+\beta_4time6+\beta_5time12+\beta_6\text{ Small vessel }+\beta_7 \text{ Large vessel or cardioembolic stroke} \\
&&+\beta_8 \text{ Small vessel}*\text{group}+\beta_9\text{ Large vessel or cardioembolic stroke}*\text{group})\\
&&-(\alpha_{k1}\text{ Small vessel}+\alpha_{k2}\text{ Large vessel or cardioembolic stroke})-b_i)]^{-1},
\end{aligned}\]
where \(k=1,\cdots,K-1\).
\[\begin{aligned}
\lambda_1(t)&=&\lambda_{01}(t)\exp(\gamma_{11}Group+\gamma_{12}\text{ Modified Rankin scale prior onset}\\
&&+\gamma_{13}\text{ Small vessel}+\gamma_{14}\text{ Large vessel or cardioembolic stroke}\\
&&+\gamma_{15}\text{ Small vessel}*group+\gamma_{16}\text{ Large vessel or cardioembolic stroke} *group+u_i)\\
\lambda_2(t)&=&\lambda_{02}(t)\exp(\gamma_{21}Group+\gamma_{22}\text{ Modified Rankin scale prior onset}\\
&&+\gamma_{23}\text{ Small vessel}+\gamma_{24}\text{ Large vessel or cardioembolic stroke}\\
&&+\gamma_{25}\text{ Small vessel}*group+\gamma_{26}\text{ Large vessel or cardioembolic stroke }*group+\nu_2u_i)
\end{aligned}\]
We first load the package and the data.
library(JMcmprsk)
set.seed(123)
data(ninds)
yread <- ninds[, c(1, 2:14)]
cread <- ninds[, c(1, 15, 16, 6, 10:14)]
cread <- unique(cread)and the other arrangements are the same with those in jmc(),
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)where NP is the list of non-proportional odds covariates and FE the
list of proportional odds covariates.
To see a concise summary of the result, we can type:
>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:
Method: Standard Guass-Hermite quadrature
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.0000The usage of function coef() is similar to those in Model 1. More
specifically, coef() can extract the coefficients of non-proportional
odds fixed effects and logit-specific intercepts. For example,
alpha <- coef(jmofit, coeff = "alpha")
>alpha
smlves_NP lvORcs_NP
[1,] 0.2724605 -0.4528214
[2,] 1.7843743 -0.1363731
theta <- coef(jmofit, coeff = "theta")
> theta
[1] -6.233618 -4.191114 3.980638The usage of function summary() is the same as in Model 1. It extracts
the point estimate, standard error, 95%CI, and \(p\)-values of the
coefficients of both submodels as demonstrated below:
> 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.9263Analogous to jmcfit, linearTest() can be used to the global
hypothesis for the longitudinal and the survival submodels.
> 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.0669According to the \(p\)-values, the hypothesis \(\beta_1=\beta_2=\cdots=\beta_9=0\) is rejected, \(\gamma_{11}=\gamma_{12}=\cdots=\gamma_{16}= \gamma_{21}=\gamma_{22}=\cdots=\gamma_{26}=0\) is rejected, but \(\alpha_{11}=\alpha_{12}= \alpha_{21}=\alpha_{22}=0\) is not rejected at the significance level of 0.05.
Similarly, linearTest() can be used to test a linear hypothesis for
non-proportional odds fixed effects in the longitudinal submodel. For
example, if we want to test \(H_0 : \alpha_{11} = \alpha_{21}\), then we
can simply type:
La <- matrix(c(1, 0, -1, 0), ncol = length(alpha), nrow = 1)
> linearTest(jmofit, coeff = "alpha", La = La)
Chisq df Pr(>|Chi|)
L*alpha=Ca 1.929563 1 0.1648It is seen that the hypothesis \(\alpha_{11} = \alpha_{21}\) is not rejected at level 0.05.
Likewise, jmo() function allows for categorical variables. Moreover,
categorical variables are allowed for setting up non-proportional odds
covariates. As an illustration, here we consider the "sex" and
"race" variables and use them as two of the non-proportional odds
covariates. Below is another example:
#Create two categorical variables and add them into yread
ID <- cread$ID
set.seed(123)
sex <- sample(c("Female", "Male"), nrow(cread), replace = TRUE)
race <- sample(c("White", "Black", "Asian", "Hispanic"), nrow(cread), replace = TRUE)
cate_var <- data.frame(ID, sex, race)
if (require(dplyr)) {
yread <- dplyr::left_join(yread, cate_var, by = "ID")
}
res2 <- jmo(yread, cread, out = "Y",
FE = c("group", "time3", "time6", "time12", "mrkprior",
"smlves", "lvORcs", "smlves.group", "lvORcs.group"), cate = c("race", "sex"),
RE = "intercept", NP = c("smlves", "lvORcs", "race", "sex"), ID = "ID",intcpt = 1,
quad.points = 20, max.iter = 10000, quiet = FALSE, do.trace = FALSE)
res2
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", "race", "sex"), ID = "ID", cate = c("race", "sex"),
intcpt = 1, quad.points = 20, max.iter = 10000, 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:
Method: Standard Guass-Hermite quadrature
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: -2271.831
Longitudinal sub-model proportional odds: Y ~ group + time3 + time6 + time12 + mrkprior + smlves +
lvORcs + smlves.group + lvORcs.group + Black + Hispanic + White + Male
Longitudinal sub-model non-proportional odds: smlves_NP + lvORcs_NP + Black_NP + Hispanic_NP +
White_NP + Male_NP
Estimate Std. Error 95% CI Pr(>|Z|)
Longitudinal:
Fixed effects:
proportional odds:
group 1.1430 0.1989 ( 0.7532, 1.5328) 0.0000
time3 2.4607 0.1963 ( 2.0758, 2.8455) 0.0000
time6 2.6310 0.1986 ( 2.2416, 3.0203) 0.0000
time12 2.8717 0.2111 ( 2.4579, 3.2854) 0.0000
mrkprior -2.3329 0.1855 (-2.6965,-1.9693) 0.0000
smlves 3.9941 0.4413 ( 3.1292, 4.8589) 0.0000
lvORcs -0.9469 0.3219 (-1.5778,-0.3160) 0.0033
smlves.group -4.3940 0.7560 (-5.8758,-2.9123) 0.0000
lvORcs.group -3.6954 0.4768 (-4.6299,-2.7608) 0.0000
Black 0.8235 0.3162 ( 0.2038, 1.4433) 0.0092
Hispanic -0.0218 0.3289 (-0.6665, 0.6229) 0.9471
White 0.0523 0.3457 (-0.6253, 0.7299) 0.8797
Male -0.3528 0.2323 (-0.8080, 0.1025) 0.1288
Non-proportional odds:
smlves_NP_2 0.3314 0.4310 (-0.5133, 1.1761) 0.4419
lvORcs_NP_2 -0.3148 0.2696 (-0.8432, 0.2136) 0.2429
Black_NP_2 0.3781 0.2936 (-0.1973, 0.9535) 0.1978
Hispanic_NP_2 -0.0303 0.3176 (-0.6528, 0.5923) 0.9241
White_NP_2 -0.3802 0.3034 (-0.9748, 0.2144) 0.2102
Male_NP_2 0.0531 0.2221 (-0.3822, 0.4884) 0.8110
smlves_NP_3 2.2743 1.0748 ( 0.1677, 4.3809) 0.0343
lvORcs_NP_3 0.0033 0.4632 (-0.9045, 0.9111) 0.9943
Black_NP_3 -0.2274 0.5419 (-1.2896, 0.8349) 0.6748
Hispanic_NP_3 -0.5070 0.5087 (-1.5040, 0.4901) 0.3190
White_NP_3 0.4205 0.5722 (-0.7010, 1.5420) 0.4624
Male_NP_3 -0.8489 0.3911 (-1.6155,-0.0824) 0.0300
Logit-specific intercepts:
theta1 -6.0565 0.2868 (-6.6186,-5.4945) 0.0000
theta2 -4.0881 0.2379 (-4.5545,-3.6217) 0.0000
theta3 4.1340 0.3437 ( 3.4602, 4.8077) 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.2815 0.2545 (-0.7802, 0.2173) 0.2687
mrkprior_1 0.6404 0.1549 ( 0.3367, 0.9440) 0.0000
smlves_1 -1.8107 0.8252 (-3.4280,-0.1934) 0.0282
lvORcs_1 0.4894 0.2450 ( 0.0092, 0.9696) 0.0458
smlves.group_1 1.2608 1.6390 (-1.9517, 4.4733) 0.4417
lvORcs.group_1 1.4503 0.4901 ( 0.4898, 2.4108) 0.0031
group_2 0.2073 0.4831 (-0.7396, 1.1542) 0.6678
mrkprior_2 0.0617 0.4343 (-0.7896, 0.9129) 0.8871
smlves_2 0.7871 0.6026 (-0.3940, 1.9683) 0.1915
lvORcs_2 -0.3266 0.5085 (-1.3233, 0.6701) 0.5207
smlves.group_2 -0.0374 1.1600 (-2.3110, 2.2362) 0.9743
lvORcs.group_2 0.0952 1.0591 (-1.9807, 2.1711) 0.9284
Association prameter:
v2 0.0036 0.1577 (-0.3056, 0.3128) 0.9818
Random effects:
sigma_b11 49.0241 5.0606 ( 39.1053, 58.9430) 0.0000
sigma_u 6.3475 1.5884 ( 3.2343, 9.4607) 0.0001
Covariance:
sigma_b1u -17.6331 0.7415 (-19.0864,-16.1797) 0.0000coef(res2, coeff = "beta")
group time3 time6 time12 mrkprior smlves lvORcs
1.14302264 2.46065107 2.63095850 2.87165209 -2.33288371 3.99407491 -0.94689649
smlves.group lvORcs.group Black Hispanic White Male
-4.39403193 -3.69535020 0.82353645 -0.02181286 0.05232005 -0.35276916Older versions of jmc() and jmo()
In the previous versions of JMcmprsk, both the previous jmc() and
jmo() functions require the longitudinal input data "yfile" to be in
a specific format regarding the order of the outcome variable and the
random and fixed effects covariates. It also requires users to create an
additional "mfile" for the longitudinal data. At the suggestions of
the reviewers, in the most recent version, we focus and develop
user-friendly versions of these functions.
However, for both package consistency and user’s convenience, we still
keep older versions of these functions in the package, and rename these
functions to jmc_0() and jmo_0(), respectively. Supporting functions
ofjmo() and jmc(), such as coef(), summary(), linearTest(),
also apply to jmc_0() and jmo_0() functions.
Here, we show the usage of jmc_0() with some simulated data and the
"lung" data used in presenting jmc() functions.
If the data are provided as files, the function jmc_0() has the
following usage:
library(JMcmprsk)
yfile=system.file("extdata", "jmcsimy.txt", package = "JMcmprsk")
cfile=system.file("extdata", "jmcsimc.txt", package = "JMcmprsk")
mfile=system.file("extdata", "jmcsimm.txt", package = "JMcmprsk")
jmc_0fit = jmc_0(p=4, yfile, cfile, mfile, point=20, do.trace = FALSE)with p the dimension of fixed effects (including the 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 program prints the iteration details. Additionally,
the option type_file controls the type of data inputs.
If data frames or matrices are provided as inputs, we set the above
type_file option as type_file = FALSE in the jmc_0() function:
library(JMcmprsk)
data(lung)
lungY <- lung[, c(2:11)]
lungC <- unique(lung[, c(1, 12, 13, 6:10)])
lungC <- lungC[, -1]
## return a vector file with the number of repeated measurements as lungM
lungM <- data.frame(table(lung$ID))
lungM <- as.data.frame(lungM[, 2])
jmc_0fit2=jmc_0(p=8, lungY, lungC, lungM, point=20, do.trace = FALSE, type_file = FALSE)Computational Complexity
To understand the computational complexity of both jmc() and jmo()
models, we carried out a variety of simulations with different sample
size and different proportions of events. However, there was no clear
trend observed between the proportions of events and running times.
Hence, only one event distribution with different sample sizes are given
here for illustration purpose. According to Figures 1 and
2, we can easily see that the run time grows much faster as
sample size increases, which implies that the computational complexity
does not follow a linear order. In this case, it will limit joint models
to handling large and even moderate sample size data. To make the joint
modeling more scalable, it is necessary to carry out a novel algorithm
to reduce its computational complexity to a linear order.
jmc()
function (from 500 to 5000). Data setup: \(p\) = 4, \(n_q = 6\), 10.4%
censoring, 51.4% risk 1, and 38.2% risk 2. The run time under each
sample size was based on one random sample.
jmo()
function (from 500 to 5000). Data setup: \(p\) = 4, \(n_q = 10\), 22.4%
censoring, 57.2% risk 1, and 20.4% risk 2. The run time under each
sample size was based on one random sample.Data Simulation
A simulation can generate datasets with exact ground truth for evaluation. Hence, the simulation of longitudinal and survival data with multiple failures associated with random effects is an important measure to assess the performance of joint modeling approaches dealing with competing risks. In JMcmprsk, simulation tools are based on the data models proposed in (Elashoff et al. 2008) and (Li et al. 2010), which can be used for testing joint models with continuous and ordinal longitudinal outcomes, respectively.
The main function for simulation data continuous longitudinal outcomes
and survival data with multiple event outcomes is called SimDataC(),
which has the following usage:
SimDataC(k_val, p1_val, p1a_val, p2_val, g_val, truebeta, truegamma,
randeffect, yfn, cfn, mfn)We briefly explain some of the important options.k_val denotes the
number of subjects in study; p1_val and p1a_valdenote the dimension
of fixed effects and random effects in longitudinal measurements,
respectively; p2_val and g_val denotes the dimension of fixed
effects and number to competing risks in survival data; truebeta and
truegamma represent the true values of fixe effects in the
longitudinal and the survival submodels, respectively. randeffect sets
the true values for random effects in longitudinal and competing risks
parts, namely in the order of \(\sigma\),\(\sigma_b\),\(\nu_2\), and
\(\sigma_u\).
The following example generates the datasets used in simulation study in (Elashoff et al. 2008):
require(JMcmprsk)
set.seed(123)
yfn="jmcsimy1.txt";
cfn="jmcsimc1.txt";
mfn="jmcsimm1.txt";
k_val=200;p1_val=4;p1a_val=1; p2_val=2;g_val=2;
truebeta=c(10,-1,1.5,0.6);truegamma=c(0.8,-1,0.5,-1); randeffect=c(5,0.5,0.5,0.5);
#writing files
SimDataC(k_val, p1_val, p1a_val, p2_val, g_val,truebeta,
truegamma, randeffect, yfn, cfn, mfn)The output of function SimDataC() contains additional censoring rate
information and newly generated files names for further usage.
$`censoring_rate`
[1] 0.21
$rate1
[1] 0.45
$rate2
[1] 0.34
$yfn
[1] "jmcsimy1.txt"
$cfn
[1] "jmcsimc1.txt"
$mfn
[1] "jmcsimm1.txt"The main function for data simulation with ordinal longitudinal outcomes
and survival data with multiple event outcomes is called SimDataO(),
the usage of which is very similar to SimDataC():
SimDataO(k_val, p1_val, p1a_val, p2_val, g_val, truebeta, truetheta,
truegamma, randeffect, yfn, cfn, mfn)All options have the same meanings as in SimDataC(), while
SimDataO() has one more option truetheta, which sets the true values
of the non-proportional odds longitudinal coefficients subset.
The following example generates the datasets used in simulation study in (Li et al. 2010):
require(JMcmprsk)
set.seed(123)
yfn="jmosimy1.txt";
cfn="jmosimc1.txt";
mfn="jmosimm1.txt";
k_val=500;p1_val=3;p1a_val=1; p2_val=2;g_val=2;
truebeta=c(-1,1.5,0.8);truetheta=c(-0.5,1);truegamma=c(0.8,-1,0.5,-1); randeffect=c(1,0.5,0.5);
#writing files
SimDataO(k_val, p1_val, p1a_val, p2_val, g_val,
truebeta, truetheta, truegamma, randeffect, yfn, cfn, mfn)The output of the above function is
$`censoring_rate`
[1] 0.218
$rate1
[1] 0.414
$rate2
[1] 0.368
$yfn
[1] "jmosimy1.txt"
$cfn
[1] "jmosimc1.txt"
$mfn
[1] "jmosimm1.txt"4 Conclusions and Future Work
In this paper, we have illustrated the capabilities of package JMcmprsk for fitting joint models of time-to-event data with competing risks for two types of longitudinal data. We also present simulation tools to generate joint model datasets under different settings. Several extensions of JMcmprsk package are planned to further expand on what is currently available. First, as the integral over the random effects becomes computationally burdensome in the case of high dimensionality, Laplace approximations or other Gauss-Hermite quadrature rules would be applied to the E-M step to speed up the computation procedure. Second, with the increasing need for predictive tools for personalized medicine, dynamic predictions for the aforementioned joint models will be added. Third, other new joint models such as joint analysis for bivariate longitudinal ordinal outcomes will be included.
5 Acknowledgements
We thank the reviewers for their insightful and constructive comments that led to significant improvements in our paper. The research of Hong Wang was partly supported by the National Social Science Foundation of China (17BTJ019). The research of Gang Li was partly supported by the National Institute of Health Grants P30 CA-16042, UL1TR000124-02, and P01AT003960.
Note
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.