lsjm | R Documentation |
This function fits complex joint models with shared random effects. The longitudinal submodel estimates longitudinal data with a mixed-effects model in which we suppose that the variance of the residual error is time-dependent and subject-specific. The survival submodel handles right-censored and left-truncated time-to-event data and competing risks. The dependence structure between the longitudinal and the survival data can be the random effects from the mixed model or the current value of the marker and/or the slope of the marker. We can also adjust on the current variance of the marker. (See below) Parameters are estimated simultaneously through a maximum likelihood method, using a Marquardt-Levenberg algorithm.
lsjm(
formFixed,
formRandom,
formGroup,
formSurv,
timeVar,
data.long,
variability_hetero = TRUE,
formFixedVar,
formRandomVar,
correlated_re = FALSE,
sharedtype = c("current value", "variability"),
hazard_baseline = "Exponential",
formSlopeFixed = NULL,
formSlopeRandom = NULL,
indices_beta_slope = NULL,
nb_pointsGK = 15,
ord.splines = 3,
competing_risk = FALSE,
formSurv_CR = NULL,
hazard_baseline_CR = "Exponential",
sharedtype_CR = c("current value", "variability"),
left_trunc = FALSE,
Time.0 = NULL,
S1 = 1000,
S2 = 5000,
nproc = 1,
clustertype = "SOCK",
maxiter = 100,
print.info = FALSE,
file = NULL,
epsa = 0.001,
epsb = 0.001,
epsd = 0.001,
binit = NULL,
Comp.Rcpp = TRUE
)
formFixed |
A formula for the fixed effects of the longitudinal submodel |
formRandom |
A formula for the random effects of the longitudinal submodel |
formGroup |
A formula which indicates the group variable |
formSurv |
A formula which indicates the variables used in the survival submodel |
timeVar |
The name of the column of time in data.long. This variable must appears in data.long |
data.long |
A dataframe with the longitudinal data |
variability_hetero |
A logical to indicate if we suppose a subject_specific variability |
formFixedVar |
A formula for the fixed effects of the variance predictor |
formRandomVar |
A formula for the random effects of the variance predictor |
correlated_re |
A logical to indicate if the random effects of the marker and the variance predictors are correlated (By default there are supposed to be independent) |
sharedtype |
char : dependence structure for survival model : "RE" (random effects) or "CV" (current value) or "CVS" (current value and slope) or "S" (slope) |
hazard_baseline |
char : baseline hazard function : "Exponential" or "Weibull" or "Splines" |
formSlopeFixed |
A formula for the fixed effects of the slope of the longitudinal submodel : the derivative of the formFixed |
formSlopeRandom |
A formula for the random effects of the slope of the longitudinal submodel : the derivative of the formRandom |
indices_beta_slope |
A vector of index indicating which beta of the formFixed formula is used in the formSlopeFixed formula |
nb_pointsGK |
the number of points for Gauss-Kronrod approximation : choice between 7 and 15. 15 by default. |
ord.splines |
A numeric, the order of splines for the baseline risk function (3 by default) |
competing_risk |
A logical indicating if the model handles with competing risks |
formSurv_CR |
In case of competing risk A formula which indicates the variables used in the survival submodel for the second event |
hazard_baseline_CR |
In case of competing risk : a character for the baseline hazard function of the second event |
sharedtype_CR |
In case of competing risk ; a character for the dependence structure |
left_trunc |
A logical indicating if the model handles with left truncated data |
Time.0 |
In case of left truncation : a vector of entry times |
S1 |
An integer : the number of QMC draws for the first step |
S2 |
An integer : the number of QMC draws for the second step |
nproc |
An integer : the number of processors for parallel computing |
clustertype |
one of the supported types from |
maxiter |
optional maximum number of iterations for the marqLevAlg iterative algorithm. |
print.info |
logical indicating if the outputs of each iteration should be written |
file |
optional character giving the name of the file where the outputs of each iteration should be written (if print.info=TRUE) |
epsa |
optional threshold for the convergence criterion based on the parameter stability. |
epsb |
optional threshold for the convergence criterion based on the objective function stability. |
epsd |
optional threshold for the relative distance to maximum. This criterion has the nice interpretation of estimating the ratio of the approximation error over the statistical error, thus it can be used for stopping the iterative process whatever the problem. |
binit |
optional initials parameters. |
Comp.Rcpp |
boolean to indicate if the computation is performed with RCPP program or R program. True by default. |
A. LONGITUDINAL SUBMODEL
The longitudinal submodel is defined by a linear mixed effects model with the residual variance which could be supposed to be time-dependent and subject-specific :
\quad\left\{\begin{array}{ll}
Y_{ij} = Y_{i}(t_{ij}) = \widetilde{Y}_i(t_{ij}) + \epsilon_{ij} = X_{ij}^{\top} \beta+Z_{ij}^{\top} b_{i}+\epsilon_{ij}, \\
\epsilon_{ij}(t_{ij}) \sim \mathcal{N}(0,\sigma_i^2(t_{ij})) \hspace{3mm} \text{with} \hspace{3mm} \log(\sigma_i(t_{ij})) = O_{ij}^{\top} \mu+M_{ij}^{\top} \tau_{i}
\end{array}
\right.
X_{ij}
, O_{ij}
, Z_{ij}
and M_{ij}
four vectors of explanatory variables for subject i
at visit j
,
respectively associated with the fixed-effect vectors \beta
and \mu
, and the subject-specific random-effect vector b_i
and \tau_i
, such as
\quad\left(\begin{array}{c}
b_{i} \\
\tau_i
\end{array}\right) \sim N\left(\left(\begin{array}{c}
0 \\
0
\end{array}\right),\left(\begin{array}{cc}
\Sigma_{b} & \Sigma_{\tau b} \\
\Sigma_{\tau b}' & \Sigma_{\tau}
\end{array}\right)\right)
Y_{i}(t_{ij}) = \tilde{Y}_i(t_{ij}) + \epsilon_{ij} = X_{ij}^{\top} \beta+Z_{ij}^{\top} b_{i}+\epsilon_{ij}
with X_{ij}
and Z_{ij}
two covariate vectors for subject i at visit j,
respectively associated with the vector of fixed effects \beta
and the vector of
subject-specific individual random effects b_i
.
The vector b_i
is assumed to be normally distributed and a specific-subject random effect on the
variance of the measure error can be added: \epsilon_{ij} \sim \mathcal{N}(0,\sigma_i^2)
and
\quad\left(\begin{array}{c}
b_{i} \\
\log \sigma_{i}
\end{array}\right) \sim \mathcal{N}\left(\left(\begin{array}{c}
0 \\
\mu_{\sigma}
\end{array}\right),\left(\begin{array}{cc}
\Sigma_{b} & 0 \\
0 & \tau_{\sigma}^{2}
\end{array}\right)\right)
B. SURVIVAL SUBMODEL
The risk function for the event $k = {1,2}$ is defined by:
\lambda_{ik}(t)=\lambda_{0k}(t) \exp \left(W_{i}^{\top} \gamma_{k}+\alpha_{1k}\tilde{y}_i(t)+\\
\alpha_{2k}\tilde{y}'_i(t)+ \alpha_{\sigma k} \sigma_i(t) \right)
with \lambda_{0k}(t)
the baseline risk function, W_{i}
a vector of baseline covariates associated with the regression coefficient \gamma_k
,
and \alpha_{1k}
, \alpha_{2k}
and \alpha_{\sigma k}
the regression coefficients associated with the current value \tilde{y}_i(t)
,
the current slope \tilde{y}'_i(t)
and the current variability \sigma_i(t)
of the marker, respectively.
Different parametric forms for the baseline risk function can be considered, such as exponential, Weibull, or, for more flexibility, a B-splines base.
A FlexVarJoint object which contains the following elements :
result
A marqLevAlg object with the results of the estimation.
table.res
The table of results : Estimation and SE
time.compute
Computation time
control
A list of control elements
#fit a joint model with competing risks and subject-specific variability
example <- lsjm(formFixed = y~visit,
formRandom = ~ visit,
formGroup = ~ID,
formSurv = Surv(time, event ==1 ) ~ 1,
timeVar = "visit",
data.long = Data_toy,
variability_hetero = TRUE,
formFixedVar =~visit,
formRandomVar =~visit,
correlated_re = TRUE,
sharedtype = c("current value", "variability"),
hazard_baseline = "Weibull",
formSlopeFixed =~1,
formSlopeRandom = ~1,
indices_beta_slope = c(2),
competing_risk = TRUE,
formSurv_CR = Surv(time, event ==2 ) ~ 1,
hazard_baseline_CR = "Weibull",
sharedtype_CR = c("current value", "variability"),
S1 = 100,
S2 = 1000,
nproc = 1,
maxiter = 100,
Comp.Rcpp = TRUE
)
summary(example)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.