# Semiparametric Joint Models for Survival and Longitudinal Data with Nonparametric Multiplicative Random Effects This function applies a maximum likelihood approach to fit the semiparametric joint models of survival and normal longitudinal data. The survival model is assumed to come from a class of transformation models, including the Cox proportional hazards model and the proportional odds model as special cases. The longitudinal process is modeled by nonparametric multiplicative random effects (NMRE) model. ```r jmodelMult(fitLME, fitCOX, data, model = 1, rho = 0, timeVarY = NULL, timeVarT = NULL, control = list(), ...) ``` ## Arguments - `fitLME`: an object inheriting from class `lme` representing a fitted nonparametric multiplicative random effects model. See Details and Note and Examples . - `fitCOX`: an object inheriting from class `coxph` representing a fitted Cox proportional hazards regression model. Specifying `x = TRUE` is required in the call to `coxph()` to include the design matrix in the object fit. See Note . - `data`: a data.frame containing all the variables included in the joint modeling. See Note . - `model`: an indicator specifying the dependency between the survival and longitudinal outcomes. Default is 1. See Details . - `rho`: a nonnegative real number specifying the transformation model you would like to fit. Default is 0, i.e. the Cox proportional hazards model. See Details . - `timeVarY`: a character string indicating the time variable in the NMRE model. See Examples . - `timeVarT`: a character string indicating the time variable in the `coxph` object. Normally it is `NULL`. See Note and Examples . - `control`: a list of control values for the estimation algorithm with components: - **tol.P**: tolerance value for convergence in the parameters with default value 1e-03. See Details . - **tol.L**: tolerance value for convergence in the log-likelihood with default value 1e-06. See Details . - **max.iter**: the maximum number of EM iterations with default value 250. - **SE.method**: a character string specifying the standard error estimation method. Default is `"PRES"`. See Details and Note . - **delta**: a positive value used for numerical differentiation in the `SE.method`. Default is 1e-05 if `"PRES"` is used and 1e-03 otherwise. See Details . - **nknot**: the number of Gauss-Hermite quadrature knots used to approximate the integrals over the random effects. Under the nonparametric multiplicative random effects model, there are only one-dimensional integrations and the default for `nknot` is 11. - ``...``: additional options to be passed to the `control` argument. ## Details The `jmodelMult` function fits joint models for survival and longitudinal data. Nonparametric multiplicative random effects models (NMRE) are assumed for the longitudinal processes. With the Cox proportional hazards model and the proportional odds model as special cases, a general class of transformation models are assumed for the survival processes. The baseline hazard functions are left unspecified, i.e. no parametric forms are assumed, thus leading to semiparametric models. For detailed model formulation, please refer to Xu, Hadjipantelis and Wang (2017). The longitudinal model (NMRE) is written as $$ Y_i(t)=\mu_i(t) + \varepsilon_i(t) = b_i\times\mathbf{B}^\top(t)\boldsymbol\gamma + \varepsilon_i(t), $$ where $\mathbf{B}(t)=(B_1(t), \cdots, B_L(t))$ is a vector of B-spline basis functions and $b_i$ is a random effect $\sim\mathcal{N}(1, \sigma_b^2)$. Note that we also allow the inclusion of baseline covariates as columns of $\mathbf{B}(t)$. If `model = 1`, then the linear predictor for the survival model is expressed as $$ \eta(t) = \mathbf{W}_i^\top(t)\boldsymbol\phi + \alpha\mu_i(t), $$ indicating that the entire longitudinal process (free of error) enters the survival model as a covariate. If other values are assigned to the `model` argument, the linear predictor for the surival model is then expressed as $$ \eta(t) = \mathbf{W}_i^\top(t)\boldsymbol\phi + \alpha b_i, $$ suggesting that the survival and longitudinal models only share the same random effect. The survival model is written as $$ \Lambda(t|\eta(t)) = G\left[\int_0^t\exp\{\eta(s)\}d\Lambda_0(s)\right], $$ where $G(x) = \log(1 + \rho x) / \rho$ with $\rho \geq 0$ is the class of logarithmic transfomrations. If `rho = 0`, then $G(x) = x$, yielding the Cox proportional hazards model. If `rho = 1`, then $G(x) = \log(1 + x)$, yielding the proportional odds model. Users could assign any nonnegative real value to `rho`. An expectation-maximization (EM) algorithm is implemented to obtain parameter estimates. The convergence criterion is either of (i) $\max \{ | \boldsymbol\theta^{(t)} - \boldsymbol\theta^{(t - 1)} | / ( | \boldsymbol\theta^{(t - 1)} | + .Machine\$double.eps \times 2 ) \} < tol.P$, or (ii) $| L(\boldsymbol\theta^{(t)}) - L(\boldsymbol\theta^{(t - 1)})| / ( | L(\theta^{(t - 1)}) | + .Machine\$double.eps \times 2 ) < tol.L$, is satisfied. Here $\boldsymbol\theta^{(t)}$ and $\boldsymbol\theta^{(t-1)}$ are the vector of parameter estimates at the $t$-th and $(t-1)$-th EM iterations, respectively; $L(\boldsymbol\theta)$ is the value of the log-likelihood function evaluated at $\boldsymbol\theta$. Users could specify the tolerance values `tol.P` and `tol.L` through the `control` argument. For standard error estimation for the parameter estimates, three methods are provided, namely `"PRES"`, `"PFDS"` and `"PLFD"` (detailed information are referred to Xu, Baines and Wang (2014)). In the `control` argument, if `SE.method = "PRES"`, numerically differentiating the profile Fisher score vector with Richardson extrapolation is applied; if `SE.method = "PFDS"`, numerically differentiating the profile Fisher score vector with forward difference is applied; if `SE.method = "PLFD"`, numerially (second) differentiating the profile likelihood with forward difference is applied. Generally, numerically differentiating a function $f(x)$ (an arbitrary function) with forward difference is expressed as $$ f^\prime(x) = \frac{f(x + \delta) - f(x)}{\delta}, $$ and that with Richardson extrapolation is expressed as $$ f^\prime(x) = \frac{f(x - 2\delta) - 8f(x - \delta) + 8f(x + \delta) - f(x + 2\delta)}{12\delta}. $$ Users could specify the value of $\delta$ through the `delta` item in the `control` argument. ## Returns See `jmodelMultObject` for the components of the fit. ## References Dabrowska, D. M. and Doksun K. A. (1988) Partial Likelihood in Transformation Models with Censored Data. **Scandinavian Journal of Statistics** 15 , 1--23. Ding, J. and Wang, J. L. (2008) Modeling longitudinal data with nonparametric multiplicative random effects jointly with survival data. **Biometrics** 64 , 546--556. Tsiatis, A. A. and Davidian, M. (2004) Joint modeling of longitudinal and time-to-event data: an overview. **Statistica Sinica** 14 , 809--834. Xu, C., Baines, P. D. and Wang, J. L. (2014) Standard error estimation using the EM algorithm for the joint modeling of survival and longitudinal data. **Biostatistics** 15 , 731--744 Xu, C., Hadjipantelis, P. Z. and Wang, J. L. (2020) Semiparametric joint modeling of survival and longitudinal data: the R package JSM. **Journal of Statistical Software** . Zeng, D. and Lin, D. (2007) Maximum likelihood estimation in semiparametric regression models with censored data. **Journal of the Royal Statistical Society: Series B** 69 , 507--564. ## Author(s) Cong Xu helenxu1112@gmail.com Pantelis Z. Hadjipantelis pantelis@ucdavis.edu ## Note 1. To fit a nonparametric multiplicative random effects model, the fixed effect in the `fitLME` object should be a matrix of B-spline basis functions (an object from the `bs` function) with the possibility of including baseline covariates and the random effect should only include a random intercept. In the `bs` function, it is a good practice to specify the boundary knots through the `Boundary.knots` argument, where the upper boundary knot is typically the longest follow-up time among all subjects. See Examples . 2. Currently, `jmodelMult()` could only handle the `fitLME` object with a simple random-effects structure (only the `pdDiag()` class). Moreover, the within-group correlation and heteroscedasticity structures in the `fitLME` object (i.e. the `correlation` and `weights` argument of `lme()`) are ignored. 3. The `data` argument in `jmodelMult()`, `lme()` and `coxph()` should be the same data frame. 4. For the `fitCOX` object, only the $W_i(t)$ in the linear predictor $\eta(t)$ for the survial model (see Details ) should be involved in the `formula` argument of `coxph{}`. Since `coxph()` uses the same data frame as `lme()` does, a time-dependent Cox model must be fitted by `coxph()` although $W_i(t)$ may only contain time-independent covariates. See Examples . 5. If $W_i(t)$ in the linear predictor $\eta(t)$ for the survial model (see Details ) does involve time-dependent covariate, then `timeVarT` must specify the name of the time variable involved (see Examples ). 6. The standard error estimates are obtained by numerical approximations which is naturally subject to numerical errors. Therefore, in extreme cases, there may be `NA` values for part of the standard error estimates. ## See Also `jmodelMultObject`, `lme`, `coxph`, `Surv`, `bs` ## Examples ```r # linear mixed-effects model fit where the fixed effect is modeled by # quadratic B-splie basis with no internal knots fitLME <- lme(log(serBilir) ~ bs(obstime, degree = 2, Boundary.knots = c(0, 15)), random = ~ 1 | ID, data = pbc) # Cox proportional hazards model fit with a single time-independent covariate fitCOX <- coxph(Surv(start, stop, event) ~ drug, data = pbc, x = TRUE) # joint model fit which assumes the Cox proportional hazards model for the survival process # and NMRE for the longitudinal process. Use 'max.iter = 25', 'nknot = 3' and # the 'PFDS' method to calculate standard error estimates as a quick toy example fitJTMult.ph <- jmodelMult(fitLME, fitCOX, pbc, timeVarY = "obstime", control = list(SE.method = 'PFDS', max.iter = 25, nknot = 3)) summary(fitJTMult.ph) ## Not run: # joint model fit with the default control fitJTMult.ph2 <- jmodelMult(fitLME, fitCOX, pbc, timeVarY = "obstime") summary(fitJTMult.ph2) # joint model fit where the survival and longitudinal processes only share # the same random effect fitJTMult.ph3 <- jmodelMult(fitLME, fitCOX, pbc, model = 2, timeVarY = "obstime") summary(fitJTMult.ph3) # joint model fit which assumes the proportional odds model for the survival process # and NMRE model for the longitudinal process fitJTMult.po <- jmodelMult(fitLME, fitCOX, pbc, rho = 1, timeVarY = "obstime") summary(fitJTMult.po) # joint model fit where the survival and longitudinal processes only share # the same random effect fitJTMult.po2 <- jmodelMult(fitLME, fitCOX, pbc, model = 2, rho = 1, timeVarY = "obstime") summary(fitJTMult.po2) # allow baseline covariates in the NMRE model for the longitudinal process fitLME2 <- lme(log(serBilir) ~ drug + bs(obstime, degree = 2, Boundary.knots = c(0, 15)), random = ~1 | ID, data = pbc) fitJTMult.ph4 <- jmodelMult(fitLME2, fitCOX, pbc, timeVarY = "obstime") summary(fitJTMult.ph4) # Cox proportional hazards model fit with a time-dependent covariate fitCOX2 <- coxph(Surv(start, stop, event) ~ drug + as.numeric(drug) : obstime, data = pbc, x = TRUE) # joint model fit in which \code{timeVarT} must be specified fitJTMult.ph5 <- jmodelMult(fitLME, fitCOX2, pbc, timeVarY = "obstime", timeVarT = 'obstime', control = list(max.iter = 300)) summary(fitJTMult.ph5) ## End(Not run) ```

jmodelMult function