Description Usage Arguments Details Value Note References See Also Examples
Shared Frailty model
Fit a shared gamma or lognormal frailty model using a semiparametric Penalized Likelihood estimation or parametric estimation on the hazard function. Lefttruncated, rightcensored data, intervalcensored data and strata (up to 6 levels) are allowed. It allows to obtain a nonparametric 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 ω_{i}, of a shared gamma frailty model for the j^{th} subject in the i^{th} group:
where λ_{0}(t) is the baseline hazard function, β 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 lognormal frailty model, we have for the j^{th} subject in the i^{th} group:
From now on, you can also consider timevarying effects covariates in your
model, see timedep
function for more details.
Joint Frailty model
Fit a joint either with gamma or lognormal frailty model for recurrent and terminal events using a penalized likelihood estimation on the hazard function or a parametric estimation. Rightcensored data and strata (up to 6 levels) for the recurrent event part are allowed. Lefttruncated 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 ω_{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 (ω_{i} for the recurrent rate and ω_{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 λ_{i} is :
where r_{0}(t) (resp. λ_{0}(t)) is the recurrent (resp. terminal) event baseline hazard function, β_{1} (resp. β_{2}) the regression coefficient vector, Z_{i}(t) the covariate vector. The random effects of frailties ω_{i} ~ Γ(1/θ,1/θ) and are iid.
The joint lognormal 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 semicompeting events. The followup 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 α parameter as there is not enough information to estimate it and thus there might be convergence problems.
In case of a lognormal 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 timevarying 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 casecontrol 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} ~ Γ(1/θ,1/θ) and the iid random effects v_{i} ~ Γ(1/η,1/η) are independent from each other. The joint model is fitted using a penalized likelihood estimation on the hazard. Rightcensored data and timevarying 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. Lefttruncated and rightcensored 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 λ_{0}(t) is the baseline hazard function, X_{ijk} denotes the covariate vector and β 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. Rightcensored data are allowed but lefttruncated 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 λ_{fi} is :
where r_{0}(resp. λ_{0}) is the recurrent (resp. terminal) event baseline hazard function, β (resp. γ) the regression coefficient vector, X_{fij}(t) the covariates vector. The random effects are ω_{f} ~ Γ(1/η,1/η) and u_{fi} ~ Γ(1/θ,1/θ).
1 2 3 4 5 6  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 = 1e3, LIMlogl = 1e3, LIMderiv = 1e3, print.times =
TRUE)

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 intervalcensored 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 AndersenGill 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 timedependent 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 intervalcensored 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 θ) and an association amongst the recurrent events (explained by the variance η). 
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 Msplines of order 4. When the user set a number of knots equals to k (n.knots=k) then the number of interior knots is (k2) and the number of splines is (k2)+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 
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 "Splinesper" using percentile with
the penalized likelihood estimation, "Piecewiseper" for piecewise constant
hazard function using percentile (not available for intervalcensored data),
"Piecewiseequi" for piecewise constant hazard function using equidistant
intervals, "Weibull" for parametric Weibull functions. Default is "Splines".
In case of 
nb.int 
Number of time intervals (between 1 and 20) for the parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 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 lognormal distribution. Default is "Gamma". Not
implemented for nested model. If 
nb.gh 
Number of nodes for the GaussianHermite 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 GaussianLaguerre 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 Bsplines. Default is 1. See 'timedep' function for more details. Not implemented for nested and joint nested frailty models. 
betaorder 
Order of the Bsplines. Default is cubic Bsplines (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 
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 ξ. 
Ksi 
Only for joint nested frailty model : input 
init.Eta 
Only for general joint and joint nested frailty models : initial value for the variance η of the frailty v_i (general joint model) and of the frailty ω_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 loglikelihood (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. 
Typical usages are for a Cox model
1  frailtyPenal(Surv(time,event)~var1+var2, data, \dots)

for a shared model
1 2  frailtyPenal(Surv(time,event)~cluster(group)+var1+var2, data,
\dots)

for a joint model
1 2 
for a joint model for clustered data
1 2 3 
for a joint model for data from nested casecontrol studies
1 2 3 
for a nested model
1 2  frailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+
var1+var2, data, \dots)

for a joint nested frailty model
1 2  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 NewtonRaphson algorithm and a steepest descent algorithm. The iterations are stopped when the difference between two consecutive loglikelihoods 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 loglikelihood is used (see Rondeau, 2003 for further details). Cubic Msplines of order 4 are used for the hazard function, and Isplines (integrated Msplines) 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 Δ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,∞[
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 θ 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 θ and the coefficient α 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 η and
ξ 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 θ of
the frailty terms and α. The variances of the second frailty term
(η) and the second coefficient ξ 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 casecontrol studies using the approach of weighted penalized maximum likelihood. For this model, only splines can be used for baseline hazards and no timevarying effects of covariates can be included. To accommodate the nested casecontrol 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).
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 loglikelihood in the semiparametric case. 
loglik 
the marginal loglikelihood 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 ("Piecewiseper", "Piecewiseequi"). 
npar 
number of parameters. 
nvar 
number of explanatory variables. 
LCV 
the approximated likelihood crossvalidation criterion in the semiparametric case (with H minus the converged Hessian matrix, and l(.) the full loglikelihood). 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 
pvalues 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 pieceiwseconstaant functions) : 1
for equidistant intervals ; 0 for intervals using percentile (note:

intcens 
Logical value. Indicator if a joint frailty model with intervalcensored data was fitted) 
theta 
variance of the gamma frailty parameter (\bold{Var}(ω_i)) 
sigma2 
variance of the lognormal frailty parameter (\bold{Var}(η_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 lognormal frailty distribution. 
BetaTpsMat 
matrix of time varyingeffects and confidence bands (the first column used for abscissa of times) 
theta_p.value 
pvalue of the Wald test for the estimated variance of the gamma frailty. 
sigma2_p.value 
pvalue of the Wald test for the estimated variance of the lognormal frailty. 
The following components are specific to joint models.
intcens 
Logical value. Indicator if a joint frailty model with intervalcensored data was fitted) 
theta 
variance of the gamma frailty parameter (\bold{Var}(ω_i)) or (\bold{Var}(u_i)) 
sigma2 
variance of the lognormal frailty parameter (\bold{Var}(η_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 α parameter was fitted 
alpha 
the coefficient α associated with the frailty parameter in the terminal hazard function. 
nbintervR 
Number of intervals (between 1 and 20) for the recurrent parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 
nbintervDC 
Number of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewiseper", "Piecewiseequi"). 
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 lognormal frailty distribution 
lineardeath.pred 
linear predictor for the terminal part : "Beta'X + alpha.log w_i" for gamma, "Beta'X + alpha.eta_i" for lognormal 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 varyingeffects and confidence bands for recurrent event (the first column used for abscissa of times of recurrence) 
BetaTpsMatDc 
matrix of time varyingeffects and confidence bands for terminal event (the first column used for abscissa of times of death) 
alpha_p.value 
pvalue of the Wald test for the estimated α. 
ncc 
Logical value whether nested casecontrol 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 
pvalue of the Wald test for the estimated variance of the cluster effect. 
eta_p.value 
pvalue 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}(ω_f)) 
alpha 
the power coefficient α associated with the frailty parameter (u_{fi}) in the terminal event hazard function. 
ksi 
the power coefficient ξ associated with the frailty parameter (ω_f) in the recurrent event hazard function. 
indic_alpha 
indicator if a joint frailty model with α parameter was fitted or not. 
indic_ksi 
indicator if a joint frailty model with ξ parameter was fitted or not. 
frailty.fam.pred 
empirical Bayes prediction of the frailty term by family. 
eta_p.value 
pvalue of the Wald test for the estimated variance of the cluster effect. 
alpha_p.value 
pvalue of the Wald test for the estimated power coefficient α. 
ksi_p.value 
pvalue of the Wald test for the estimated power coefficient ξ. 
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.
I. Jazic, S. Haneuse, B. French, G. MacGrogan, and V. Rondeau. Design and analysis of nested casecontrol 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), 152.
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, 128.
Y. Mazroui, S. MathoulinPelissier, 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, 1112, 11621176.
V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the dependance between clustered times to tumour progression and deaths: A metaanalysis of chemotherapy in head and neck cancer. Statistical methods in medical research 897, 119.
V. Rondeau, S. MathoulinPellissier, H. JacqminGadda, 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, 708721.
V. Rondeau, L. Filleul, P. Joly (2006). Nested frailty models using maximum penalized likelihood estimation. Statistics in Medicine, 25, 40364052.
V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gammafrailty model. Lifetime Data Analysis 9, 139153.
C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in survival analysis. Biometrics 47, 461466.
D. Marquardt (1963). An algorithm for leastsquares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431441.
SurvIC
, cluster
,
subcluster
, terminal
, num.id
,
timedep
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148  ## Not run:
### 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)
# lognormal 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)
# Gaptime
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")
# lognormal 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 semicompeting 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))
### 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)
JoiNes_GapSpline < 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.