JMFit: A SAS Macro for Joint Models of Longitudinal and Survival Data

Joint models for longitudinal and survival data now have a long history of being used in clinical trials or other studies in which the goal is to assess a treatment effect while accounting for a longitudinal biomarker such as patient-reported outcomes or immune responses. Although software has been developed for fitting the joint model, no software packages are currently available for simultaneously fitting the joint model and assessing the fit of the longitudinal component and the survival component of the model separately as well as the contribution of the longitudinal data to the fit of the survival model. To fulfill this need, we develop a SAS macro, called JMFit. JMFit implements a variety of popular joint models and provides several model assessment measures including the decomposition of AIC and BIC as well as ΔAIC and ΔBIC recently developed in Zhang et al. (2014). Examples with real and simulated data are provided to illustrate the use of JMFit.


Introduction
The joint analysis of longitudinal and time-to-event outcomes has been widely published in statistical journals. One popular approach in joint modeling of longitudinal and survival data is based on shared random effects, where the longitudinal model and survival model share common random effects and these random effects then induce correlation between the longitudinal and survival components of the model. This family of joint models is also called the "shared parameter models" (SPMs). There are two basic formulations of SPMs. The first is the "time trajectory model", denoted by SPM1, where one essentially substitutes the polynomial time trajectory function from the longitudinal model into the hazard function of the survival model, and in this case, the trajectory function acts like a time-varying covariate in the survival model. The second formulation, denoted by SPM2, is to directly include the random effects as covariates in the survival model. There are several R packages available in fitting joint models based on shared random effects, including JM (Rizopoulos 2012), JMbayes (Rizopoulos 2014), and joineR (Philipson et al. 2012). There is also a Stata module stjm (Crowther 2012;Crowther et al. 2013), which estimates shared random effects models. In addition, another R package, lcmm (Proust-Lima et al. 2014), estimates joint models based on shared latent classes.
One important issue in the joint modeling of longitudinal and survival data concerns the separate contribution of the model components to the overall goodness-of-fit of the joint model. Recently, Zhang et al. (2014) derived a novel decomposition of the AIC and BIC criteria into additive components that will allow us to assess the goodness of fit for each component of the joint model. Such a decomposition leads to the development of ΔAIC and ΔBIC, which quantify the change of AIC and BIC in fitting the survival data with and without using the longitudinal data. Thus, ΔAIC and ΔBIC can be used to determine the importance of the longitudinal data relative to the model fit of the survival data. In addition, ΔAIC and ΔBIC are also very useful in assessing whether a linear trajectory or quadratic trajectory is more suitable and also facilitating a direct comparison between SPM1's and SPM2's. These measures will help the data analyst in not only assessing each component of the joint model but also in determining the contribution of the longitudinal measures to the fit of the survival data. These newly developed model assessment criteria are not available in any of these packages or module mentioned before. We mention here that the methodology for ΔAIC and ΔBIC was fully developed in Zhang et al. (2014), but our goal here is the novel implementation of this methodology into user-friendly software along with a class of joint models for jointly analyzing longitudinal and time-to-event data. This paper introduces JMFit, a SAS macro, that will allow us to fit the SPM1, SPM2, timevarying covariates, and two-stage models as well as to assess the goodness-of-fit of each of the longitudinal and survival components in the joint model. A detailed analysis of the longitudinal and survival data from a cancer clinical trial as well as an analysis of the simulated data are carried out to illustrate the functionality of JMFit. A detailed description of JMFit is given in Appendix A.
is right-censored for the ith subject. We further let x i (t) and z i denote a p L -dimensional vector of time-dependent covariates and a p S -dimensional vector of baseline covariates, respectively. The joint model for (y i , t i ) consists of the longitudinal component and the survival component.
For the longitudinal component, a mixed effects regression model is assumed for y i (t), which takes the form: (1) where is a polynomial vector of order q for j = 1, . . . , m i , θ i is a (q + 1)-dimensional vector of random effects, and γ is a p-dimensional vector of regression coefficients. In (1), we further assume θ i ~ N(θ, Ω), where θ is the (q +1)dimensional vector of the overall effects, Ω is a (q + 1) × (q + 1) positive definite covariance matrix with the lower triangle consisting of {Ω 00 , Ω 10 , Ω 11 , . . . , Ω qq }, ε i (a ij ) ~ N(0, σ 2 ), and θ i and ε i (a ij ) are independent. We note that in (1), if q = 1, and represents a linear trajectory, and if q = 2, and leads to a quadratic trajectory.
For the failure time t i , we assume that the hazard function takes the form (2) or (3) where λ 0 (t) is the baseline hazard function, β is a one-dimensional regression coefficient in (2) and β is a (q + 1)-dimensional vector of the regression coefficients in (3). Note that in (2) or (3), θ i and g(t) are the parameters and the functions from the longitudinal component of the joint model in (1) while λ 0 , β (or β), and α are the only parameters pertaining to the survival component. As shown in Section 3.1, β (or β) controls the association between the longitudinal marker and the time-to-event. A value of β = 0 (or β = 0) implies no association between the longitudinal marker and the time-to-event. The joint model with hazard function specified in (2) is the trajectory model, denoted by SPM1, while the one with hazard function given by (3) is denoted by SPM2. Under SPM1, a positive value of β implies that a

The construction of the piecewise constant baseline hazard function
Assuming λ 0 (t) to be a piecewise constant baseline hazard function, we partition the time axis into J intervals with . Then we assign a constant baseline hazard to each of the J intervals, that is, Let be the n* event times of the t i 's, where . We consider four algorithms to construct the 's.
Step 2: Let n j = [p j n*], which is the integer part of p j n*.
The ESQP is a popular approach to construct the piecewise constant hazard function, which is discussed in Ibrahim et al. (2001, Chapter 5) and also implemented in the R package JM developed by Rizopoulos (2010). We note that is the p j th quantile of the 's and (5) is implemented in the SAS UNIVARIATE procedure as the default option for computing quantiles. We also note that the ESQP algorithm does not yield nested partitions in the sense that is not necessarily a subset of when J 1 < J 2 . In order to construct partitions, we propose the following three bi-sectional quantile partition algorithms.
Step 4: Use Steps 2 and 3 of Algorithm 1 to compute .

Steps 3 and 4:
The same as Algorithm 2.
Steps 3 and 4: The same as Algorithm 2.
If there are ties in , say , the interval ( , ] is undefined. Thus we only use distinct values of the . Let J t denote the number of ties in . Then the number of distinct intervals reduces to J − J t . Figure 1 shows how the partition intervals are constructed based on LBSQP, MBSQP, and RBSQP. Notice that when J = 2 K , K = 1, 2, . . . , ESQP, LBSQP, MBSQP, and RBSQP yield the same partition. LBSQP, MBSQP, and RBSQP are desirable when there are more events at the beginning, in the middle, and at the end of the follow-up period, respectively. Another advantage of LBSQP, MBSQP, and RBSQP is that the resulting partitions are nested and, hence, the log-likelihood of the joint model increases in J when the longitudinal component remains fixed.

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript and L 0 (λ, α|t i , δ i , z i ) is defined by (11). The ΔAIC or ΔBIC measure the gain of the fit in the survival component due to the longitudinal data with a penalty for the additional parameters in the survival component of the joint model. The model with a large value of ΔAIC (ΔBIC) is more preferred.

Design
The SAS macro JMFit has been developed to assess model fit in joint models of longitudinal and survival data. In fact, it can fit five models, including the two types of joint models with linear and quadratic trajectories, as well as the time-varying covariates model. The macro JMFit consists of five submacros, SPM1L, SPM1Q, SPM2L, SPM2Q, and TVC, corresponding to the five models, respectively. The MODEL argument of JMFit specifies one of the following five models to be fitted: SPM1L: SPM1 with Linear trajectory. The hazard function has the form SPM1Q: SPM1 with Quadratic trajectory. The hazard function has the form SPM2L: SPM2 with Linear trajectory. The hazard function has the form where θ 0i and θ 1i are subject-level random intercept and random slope.
TVC: Time-Varying Covariates model. The hazard function has the form where y i (t) = y i (a ij ) for a ij ≤ t < a i,j+1 for j = 1, . . . , m i , where a i,m i +1 = ∞.
Note that the TVC model is a"non-joint"model, and the use of this model has great potential for bias (Fisher and Lin 1999). Zhang et al. Page 8 We provide the two versions of SPM1L, SPM1Q, SPM2L, and SPM2Q with the TS argument. If TS is missing or equal to 0, the joint model will be fit; while "TS=1" yields the corresponding two-stage model. Similar to the method in Tsiatis et al. (1995), (i) we first fit the linear mixed model specified in (1) to the longitudinal data alone and then obtain the estimates of θ i , denoted by ; and (ii) we replace θ i in (2) or (3) with the estimate at the second stage. The only difference between the joint model and the corresponding two-stage model is that θ i in (2) or (3) is replaced with the estimate . This two-stage approach may potentially lead to biased and inefficient estimates (Ibrahim et al. 2010).

Implementation details
If the observed longitudinal measures are sparse, the full trajectories of longitudinal measures might not be well estimated. For example, in the case of fitting a quadratic trajectory, the sign of the estimated second-order coefficient could be incorrect if the longitudinal measures were observed only within the first half of the follow-up period, leading to incorrect extrapolation when the observed progression time was far beyond the time of the last observed longitudinal measure.
Let t max,i = max 1≤jm i {a ij }. When t > t max,i , y i (t) is never observed and no longitudinal data are available to estimate the trajectory for t > t max,i . Under SPM1, the extrapolation of the trajectory beyond t max,i may lead to a survival component of the joint model that fits the survival data poorly. In addition, such an extrapolation also causes a severe convergence problem in the SAS NLMIXED (SAS Institute Inc. 2011b) procedure especially when t max,i < < t i for many subjects. To circumvent these issues, for SPM1, we modify the hazard function in (2) as (14) or (15) where , , w(∈ [0,1]) is the proportion of max(t i − t max,i ,0), and τ = max 1≤i≤n {t i }, which is the last follow-up survival time. If w = 0 (default), t max,i will be the starting point of the modified extrapolation of the trajectory; while w = 1 implies that the trajectory extends to t i with no t max,i adjustment. This modification can also be applied to the TS model corresponding to SPM1L or SPM1Q. There are two arguments called TMAXI and WEIGHT. If TMAXI is missing or equal to 0, no t max,i adjustment will be applied. If TMAXI=1, the t max,i adjustment based on hazard function given in (14) will be applied with weight given in the WEIGHT argument, that is, the trajectory will become flat after . If TMAXI=2, the t max,i adjustment based on hazard function given in (15) will be applied with weight given in the WEIGHT argument, that is, starting at , the trajectory will linearly go down to 0 at the last follow-up survival time. An "optimal" choice of weight w in conjunction with TMAXI option may be determined by either AIC Surv|Long or BIC Surv|Long . The purpose of defining t max,i is that we create an extrapolation of the longitudinal measures so that the trajectory function can be well estimated. This extrapolation is needed when there are few longitudinal measures at later points. We note here that the t max,i method corresponds to a "Prediction Carried Forward" approach, which may potentially induce bias in the estimation. We also note that this is not an issue for the SPM2, in which the hazard function is independent of g(t).
The OPTIONS argument allows users to specify options (e.g., integration method, optimization technique, and convergence criteria) that are available in the PROC NLMIXED statement. If OPTIONS is missing, JMFit will use adaptive Gaussian quadrature to approximate the integral of the likelihood over the random effects, perform a quasi-Newton optimization, and apply a relative gradient convergence criterion of 10 −8 . All the methods mentioned above are default methods in PROC NLMIXED. The Riemann integral is used to compute the cumulative hazard function for the trajectory models. Each time interval is divided into 200 subintervals.
A big challenge in fitting joint models using the SAS NLMIXED procedure is convergence.
Poor initial values may lead to the failure of convergence in NLMIXED. To address this issue, we first fit the longitudinal data alone using the SAS MIXED procedure to obtain the estimates, , , , , , for the parameters in the longitudinal component of the joint model.

The EMPHACIS data
To illustrate how JMFit works, we use the data from a multicenter, randomized, single-blind, EMPHACIS lung cancer clinical trial (Evaluation of MTA in Mesothelioma in a Phase III Study with Cisplatin). The study drug was a multi-targeted antifolate (MTA) pemetrexed which was given in combination with cisplatin (the PEM/Cis arm), and the active-treatment comparator was cisplatin alone (the Cis arm). The treatment for both arms was structured as six 21-day cycles of therapy; patients receiving treatment benefit could receive additional cycles based on investigator discretion.
Malignant pleural mesothelioma is a rapidly progressing and highly symptomatic malignancy with a median survival time of 6 to 9 months. Accordingly, patient-reported assessments are important for evaluation of disease progression and patients' response to therapy. In oncology, the patients' importance ratings on the magnitude of progression-free survival improvement have been shown to depend on the severity of disease-related symptoms (Bridges et al. 2012). We analyzed the disease-specific patient-reported Lung Cancer Symptom Scales (LCSS) to evaluate the patient-level association of five instrument items, anorexia, cough, dyspnea, fatigue, and pain, with progression-free survival using EMPHACIS trial data.
We consider the same subset of the EMPHACIS data as in Zhang et al. (2014), which consists of 425 patients with at least one post-baseline value of each longitudinal measure and seven binary covariates, including therapy (1=pemetrexed/cisplatin and 0=cisplatin alone), race (1=white and 0=others), gender (1=male and 0=female), age (1='≥ 65' and 0='< 65'), Karnofsky (1='90-100' and 0='<90'), stage (1=stage I/II and 0=stage III/IV), and bf (1=full vitamin supplementation and 0=other). The detailed description of this study can be found in Zhang et al. (2014). By applying joint models in this study, we compare the longitudinal LCSS symptoms in terms of their contribution to the overall fit of survival data via ΔAIC and ΔBIC, which are computed using JMFit.
Suppose we fit the trajectory model with a linear trajectory and J = 3 for pain. Then we create two data sets named as "lcss_long_pain" and "lcss_surv_pain" for the macro JMFit's LONG and SURV options, respectively. The MODEL option is set to "SPM1L" and TS is set to 0. We assign 1 to TMAXI for t max,i adjustment with WEIGHT=0. The NPIECES option is set to 3 and the PARTITION option is set to 2 for the LBSQP algorithm. The rtf file "Output for pain under SPM1L with J = 3 (Partition=2).rtf" is generated by JMFit, which lists five tables. Table "Number of Subjects" given in Figure 2 shows that there are 425 patients in data sets lcss_long_pain and lcss_surv_pain, respectively, and 425 subjects are used, implying that the id's of the subjects in these two data sets match. Figure 3 shows the"Fit Statistics"table, which consists of the log likelihood, AIC Long (BIC Long ), AIC Surv (BIC Surv ), and ΔAIC (ΔBIC). Inside the JMFit macro, PROC NLMIXED provides the log likelihood, AIC and BIC, and PROC IML (SAS Institute Inc. 2011a) is used to compute AIC Long and BIC Long .
Next, the table, titled "Survival Parameter Estimates (Survival Alone)" shown in Figure 4, is obtained by fitting the survival data alone. For this table, the estimate, the standard error (SE), the degrees of freedom (DF), the t-value, the p-value, and the 95% confidence interval (CI) are all shown for each parameter. In addition, the gradient of the negative log-likelihood function is displayed, which can be used to check convergence of PROC NLMIXED. A small gradient implies better convergence. The largest absolute value of the gradients shown in Figure 4 is 0.00322, indicating good convergence of PROC NLMIXED.
PROC NLMIXED produces the table titled "Parameter Estimates" in Figure 5, which consists of three subtables: "Covariance Parameter Estimates", "Longitudinal Parameter Estimates", and "Survival Parameter Estimates". For each parameter in this table, the estimate, SE, DF, t value, p value, 95% CI, and gradient are all provided. The Subtable "Covariance Parameter Estimates" consists of the estimates of the three lower-triangle elements (Ω 00 , Ω 10 , Ω 11 ) of the random-effects covariance matrix as well as the estimate of the standard deviation (σ) for the error term. The estimates of the coefficients of the variables from the longitudinal component are shown in Subtable "Longitudinal Parameter Estimates". In Subtable "Survival Parameter Estimates", "Param/Var" lists all the names of the variables as well as the parameters from the survival component. In this subtable, log λ 1 , log λ 2 , and log λ 3 are the natural logarithms of the piecewise baseline hazards for the three intervals, respectively. The hazard ratios of the covariates and β as well as the estimates of λ are shown in Figure 6. From Figure 5, we see that therapy has p-values of 0.1451 and <.0001 in the longitudinal and survival submodels, respectively, indicating that therapy is not statistically significant at the 0.05 level in the longitudinal submodel and it is highly significant in the survival submodel. Also, from Figure 6, the estimated hazard ratio of therapy in the survival model is 0.6272 with a 95% CI of (0.508, 0.775), implying that the treatment reduces the risk of disease progression by about 37%. From Figure 5, we also see that (i) karnofsky is highly significant in both the longitudinal and survival submodels and (ii) the other significant variables in the survival submodel include the coefficient β corresponding to the linear trajectory and stage. A significant coefficient β indicates that pain is highly associated with disease progression. From the "Fit Statistics" tables in Figure 3 and Figure 7, we see that pain had the largest values of ΔAIC and ΔBIC and cough had the smallest values of ΔAIC and ΔBIC, which indicate that pain led to the most gain in fitting the PFS data while cough had the least contribution to the fit of the PFS data.
Since ΔAIC = AIC Surv,0 − AIC Surv|Long , both AIC Surv,0 and AIC Surv|Long change when J changes. Thus, the "best J" based on ΔAIC may not necessarily lead to the best survival submodel or the best joint model. As an illustration, Table 1 shows the values of AIC Surv|Long and ΔAIC for pain under SPM1L with different partition algorithms and different J's. From Table 1, we see that (i) J = 7 (or J = 6) is the "best" choice according to ΔAIC; (ii) J = 2 is the "best" one based on AIC Surv|Long ; and (iii) the values of AIC Surv|Long are 2175.29, 2175.29, and 2175.60 for LBSQP, MBSQP, and RBSQP when J = 7, which are much larger than those when J = 2. Since our goal is to fit the joint model to both the longitudinal data and the survival data, it is more appropriate to use AIC Surv|Long to determine the number of intervals for the piecewise constant baseline hazard function. Table 2 shows the values of AIC Surv|Long for cough and pain under SPM1L with t max,i adjustments using different w's. We can see from Table 2 that, among the five values of w, for both t max,i adjustments, cough and pain have the smallest values of AIC Surv|Long at w = 0.25 and 0, respectively. For cough, the two t max,i adjustments result in similar values of AIC Surv|Long ; and for pain, TMAXI=2 slightly outperforms TMAXI=1. In addition, both cough and pain have the largest values of AIC Surv|Long at w = 1, indicating that the full trajectories were not well estimated.

A simulated data example
Since we do not have permission to distribute the EMPHACIS dataset used in Section 4.1 publicly, we consider a simulated data example, in which the simulated data can be downloaded directly from the journal's website.
The failure time and censoring indicator were calculated as and , where 1{A} denotes the indicator function such that 1{A} = 1 if A is true and 0 otherwise. This survival dataset is denoted by D Surv .
We also generated three additional sets of longitudinal data, with longitudinal trajectories simulated from where (i) ; (ii) ; and (iii) . Then the longitudinal data were generated from a N(μ ℓi (a ij ), σ 2 ) distribution with σ = 0.54 for ℓ = 1, 2, 3. These three additional sets of longitudinal data were coupled with the same survival times as in D Long + D Surv to form three additional data sets. These resulting data sets are denoted by D Long1 + D Surv , D Long2 + D Surv , and D Long3 + D Surv . This simulation setting is similar to Simulation III in Zhang et al. (2014) except for the six additional covariates. Figure 8 shows the fit statistics for D Long + D Surv , D Long1 + D Surv , D Long2 + D Surv , and D Long3 + D Surv , respectively, using JMFit. We see from Figure 8 that D Long + D Surv has the largest values of ΔAIC and ΔBIC.
The rest of the output for D Long + D Surv are provided in Appendix B. The output for D Longℓ +D Surv (ℓ = 1, 2, 3) are omitted here for brevity.

Concluding remarks
The JMFit SAS macro fits the joint models for longitudinal and survival data. The piecewise exponential constant hazard model is assumed for the baseline hazard function. The timeaxis is partitioned into J intervals, which are constructed by four algorithms, namely, ESQP, LBSQP, MBSQP, and RBSQP. JMFit allows users to specify the number of intervals J and the partition method. Five models, including SPM1L, SPM1Q, SPM2L, SPM2Q, and TVC, are implemented in JMFit. This SAS macro also computes AIC, BIC, ΔAIC, ΔBIC, and the estimates of the parameters in the joint model. The computational time of JMFit depends on which of SPM1L, SPM1Q, SPM2L, SPM2Q, or TVC is chosen and how big the dataset is.
For the example given in Section 4.1, it took 42 seconds to fit SPM1L with J = 3 based on LBSQP for pain on a Dell PC with an Intel i5 processor, 3.30 GHz CPU, and 8 GB of memory. On the same PC, it only took 10 to 11 seconds to fit SPM2L to each of the simulated datasets illustrated in Section 4.2.
The JMFit SAS macro provides two versions for SPM1L, SPM1Q, SPM2L, and SPM2Q with the TS argument, with "TS=1" yielding the corresponding two-stage model. Comparing these models with TS missing or equal to 0, for the two-stage models, (i) we first fit (1) to the longitudinal data alone and obtain the estimates of θ i , denoted by ; and (ii) we then use in (2) and (3). The TS models are also substantially different than the TVC model since the TS models do not require the LOCF assumption.
The current version of the JMFit SAS macro only fits linear and quadratic models for the longitudinal outcome and the piecewise constant baseline hazard function for the survival submodel. In the joint modeling framework, other dependence structures, such as dependence through the derivatives of the trajectory function or interactions with covariates as well as spline approximations to the baseline hazard could be assumed. In addition, other trajectory functions may be more appropriate to model the time effect on the longitudinal outcomes in certain applications. These additional features could be built in the JMFit macro in a future release.
Finally, we note that the EMPHACIS dataset used in this paper is proprietary and cannot be distributed publicly. However, the simulated datasets D Surv , D Long , D Long1 , D Long2 , and D Long3 are available for downloading from the journal website. be placed after column 3, which can be enforced in SAS by using the "retain" command. Note that X L1 , ..., X Lp can be time-dependent or baseline covariates. Required.
SURV: Data set with the first three columns, SID, survival time (T), censoring indicator (delta) (1 = death and 0 = censored), and additional columns for covariates (X S1 -X Sq ), where SID, T, and delta should be arranged in the first, second, and third columns, and X S1 , . . . , X Sq should be placed after column 3. Required.
MODEL: Model specification. Required. One of
TS: Indicates whether to implement the model specified in the MODEL argument or the corresponding two-stage model. If 0 (default), the model specified in the MODEL argument will be fit. If 1, the corresponding two-stage model will be fit instead. It only works for SPM1L, SPM1Q, SPM2L, and SPM2Q.
TMAXI: t max,i adjustment to the model specified in the MODEL argument. If 0 (default), no t max,i adjustment will be applied. If 1, the trajectory will become flat after . If 2, starting at , the trajectory will linearly go down to 0 at the last follow-up survival time. It only works for SPM1L and SPM1Q.
WEIGHT: The proportion (∈ [0, 1]) of max(t i − t max,i , 0). If 0 (default), the starting point of the modified extrapolation of the trajectory is t max,i . If 1, the trajectory extends to t i with no t max,i adjustment. It only works when TMAXI=1 or TMAXI=2. NPIECES: Number of intervals J (≥ 1) for the piecewise constant baseline hazard function.

Required.
PARTITION: Algorithm for constructing the partition of the time axis. Required. One of (i) 1: Equally-Spaced Quantile Partition (ESQP).
OPTIONS: Allows users to specify options that are available in the PROC NLMIXED statement. For example, "OPTIONS=%str(QPOINTS=5 TECH=CONGRA ABSGCONV=0.0001)" specifies Gaussian quadrature with five quadrature points for approximating the integral of the likelihood over the random effects, the conjugate-gradient optimization, and an absolute gradient convergence criterion of 0.0001.
INITIAL: Allows users to set their own initial values. JMFit will automatically generate the starting values for the model parameters and these initial values will be stored in the data set "_initial". Since the order of the parameters is very important when calculating AIC Long and BIC Long , users are recommended to change the initial values in "_initial" and then rename "_initial".
OUTPUT: Name of the output rich text file (rtf). One can also specify the directory in which the file will be put. For example, the output file named "myoutput" will be stored in "C: \ . . . " by "%JMFit(..., OUTPUT = C: \. . . \myoutput);" If OUTPUT is not specified, the file will be indexed by the name of Y from LONG and the model's name.
Note 1: (i) The name of the SID variable in LONG should be the same as that of the SID variable in SURV; (ii) A and T should be in the same unit of time (month preferred); (iii) the categorical covariates must be coded as dummy variables; and (iv) the SAS macro allows for any numbers of covariates for both components of the joint model and the covariates for the longitudinal component can be totally different from those for the survival component.
Note 2: (i) "ERROR: Not enough memory to generate code." This might occur if J is too big; (ii) a too long path for the OUTPUT may lead to an error; (iii) the macro is assuming "options validvarname=v7;" for valid variable names that can be processed in SAS; and (iv) the calculations of AIC Long and BIC Long require the IML Procedure. The number of subjects for D Long + D Surv .  The hazard ratios and λ estimates for D Long + D Surv under SPM2L.   An illustration of three bi-sectional partition methods.