frailtyPenal {frailtypack}

Fit a Shared, Joint or Nested Frailty model

Description

Shared Frailty model

Fit 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 model

Fit 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 model

Data should be ordered according to cluster and subcluster

Fit 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 Model

Fit 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)

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 VALUES

The 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 DESIGN

It 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).

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.

Note

From a prediction aim, we recommend you to input a data sorted by the group variable with numerical numbers from 1 to n (number of groups). In case of a nested model, we recommend you to input a data sorted by the group variable then sorted by the subgroup variable both with numerical numbers from 1 to n (number of groups) and from 1 to m (number or subgroups). "kappa" and "n.knots" are the arguments that the user have to change if the fitted model does not converge. "n.knots" takes integer values between 4 and 20. But with n.knots=20, the model would take a long time to converge. So, usually, begin first with n.knots=7, and increase it step by step until it converges. "kappa" only takes positive values. So, choose a value for kappa (for instance 10000), and if it does not converge, multiply or divide this value by 10 or 5 until it converges.

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.

See Also

SurvIC, cluster, subcluster, terminal, num.id, timedep

Examples

###---  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))


## Not run: 
###--- 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")

## End(Not run)