frailtyPenal: Fit a Shared, Joint or Nested Frailty model

Description

Shared Frailty modelFit a shared gamma or log-normal frailty model using a semiparametric Penalized Likelihood estimation or parametric estimation on the hazard function. Left-truncated, right-censored data, interval-censored data and strata (up to 6 levels) are allowed. It allows to obtain a non-parametric smooth hazard of survival function. This approach is different from the partial penalized likelihood approach of Therneau et al.The hazard function, conditional on the frailty term $\omega$_i, of a shared gamma frailty model for the j^th subject in the i^th group: where $\lambda$₀(t) is the baseline hazard function, $\beta$ the vector of the regression coefficient associated to the covariate vector $Z$_ij for the j^th individual in the i^th group.Otherwise, in case of a shared log-normal frailty model, we have for the j^th subject in the i^th group: From now on, you can also consider time-varying effects covariates in your model, see timedep function for more details.Joint Frailty modelFit a joint either with gamma or log-normal frailty model for recurrent and terminal events using a penalized likelihood estimation on the hazard function or a parametric estimation. Right-censored data and strata (up to 6 levels) for the recurrent event part are allowed. Left-truncated data is not possible. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring.There is two kinds of joint frailty models that can be fitted with frailtyPenal :- The first one (Rondeau et al. 2007) includes a common frailty term to the individuals $\omega$_i for the two rates which will take into account the heterogeneity in the data, associated with unobserved covariates. The frailty term acts differently for the two rates ($\omega$_i for the recurrent rate and $\omega$_i^\eqn{\alpha} for the death rate). The covariates could be different for the recurrent rate and death rate.For the j^th recurrence (j=1,...,n_i) and the i^th subject (i=1,...,G), the joint gamma frailty model for recurrent event hazard function $r$_ij(.) and death rate $\lambda$_i is :where $r$₀(t) (resp. $\lambda$₀(t)) is the recurrent (resp. terminal) event baseline hazard function, $\beta$₁ (resp. $\beta$₂) the regression coefficient vector, $Z$_i(t) the covariate vector. The random effects of frailties $\omega$_i ~ $\Gamma$(1/$\theta$,1/$\theta$) and are iid.The joint log-normal frailty model will be : - The second one (Rondeau et al. 2011) is quite similar but the frailty term is common to the individuals from a same group. This model is useful for the joint modelling two clustered survival outcomes. This joint models have been developed for clustered semi-competing events. The follow-up of each of the two competing outcomes stops when the event occurs. In this case, j is for the subject and i for the cluster.It should be noted that in these models it is not recommended to include $\alpha$ parameter as there is not enough information to estimate it and thus there might be convergence problems.In case of a log-normal distribution of the frailties, we will have : This joint frailty model can also be applied to clustered recurrent events and a terminal event (example on "readmission" data below).From now on, you can also consider time-varying effects covariates in your model, see timedep function for more details.There is a possibility to use a weighted penalized maximum likelihood approach for nested case-control design, in which risk set sampling is performed based on a single outcome (Jazic et al., Submitted).General Joint Frailty model Fit a general joint frailty model for recurrent and terminal events considering two independent frailty terms. The frailty term $u$_i represents the unobserved association between recurrences and death. The frailty term $v$_i is specific to the recurrent event rate. Thus, the general joint frailty model is:where the $iid$ random effects $u$_i ~ $\Gamma$(1/$\theta$,1/$\theta$) and the $iid$ random effects $v$_i ~ $\Gamma$(1/$\eta$,1/$\eta$) are independent from each other. The joint model is fitted using a penalized likelihood estimation on the hazard. Right-censored data and time-varying covariates $Z$_i(t) are allowed.Nested Frailty modelData should be ordered according to cluster and subclusterFit a nested frailty model using a Penalized Likelihood on the hazard function or using a parametric estimation. Nested frailty models allow survival studies for hierarchically clustered data by including two iid gamma random effects. Left-truncated and right-censored data are allowed. Stratification analysis is allowed (maximum of strata = 2). The hazard function conditional on the two frailties $v$_i and $w$_ij for the k^th individual of the j^th subgroup of the i^th group is :where $\lambda$₀(t) is the baseline hazard function, $X$_ijk denotes the covariate vector and $\beta$ the corresponding vector of regression parameters.Joint Nested Frailty ModelFit a joint model for recurrent and terminal events using a penalized likelihood on the hazard functions or a parametric estimation. Right-censored data are allowed but left-truncated data and stratified analysis are not allowed.Joint nested frailty models allow studying, jointly, survival processes of recurrent and terminal events for hierarchically clustered data, by considering the terminal event as an informative censoring and by including two iid gamma random effects.The joint nested frailty model includes two shared frailty terms, one for the subgroup ($u$_fi) and one for the group ($w$_f) into the hazard functions. This random effects account the heterogeneity in the data, associated with unobserved covariates. The frailty terms act differently for the two rates ($u$_fi, $w$_f^\eqn{\xi} for the recurrent rate and $u$_fi^\eqn{\alpha}, $w$_i for the terminal event rate). The covariates could be different for the recurrent rate and death rate.For the j^th recurrence (j = 1, ..., $n$_i) of the i^th individual (i = 1, ..., $m$_f) of the $f$^th group (f = 1,..., n), the joint nested gamma frailty model for recurrent event hazard function $r$_fij(.) and for terminal event hazard function $\lambda$_fi is :where $r$₀(resp. $\lambda$₀) is the recurrent (resp. terminal) event baseline hazard function, $\beta$ (resp. $\gamma$) the regression coefficient vector, $X$_fij(t) the covariates vector. The random effects are $\omega$_f ~ $\Gamma$(1/$\eta$,1/$\eta$) and $u$_fi ~ $\Gamma$(1/$\theta$,1/$\theta$).

Usage

frailtyPenal(formula, formula.terminalEvent, data, recurrentAG = FALSE,
cross.validation = FALSE, jointGeneral,n.knots, kappa, maxit = 300, hazard =
"Splines", nb.int, RandDist = "Gamma", nb.gh, nb.gl, betaknots = 1, betaorder = 3,
initialize = TRUE, init.B, init.Theta, init.Alpha, Alpha, init.Ksi, Ksi,
init.Eta, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times =
TRUE)

Value

The following components are included in a 'frailtyPenal' object for each model.

b: sequence of the corresponding estimation of the coefficients for the hazard functions (parametric or semiparametric), the random effects variances and the regression coefficients.
call: The code used for the model.
formula: the formula part of the code used for the model.
coef: the regression coefficients.
cross.Val: Logical value. Is cross validation procedure used for estimating the smoothing parameters in the penalized likelihood estimation?
DoF: Degrees of freedom associated with the "kappa".
groups: the maximum number of groups used in the fit.
kappa: A vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components.
loglikPenal: the complete marginal penalized log-likelihood in the semiparametric case.
loglik: the marginal log-likelihood in the parametric case.
n: the number of observations used in the fit.
n.events: the number of events observed in the fit.
n.iter: number of iterations needed to converge.
n.knots: number of knots for estimating the baseline functions in the penalized likelihood estimation.
n.strat: number of stratum.
varH: the variance matrix of all parameters before positivity constraint transformation. Then, the delta method is needed to obtain the estimated variance parameters. That is why some variances don't match with the printed values at the end of the model.
varHIH: the robust estimation of the variance matrix of all parameters.
x: matrix of times where both survival and hazard function are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times.
lam: array (dim=3) of hazard estimates and confidence bands.
surv: array (dim=3) of baseline survival estimates and confidence bands.
median: The value of the median survival and its confidence bands. If there are two stratas or more, the first value corresponds to the value for the first strata, etc.
nbintervR: Number of intervals (between 1 and 20) for the parametric hazard functions ("Piecewise-per", "Piecewise-equi").
npar: number of parameters.
nvar: number of explanatory variables.
LCV: the approximated likelihood cross-validation criterion in the semiparametric case (with H minus the converged Hessian matrix, and l(.) the full log-likelihood).$$LCV=\frac{1}{n}(trace(H^{-1}_{pl}H) - l(.))$$
AIC: the Akaike information Criterion for the parametric case.$$AIC=\frac{1}{n}(np - l(.))$$
n.knots.temp: initial value for the number of knots.
shape.weib: shape parameter for the Weibull hazard function.
scale.weib: scale parameter for the Weibull hazard function.
martingale.res: martingale residuals for each cluster.
martingaleCox: martingale residuals for observation in the Cox model.
Frailty: Logical value. Was model with frailties fitted ?
frailty.pred: empirical Bayes prediction of the frailty term (ie, using conditional posterior distributions).
frailty.var: variance of the empirical Bayes prediction of the frailty term (only for gamma frailty models).
frailty.sd: standard error of the frailty empirical Bayes prediction (only for gamma frailty models).
global_chisq: a vector with the values of each multivariate Wald test.
dof_chisq: a vector with the degree of freedom for each multivariate Wald test.
global_chisq.test: a binary variable equals to 0 when no multivariate Wald is given, 1 otherwise.
p.global_chisq: a vector with the p_values for each global multivariate Wald test.
names.factor: Names of the "as.factor" variables.
Xlevels: vector of the values that factor might have taken.
contrasts: type of contrast for factor variable.
beta_p.value: p-values of the Wald test for the estimated regression coefficients.

The following components are specific to shared models.

equidistant: Indicator for the intervals used the estimation of baseline hazard functions (for splines or pieceiwse-constaant functions) : 1 for equidistant intervals ; 0 for intervals using percentile (note: equidistant = 2 in case of parametric estimation using Weibull distribution).
intcens: Logical value. Indicator if a joint frailty model with interval-censored data was fitted)
theta: variance of the gamma frailty parameter $(\bold{Var}(\omega_i))$
sigma2: variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$
linear.pred: linear predictor: uses simply "Beta'X" in the cox proportional hazard model or "Beta'X + log w_i" in the shared gamma frailty models, otherwise uses "Beta'X + w_i" for log-normal frailty distribution.
BetaTpsMat: matrix of time varying-effects and confidence bands (the first column used for abscissa of times)
theta_p.value: p-value of the Wald test for the estimated variance of the gamma frailty.
sigma2_p.value: p-value of the Wald test for the estimated variance of the log-normal frailty.

The following components are specific to joint models.

intcens: Logical value. Indicator if a joint frailty model with interval-censored data was fitted)
theta: variance of the gamma frailty parameter $(\bold{Var}(\omega_i))$ or $(\bold{Var}(u_i))$
sigma2: variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$ or $(\bold{Var}(v_i))$
eta: variance of the second gamma frailty parameter in general joint frailty models $(\bold{Var}(v_i))$
indic_alpha: indicator if a joint frailty model with $\alpha$ parameter was fitted
alpha: the coefficient $\alpha$ associated with the frailty parameter in the terminal hazard function.
nbintervR: Number of intervals (between 1 and 20) for the recurrent parametric hazard functions ("Piecewise-per", "Piecewise-equi").
nbintervDC: Number of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewise-per", "Piecewise-equi").
nvar: A vector with the number of covariates of each type of hazard function as components.
nvarRec: number of recurrent explanatory variables.
nvarEnd: number of death explanatory variables.
noVar1: indicator of recurrent explanatory variables.
noVar2: indicator of death explanatory variables.
xR: matrix of times where both survival and hazard function are estimated for the recurrent event. By default seq(0,max(time),length=99), where time is the vector of survival times.
xD: matrix of times for the terminal event.
lamR: array (dim=3) of hazard estimates and confidence bands for recurrent event.
lamD: the same value as lamR for the terminal event.
survR: array (dim=3) of baseline survival estimates and confidence bands for recurrent event.
survD: the same value as survR for the terminal event.
martingale.res: martingale residuals for each cluster (recurrent).
martingaledeath.res: martingale residuals for each cluster (death).
linear.pred: linear predictor: uses "Beta'X + log w_i" in the gamma frailty model, otherwise uses "Beta'X + eta_i" for log-normal frailty distribution
lineardeath.pred: linear predictor for the terminal part : "Beta'X + alpha.log w_i" for gamma, "Beta'X + alpha.eta_i" for log-normal frailty distribution
Xlevels: vector of the values that factor might have taken for the recurrent part.
contrasts: type of contrast for factor variable for the recurrent part.
Xlevels2: vector of the values that factor might have taken for the death part.
contrasts2: type of contrast for factor variable for the death part.
BetaTpsMat: matrix of time varying-effects and confidence bands for recurrent event (the first column used for abscissa of times of recurrence)
BetaTpsMatDc: matrix of time varying-effects and confidence bands for terminal event (the first column used for abscissa of times of death)
alpha_p.value: p-value of the Wald test for the estimated $\alpha$.
ncc: Logical value whether nested case-control design with weights was used for the joint model.

The following components are specific to nested models.

alpha: variance of the cluster effect $(\bold{Var}(v_{i}))$
eta: variance of the subcluster effect $(\bold{Var}(w_{ij}))$
subgroups: the maximum number of subgroups used in the fit.
frailty.pred.group: empirical Bayes prediction of the frailty term by group.
frailty.pred.subgroup: empirical Bayes prediction of the frailty term by subgroup.
linear.pred: linear predictor: uses "Beta'X + log v_i.w_ij".
subgbyg: subgroup by group.
n.strat: A vector with the number of covariates of each type of hazard function as components.
alpha_p.value: p-value of the Wald test for the estimated variance of the cluster effect.
eta_p.value: p-value of the Wald test for the estimated variance of the subcluster effect.

The following components are specific to joint nested frailty models.

theta: variance of the subcluster effect $(\bold{Var}(u_{fi}))$
eta: variance of the cluster effect $(\bold{Var}(\omega_f))$
alpha: the power coefficient $\alpha$ associated with the frailty parameter ($u_{fi}$) in the terminal event hazard function.
ksi: the power coefficient $\xi$ associated with the frailty parameter ($\omega_f$) in the recurrent event hazard function.
indic_alpha: indicator if a joint frailty model with $\alpha$ parameter was fitted or not.
indic_ksi: indicator if a joint frailty model with $\xi$ parameter was fitted or not.
frailty.fam.pred: empirical Bayes prediction of the frailty term by family.
eta_p.value: p-value of the Wald test for the estimated variance of the cluster effect.
alpha_p.value: p-value of the Wald test for the estimated power coefficient $\alpha$.
ksi_p.value: p-value of the Wald test for the estimated power coefficient $\xi$.

Arguments

formula: a formula object, with the response on the left of a $\sim$ operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function like in survival package. In case of interval-censored data, the response must be an object as returned by the 'SurvIC' function from this package. Interactions are possible using * or :.
formula.terminalEvent: only for joint and joint nested frailty models : a formula object, only requires terms on the right to indicate which variables are modelling the terminal event. Interactions are possible using * or :.
data: a 'data.frame' with the variables used in 'formula'.
recurrentAG: Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent covariates. The default is FALSE.
cross.validation: Logical value. Is cross validation procedure used for estimating smoothing parameter in the penalized likelihood estimation? If so a search of the smoothing parameter using cross validation is done, with kappa as the seed. The cross validation is not implemented for several strata, neither for interval-censored data. The cross validation has been implemented for a Cox proportional hazard model, with no covariates. The default is FALSE.
jointGeneral: Logical value. Does the model include two independent random effects? If so, this will fit a general joint frailty model with an association between the recurrent events and a terminal event (explained by the variance $\theta$) and an association amongst the recurrent events (explained by the variance $\eta$).
n.knots: integer giving the number of knots to use. Value required in the penalized likelihood estimation. It corresponds to the (n.knots+2) splines functions for the approximation of the hazard or the survival functions. We estimate I or M-splines of order 4. When the user set a number of knots equals to k (n.knots=k) then the number of interior knots is (k-2) and the number of splines is (k-2)+order. Number of knots must be between 4 and 20. (See Note)
kappa: positive smoothing parameter in the penalized likelihood estimation. In a stratified shared model, this argument must be a vector with kappas for both strata. In a stratified joint model, this argument must be a vector with kappas for both strata for recurrent events plus one kappa for terminal event. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit (penalized log likelihood). To obtain an initial value for kappa, a solution is to fit the corresponding shared frailty model using cross validation (See cross.validation). We advise the user to identify several possible tuning parameters, note their defaults and look at the sensitivity of the results to varying them. Value required. (See Note).
maxit: maximum number of iterations for the Marquardt algorithm. Default is 300
hazard: Type of hazard functions: "Splines" for semiparametric hazard functions using equidistant intervals or "Splines-per" using percentile with the penalized likelihood estimation, "Piecewise-per" for piecewise constant hazard function using percentile (not available for interval-censored data), "Piecewise-equi" for piecewise constant hazard function using equidistant intervals, "Weibull" for parametric Weibull functions. Default is "Splines". In case of jointGeneral = TRUE or if a joint nested frailty model is fitted, only hazard = "Splines" can be chosen.
nb.int: Number of time intervals (between 1 and 20) for the parametric hazard functions ("Piecewise-per", "Piecewise-equi"). In a joint model, you need to specify a number of time interval for both recurrent hazard function and the death hazard function (vector of length 2).
RandDist: Type of random effect distribution: "Gamma" for a gamma distribution, "LogN" for a log-normal distribution. Default is "Gamma". Not implemented for nested model. If jointGeneral = TRUE or if a joint nested frailty model is fitted, the log-normal distribution cannot be chosen.
nb.gh: Number of nodes for the Gaussian-Hermite quadrature. It can be chosen among 5, 7, 9, 12, 15, 20 and 32. The default is 20 if hazard = "Splines", 32 otherwise.
nb.gl: Number of nodes for the Gaussian-Laguerre quadrature. It can be chosen between 20 and 32. The default is 20 if hazard = "Splines", 32 otherwise.
betaknots: Number of inner knots used for the estimation of B-splines. Default is 1. See 'timedep' function for more details. Not implemented for nested and joint nested frailty models.
betaorder: Order of the B-splines. Default is cubic B-splines (order = 3). See 'timedep' function for more details. Not implemented for nested and joint nested frailty models.
initialize: Logical value, only for joint nested frailty models. Option TRUE indicates fitting an appropriate standard joint frailty model (without group effect, only the subgroup effect) to provide initial values for the joint nested model. Default is TRUE.
init.B: A vector of initial values for regression coefficients. This vector should be of the same size as the whole vector of covariates with the first elements for the covariates related to the recurrent events and then to the terminal event (interactions in the end of each component). Default is 0.1 for each (for Cox and shared model) or 0.5 (for joint and joint nested frailty models).
init.Theta: Initial value for variance of the frailties.
init.Alpha: Only for joint and joint nested frailty models : initial value for parameter alpha.
Alpha: Only for joint and joint nested frailty model : input "None" so as to fit a joint model without the parameter alpha.
init.Ksi: Only for joint nested frailty model : initial value for parameter $\xi$.
Ksi: Only for joint nested frailty model : input "None" indicates a joint nested frailty model without the parameter $\xi$.
init.Eta: Only for general joint and joint nested frailty models : initial value for the variance $\eta$ of the frailty $v_i$ (general joint model) and of the frailty $\omega_i$ (joint nested frailty model).
LIMparam: Convergence threshold of the Marquardt algorithm for the parameters (see Details), $10^{-3}$ by default.
LIMlogl: Convergence threshold of the Marquardt algorithm for the log-likelihood (see Details), $10^{-3}$ by default.
LIMderiv: Convergence threshold of the Marquardt algorithm for the gradient (see Details), $10^{-3}$ by default.
print.times: a logical parameter to print iteration process. Default is TRUE.

Details

Typical usages are for a Cox model

frailtyPenal(Surv(time,event)~var1+var2, data, \dots)

for a shared model

frailtyPenal(Surv(time,event)~cluster(group)+var1+var2, data,
\dots)

for a joint model

frailtyPenal(Surv(time,event)~cluster(group)+var1+var2+
var3+terminal(death), formula.terminalEvent=~ var1+var4, data, \dots)

for a joint model for clustered data

frailtyPenal(Surv(time,event)~cluster(group)+num.id(group2)+
var1+var2+var3+terminal(death), formula.terminalEvent=~var1+var4, data,
\dots)

for a joint model for data from nested case-control studies

frailtyPenal(Surv(time,event)~cluster(group)+num.id(group2)+
var1+var2+var3+terminal(death)+wts(wts.ncc),
formula.terminalEvent=~var1+var4, data, \dots)

for a nested model

frailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+
var1+var2, data, \dots)

for a joint nested frailty model

frailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+
var1+var2++terminal(death), formula.terminalEvent=~var1+var4, data, \dots)

The estimated parameter are obtained using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm. The iterations are stopped when the difference between two consecutive log-likelihoods was small $(<10^{-3})$, the estimated coefficients were stable (consecutive values $(<10^{-3})$, and the gradient small enough $(<10^{-3})$. When frailty parameter is small, numerical problems may arise. To solve this problem, an alternative formula of the penalized log-likelihood is used (see Rondeau, 2003 for further details). Cubic M-splines of order 4 are used for the hazard function, and I-splines (integrated M-splines) are used for the cumulative hazard function.The inverse of the Hessian matrix is the variance estimator and to deal with the positivity constraint of the variance component and the spline coefficients, a squared transformation is used and the standard errors are computed by the $\Delta$-method (Knight & Xekalaki, 2000). The smooth parameter can be chosen by maximizing a likelihood cross validation criterion (Joly and other, 1998). The integrations in the full log likelihood were evaluated using Gaussian quadrature. Laguerre polynomials with 20 points were used to treat the integrations on $[0,\infty[$INITIAL VALUESThe splines and the regression coefficients are initialized to 0.1. In case of shared model, the program fits, firstly, an adjusted Cox model to give new initial values for the splines and the regression coefficients. The variance of the frailty term $\theta$ is initialized to 0.1. Then, a shared frailty model is fitted.In case of a joint frailty model, the splines and the regression coefficients are initialized to 0.5. The program fits an adjusted Cox model to have new initial values for the regression and the splines coefficients. The variance of the frailty term $\theta$ and the coefficient $\alpha$ associated in the death hazard function are initialized to 1. Then, it fits a joint frailty model.In case of a general joint frailty model we need to initialize the jointGeneral logical value to TRUE.In case of a nested model, the program fits an adjusted Cox model to provide new initial values for the regression and the splines coefficients. The variances of the frailties are initialized to 0.1. Then, a shared frailty model with covariates with only subgroup frailty is fitted to give a new initial value for the variance of the subgroup frailty term. Then, a shared frailty model with covariates and only group frailty terms is fitted to give a new initial value for the variance of the group frailties. In a last step, a nested frailty model is fitted.In case of a joint nested model, the splines and the regression coefficients are initialized to 0.5 and the variances of the frailty terms $\eta$ and $\xi$ are initialized to 1. If the option 'initialize' is TRUE, the program fits a joint frailty model to provide initial values for the splines, covariates coefficients, variance $\theta$ of the frailty terms and $\alpha$. The variances of the second frailty term ($\eta$) and the second coefficient $\xi$ are initialized to 1. Then, a joint nested frailty model is fitted.NCC DESIGNIt is possible to fit a joint frailty model for data from nested case-control studies using the approach of weighted penalized maximum likelihood. For this model, only splines can be used for baseline hazards and no time-varying effects of covariates can be included. To accommodate the nested case-control design, the formula for the recurrent events should simply include the special term wts(wts.ncc), where wts.ncc refers to a column of prespecified weights in the data set for every observation. For details, see Jazic et al., Submitted (available on request from the package authors).

References

I. Jazic, S. Haneuse, B. French, G. MacGrogan, and V. Rondeau. Design and analysis of nested case-control studies for recurrent events subject to a terminal event. Submitted.

A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels and V. Rondeau (2017). Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event. Journal of Statistical Software 81(3), 1-52.

V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametric estimation. Journal of Statistical Software 47, 1-28.

Y. Mazroui, S. Mathoulin-Pelissier, P. Soubeyranb and V. Rondeau (2012) General joint frailty model for recurrent event data with a dependent terminalevent: Application to follicular lymphoma data. Statistics in Medecine, 31, 11-12, 1162-1176.

V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the dependance between clustered times to tumour progression and deaths: A meta-analysis of chemotherapy in head and neck cancer. Statistical methods in medical research 897, 1-19.

V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P. Soubeyran (2007). Joint frailty models for recurring events and death using maximum penalized likelihood estimation:application on cancer events. Biostatistics 8,4, 708-721.

V. Rondeau, L. Filleul, P. Joly (2006). Nested frailty models using maximum penalized likelihood estimation. Statistics in Medicine, 25, 4036-4052.

V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.

C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in survival analysis. Biometrics 47, 461-466.

D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.

Examples

Run this code



# \donttest{

###---  COX proportional hazard model (SHARED without frailties) ---###
###---  estimated with penalized likelihood ---###

data(kidney)
frailtyPenal(Surv(time,status)~sex+age,
n.knots=12,kappa=10000,data=kidney)

###---  Shared Frailty model  ---###

frailtyPenal(Surv(time,status)~cluster(id)+sex+age,
n.knots=12,kappa=10000,data=kidney)

#-- with an initialisation of regression coefficients

frailtyPenal(Surv(time,status)~cluster(id)+sex+age,
n.knots=12,kappa=10000,data=kidney,init.B=c(-1.44,0))

#-- with truncated data

data(dataNested)

frailtyPenal(Surv(t1,t2,event) ~ cluster(group),
data=dataNested,n.knots=10,kappa=10000,
cross.validation=TRUE,recurrentAG=FALSE)

#-- stratified analysis

data(readmission)
frailtyPenal(Surv(time,event)~cluster(id)+dukes+strata(sex),
n.knots=10,kappa=c(10000,10000),data=readmission)

#-- recurrentAG=TRUE

frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex+dukes+
charlson,data=readmission,n.knots=6,kappa=1e5,recurrentAG=TRUE)

#-- cross.validation=TRUE

frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex+dukes+
charlson,data=readmission,n.knots=6,kappa=5000,recurrentAG=TRUE,
cross.validation=TRUE)

#-- log-normal distribution

frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex+dukes+
charlson,data=readmission,n.knots=6,kappa=5000,recurrentAG=TRUE,
RandDist="LogN")

###--- Joint Frailty model (recurrent and terminal events) ---###

data(readmission)
#-- Gap-time
modJoint.gap <- frailtyPenal(Surv(time,event)~cluster(id)+sex+dukes+charlson+
terminal(death),formula.terminalEvent=~sex+dukes+charlson,
data=readmission,n.knots=14,kappa=c(9.55e+9,1.41e+12),
recurrentAG=FALSE)

#-- Calendar time
modJoint.calendar <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+
sex+dukes+charlson+terminal(death),formula.terminalEvent=~sex
+dukes+charlson,data=readmission,n.knots=10,kappa=c(9.55e9,1.41e12),
recurrentAG=TRUE)

#-- without alpha parameter
modJoint.gap <- frailtyPenal(Surv(time,event)~cluster(id)+sex+dukes+charlson+
terminal(death),formula.terminalEvent=~sex+dukes+charlson,
data=readmission,n.knots=10,kappa=c(9.55e9,1.41e12),
recurrentAG=FALSE,Alpha="None")

#-- log-normal distribution

modJoint.log <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex
+dukes+charlson+terminal(death),formula.terminalEvent=~sex
+dukes+charlson,data=readmission,n.knots=10,kappa=c(9.55e9,1.41e12),
recurrentAG=TRUE,RandDist="LogN")

###--- Joint frailty model for NCC data ---###
data(dataNCC)
modJoint.ncc <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+cov1
+cov2+terminal(death)+wts(ncc.wts), formula.terminalEvent=~cov1+cov2,
data=dataNCC,n.knots=8,kappa=c(1.6e+10, 5.0e+03),recurrentAG=TRUE, RandDist="LogN") 


###--- Joint Frailty model for clustered data ---###

#-- here is generated cluster (5 clusters)
readmission <- transform(readmission,group=id%%5+1)

#-- exclusion all recurrent events --#
#--  to obtain framework of semi-competing risks --#
readmission2 <- subset(readmission, (t.start == 0 & event == 1) | event == 0)

joi.clus.gap <- frailtyPenal(Surv(time,event)~cluster(group)+
num.id(id)+dukes+charlson+sex+chemo+terminal(death),
formula.terminalEvent=~dukes+charlson+sex+chemo,
data=readmission2,recurrentAG=FALSE, n.knots=8,
kappa=c(1.e+10,1.e+10) ,Alpha="None")


###--- General Joint model (recurrent and terminal events) 
###--- with 2 covariates ---###

data(readmission)
modJoint.general <- frailtyPenal(Surv(time,event) ~ cluster(id) + dukes +
charlson + sex  + chemo + terminal(death), 
formula.terminalEvent = ~ dukes + charlson + sex + chemo,
data = readmission, jointGeneral = TRUE,  n.knots = 8, 
kappa = c(2.11e+08, 9.53e+11))


###--- Nested Frailty model ---###

##***** WARNING *****##
# Data should be ordered according to cluster and subcluster

data(dataNested)
modClu <- frailtyPenal(Surv(t1,t2,event)~cluster(group)+
subcluster(subgroup)+cov1+cov2,data=dataNested,
n.knots=8,kappa=50000)

modClu.str <- frailtyPenal(Surv(t1,t2,event)~cluster(group)+
subcluster(subgroup)+cov1+strata(cov2),data=dataNested,
n.knots=8,kappa=c(50000,50000))

# }
if (FALSE) {
###--- Joint Nested Frailty model ---###

#-- here is generated cluster (30 clusters)
readmissionNested <- transform(readmission,group=id%%30+1)

modJointNested_Splines <- frailtyPenal(formula = Surv(t.start, t.stop, event) 
~ subcluster(id) + cluster(group) + dukes + terminal(death), 
formula.terminalEvent = ~dukes, data = readmissionNested, recurrentAG = TRUE, 
n.knots = 8, kappa = c(9.55e+9, 1.41e+12), initialize = TRUE)

modJointNested_Weib <- frailtyPenal(Surv(t.start,t.stop,event)~subcluster(id)
+cluster(group)+dukes+ terminal(death),formula.terminalEvent=~dukes, 
hazard = ('Weibull'), data=readmissionNested,recurrentAG=TRUE, initialize = FALSE)

JoiNesGapSpline <- frailtyPenal(formula = Surv(time, event) 
~ subcluster(id) + cluster(group) + dukes + terminal(death), 
formula.terminalEvent = ~dukes, data = readmissionNested, 
recurrentAG = FALSE, n.knots = 8, kappa = c(9.55e+9, 1.41e+12), 
initialize = TRUE, init.Alpha = 1.091, Ksi = "None")
}

Run the code above in your browser using DataLab