Package 'frailtypack'

General Frailty models: shared, joint and nested frailty models with prediction; Evaluation of Failure-Time Surrogate Endpoints

Description

Frailtypack fits several classes of frailty models using a penalized likelihood estimation on the hazard function but also a parametric estimation.

Details

The following several classes of frailty models using a penalized likelihood estimation on the hazard function but also a parametric estimation can be fit using this R package:

1) A shared frailty model (with gamma or log-normal frailty distribution) and Cox proportional hazard model. Clustered and recurrent survival times can be studied.

2) Additive frailty models for proportional hazard models with two correlated random effects (intercept random effect with random slope).

3) Nested frailty models for hierarchically clustered data (with 2 levels of clustering) by including two iid gamma random effects.

4) Joint frailty models in the context of the joint modelling for recurrent events with terminal event for clustered data or not. A joint frailty model for two semi-competing risks and clustered data is also proposed.

5) Joint general frailty models in the context of the joint modelling for recurrent events with terminal event data with two independent frailty terms.

6) Joint Nested frailty models in the context of the joint modelling for recurrent events with terminal event, for hierarchically clustered data (with two levels of clustering) by including two iid gamma random effects.

7) Multivariate joint frailty models for two types of recurrent events and a terminal event.

8) Joint models for longitudinal data and a terminal event.

9) Trivariate joint models for longitudinal data, recurrent events and a terminal event.

10) Joint frailty models for the validation of surrogate endpoints in multiple randomized clinical trials with failure-time and/or longitudinal endpoints with the possibility to use a mediation analysis model.

11) Conditional and Marginal two-part joint models for longitudinal semicontinuous data and a terminal event.

12) Joint frailty-copula models for the validation of surrogate endpoints in multiple randomized clinical trials with failure-time endpoints.

13) Generalized shared and joint frailty models for recurrent and terminal events. Proportional hazards (PH), additive hazard (AH), proportional odds (PO) and probit models are available in a fully parametric framework.

14) Competing Joint Frailty Model: A single type of recurrent event and two terminal events.

The package includes concordance measures for Cox proportional hazards models and for shared frailty models. Now, you can also consider time-varying covariates effects in Cox, shared and joint frailty models (1-5). Some of the Fortran routines in the package can speed-up computation time by making use of parallelization through OpenMP. Moreover, the package can be used with its shiny application, in a local mode or by following the link below.

Package:	frailtypack
Type:	Package
Version:	3.6.5
Date:	2024-12-09
License:	GPL (>= 2.0)
LazyLoad:	no

Author(s)

Virginie Rondeau, Juan R. Gonzalez, Yassin Mazroui, Audrey Mauguen, Amadou Diakite, Alexandre Laurent, Myriam Lopez, Agnieszka Krol, Casimir L. Sofeu, Denis Rustand, Quentin Le Coent, Lacey Etzkorn and Romain Pierlot

References

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 Virginie 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 and J. R. Gonzalez (2005). Frailtypack: A computer program for the analysis of correlated failure time data using penalized likelihood estimation. Computer Methods and Programs in Biomedicine 80, 2, 154-164.

V. Rondeau, S. Michiels, B. Liquet, and J. P. Pignon (2008). Investigating trial and treatment heterogeneity in an individual patient data meta-analysis of survival data by mean of the maximum penalized likelihood approach. Statistics in Medecine, 27, 1894-1910.

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, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.

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

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

Examples

## Not run: 

###--- Additive model with 1 covariate ---###

data(dataAdditive)
modAdd <- additivePenal(Surv(t1,t2,event)~
cluster(group)+var1+slope(var1),
correlation=TRUE,data=dataAdditive,
n.knots=8,kappa=10000,hazard="Splines")

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

data(readmission)
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(100,100),
recurrentAG=FALSE,hazard="Splines")

###--- 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 model (or hierarchical model) with 2 covariates ---###

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

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

###--- Semiparametric Shared model ---###

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

###--- Parametric Shared model ---###

data(readmission)
sha.p <- frailtyPenal(Surv(t.start,t.stop,event)~
cluster(id)+sex+dukes+charlson,
data=readmission,recurrentAG=TRUE,
hazard="Piecewise-per",nb.int=6)

###--- Joint model for longitudinal ---###
###--- data and a terminal event ---###

data(colorectal)
data(colorectalLongi)

# Survival data preparation - only terminal events 
colorectalSurv <- subset(colorectal, new.lesions == 0)

model.weib.RE <- longiPenal(Surv(time1, state) ~ age + treatment + who.PS 
+ prev.resection, tumor.size ~  year * treatment + age + who.PS ,
colorectalSurv,	data.Longi = colorectalLongi, 
random = c("1", "year"), id = "id", link = "Random-effects", 
left.censoring = -3.33, hazard = "Weibull")

###--- Trivariate joint model for longitudinal ---###
###--- data, recurrent and terminal events ---###

data(colorectal)
data(colorectalLongi)

# (computation takes around 40 minutes)

model.spli.RE.cal <-trivPenal(Surv(time0, time1, new.lesions) ~ cluster(id)
+ age + treatment + who.PS +  terminal(state),
formula.terminalEvent =~ age + treatment + who.PS + prev.resection, 
tumor.size ~ year * treatment + age + who.PS, data = colorectal, 
data.Longi = colorectalLongi, random = c("1", "year"), id = "id", 
link = "Random-effects", left.censoring = -3.33, recurrentAG = TRUE,
n.knots = 6, kappa=c(0.01, 2), method.GH="Pseudo-adaptive",
n.nodes=7, init.B = c(-0.07, -0.13, -0.16, -0.17, 0.42, #recurrent events covariates
-0.23, -0.1, -0.09, -0.12, 0.8, -0.23, #terminal event covariates
3.02, -0.30, 0.05, -0.63, -0.02, -0.29, 0.11, 0.74)) #biomarker covariates


##---Surrogacy evaluation based on ganerated data with a combination 
##of Monte Carlo and classical Gaussian Hermite integration.
## (Computation takes around 5 minutes)

# Generation of data to use 
data.sim <- jointSurrSimul(n.obs=600, n.trial = 30,cens.adm=549.24, 
         alpha = 1.5, theta = 3.5, gamma = 2.5, zeta = 1, sigma.s = 0.7, 
         sigma.t = 0.7, rsqrt = 0.8, betas = -1.25, betat = -1.25, 
         full.data = 0, random.generator = 1, seed = 0, nb.reject.data = 0)

# Joint surrogate model estimation
joint.surro.sim.MCGH <- jointSurroPenal(data = data.sim, int.method = 2, 
                   nb.mc = 300, nb.gh = 20)


## End(Not run)

---

Fit an Additive Frailty model using a semiparametric penalized likelihood estimation or a parametric estimation

Description

Fit an additive frailty model using a semiparametric penalized likelihood estimation or a parametric estimation. The main issue in a meta-analysis study is how to take into account the heterogeneity between trials and between the treatment effects across trials. Additive models are proportional hazard model with two correlated random trial effects that act either multiplicatively on the hazard function or in interaction with the treatment, which allows studying for instance meta-analysis or multicentric datasets. Right-censored data are allowed, but not the left-truncated data. A stratified analysis is possible (maximum number of strata = 2). This approach is different from the shared frailty models.

In an additive model, the hazard function for the j^th subject in the i^th trial with random trial effect u_i as well as the random treatment-by-trial interaction v_i is:

where $\lambda$₀(0) is the baseline hazard function, $\beta$_k the fixed effect associated to the covariate X_ijk (k=1,..,p), $\beta$₁ is the treatment effect and X_ij1 the treatment variable. $\rho$ is the corresponding correlation coefficient for the two frailty terms.

Usage

additivePenal(formula, data, correlation = FALSE, recurrentAG =
FALSE, cross.validation = FALSE, n.knots, kappa, maxit = 350, hazard =
"Splines", nb.int, LIMparam = 1e-4, LIMlogl = 1e-4, 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. The slope() function is required. Interactions are possible using * or :.
data	a 'data.frame' with the variables used in 'formula'.
correlation	Logical value. Are the random effects correlated? If so, the correlation coefficient is estimated. The default is FALSE.
recurrentAG	Always FALSE for additive models (left-truncated data are not allowed).
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 two strata. The default is FALSE.
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. Number of knots must be between 4 and 20. (See Note)
kappa	positive smoothing parameter in the penalized likelihood estimation. In a stratified additive model, this argument must be a vector with kappas for both strata. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit. 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 Marquardtt algorithm. Default is 350
hazard	Type of hazard functions: "Splines" for semiparametric hazard functions with the penalized likelihood estimation, "Piecewise-per" for piecewise constant hazards functions using percentile, "Piecewise-equi" for piecewise constant hazard functions using equidistant intervals, "Weibull" for parametric Weibull functions. Default is "Splines".
nb.int	Number of intervals (between 1 and 20) for the parametric hazard functions ("Piecewise-per", "Piecewise-equi").
LIMparam	Convergence threshold of the Marquardt algorithm for the parameters (see Details), $10^{-4}$ by default.
LIMlogl	Convergence threshold of the Marquardt algorithm for the log-likelihood (see Details), $10^{-4}$ 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

The estimated parameter are obtained by maximizing the penalized log-likelihood or by a simple log-likelihood (in the parametric case) using the robust Marquardtt algorithm (Marquardtt,1963). The parameters are initialized with values obtained with Cox proportional hazard model. The iterations are stopped when the difference between two consecutive loglikelhoods was small $(<10^{-4})$, the estimated coefficients were stable (consecutive values $(<10^{-4})$, and the gradient small enough $(<10^{-3})$. To be sure of having a positive function at all stages of the algorithm, the spline coefficients were reparametrized to be positive at each stage. The variance space of the two random effects is reduced, so the variances are positive, and the correlation coefficient values are constrained to be between -1 and 1. The marginal log-likelihood depends on integrations that are approximated by using the Laplace integration technique with a first order approximation. The smoothing parameter can be fixed or estimated by maximizing likelihood cross-validation criterion. The usual squared Wald statistic was modified to a mixture of two $\chi^2$ distribution to get significance test for the variance of the random effects.

INITIAL VALUES

The splines and the regression coefficients are initialized to 0.1. An adjusted Cox model is fitted, it provides new initial values for the splines coefficients and the regression coefficients. The variances of the frailties are initialized to 0.1. Then an additive frailty model with independent frailties is fitted. At last, an additive frailty model with correlated frailties is fitted.

Value

An additive model or more generally an object of class 'additivePenal'. Methods defined for 'additivePenal' objects are provided for print, plot and summary.

b	sequence of the corresponding estimation of the splines coefficients, the random effects variances and the regression coefficients.
call	The code used for fitting the model.
coef	the regression coefficients.
cov	covariance between the two frailty terms $(\bold{cov}(u_i,v_i))$
cross.Val	Logical value. Is cross validation procedure used for estimating the smoothing parameters in the penalized likelihood estimation?
correlation	Logical value. Are the random effects correlated?
DoF	degrees of freedom associated with the "kappa".
formula	the formula part of the code used for the model.
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.
n.strat	number of stratum.
rho	the corresponding correlation coefficient for the two frailty terms.
sigma2	Variance for the random intercept (the random effect associated to the baseline hazard functions).
tau2	Variance for the random slope (the random effect associated to the treatment effect across trials).
varH	the variance matrix of all parameters before positivity constraint transformation (Sigma2, Tau2, the regression coefficients and the spline coefficients). Then after, the delta method is needed to obtain the estimated variance parameters.
varHIH	the robust estimation of the variance matrix of all parameters (Sigma2, Tau2, the regression coefficients and the spline coefficients).
varSigma2	The variance of the estimates of "sigma2".
varTau2	The variance of the estimates of "tau2".
varcov	Variance of the estimates of "cov".
x	matrix of times where both survival and hazard functions 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.
type.of.hazard	Type of hazard functions (0:"Splines", "1:Piecewise", "2:Weibull").
type.of.Piecewise	Type of Piecewise hazard functions (1:"percentile", 0:"equidistant").
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.
noVar	indicator of explanatory variable.
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.
frailty.pred	empirical Bayes prediction of the first frailty term.
frailty.pred2	empirical Bayes prediction of the second frailty term.
linear.pred	linear predictor: uses simply "Beta'X + u_i + v_i * X_1" in the additive 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.

Note

"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

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.

V. Rondeau, S. Michiels, B. Liquet, and J. P. Pignon (2008). Investigating trial and treatment heterogeneity in an individual patient data meta-analysis of survival data by mean of the maximum penalized likelihood approach. Statistics in Medecine, 27, 1894-1910.

See Also

slope

Examples

###--- Additive model with 1 covariate ---###

data(dataAdditive)

modAdd <- additivePenal(Surv(t1,t2,event)~cluster(group)+
var1+slope(var1),correlation=TRUE,data=dataAdditive,
n.knots=8,kappa=10000)

#-- Var1 is boolean as a treatment variable

---

Breast Cosmesis Data

Description

The often used data set for interval-censored data, described and given in full in Finkelstein and Wolfe (1985). It involves 94 breast cancer patients who were randomized to either radiation therapy with chemotherapy or radiation therapy alone. The outcome is time until the onset of breast retraction which is interval-censored between the last clinic visit before the event was observed and the first visit when the event was observed. Patients without breast retraction were right-censored.

Usage

data(bcos)

Format

A data frame with 94 observations and 3 variables:

left

left end point of the breast retraction interval

right

right end point of the breast retraction interval

treatment

type of treatment received

Source

Finkelstein, D.M. and Wolfe, R.A. (1985). A semiparametric model for regression analysis of interval-censored failure time data. Biometrics 41, 731-740.

---

Identify clusters

Description

This is a special function used in the context of the models for grouped data. It identifies correlated groups of observations defined by using 'cluster' function, and is used of 'frailtyPenal' formula for fitting univariate and joint models.

Usage

cluster(x)

Arguments

x	A character, factor, or numeric variable which is supposed to indicate the variable group

Value

x	A variable identified as a cluster

See Also

frailtyPenal

Examples

data(readmission)
modSha <- frailtyPenal(Surv(time,event)~as.factor(dukes)+cluster(id),
n.knots=10,kappa=10000,data=readmission,hazard="Splines")

print(modSha)

---

Concordance measures in shared frailty and Cox proportional hazard models

Description

Compute concordance probability estimation for Cox proportional hazard or shared frailty models in case of grouped data (Mauguen et al. 2012). Concordance is given at different levels of comparison, taking into account the cluster membership: between-groups, within-groups and an overall measure, being a weighted average of the previous two. Can also compute the c-index (Harrell et al. 1996) at these three levels. It is possible to exclude tied pairs from concordance estimation (otherwise, account for 1/2).

Usage

Cmeasures(fitc, ties = 1, marginal = 0, cindex = 0, Nboot = 0, tau = 0,
data.val)

Arguments

fitc	A frailtyPenal object, for a shared frailty model. If the fit is a Cox model, no clustering membership is taken into account and only marginal concordance probability estimation is provided. Only an overall measure is given, where all patients are compared two by two. If a counting process formulation is used to performed the fit, with 't.start' and 't.stop', the gap-times (t.stop-t.start) are used in the concordance estimation.
ties	Indicates if the tied pairs on prediction value must be included (ties=1) or excluded (ties=0) from the concordance estimation. Default is ties=1. When included, tied pairs account for 1/2 in the concordance.
marginal	Indicates if the concordance based on marginal predictions must be given (marginal=1) in addition to conditional ones or not (marginal=0). Marginal predictions do not include the frailty estimation in the linear predictor computation: uses "‘Beta’X"' instead of "Beta'X + log z_i". Default is marginal=0.
cindex	Indicates if the c-index (Harrell et al. 1996) must be computed (cindex=1) in addition to the concordance probability estimation or not (cindex=0). C-index is also given at the three comparison levels (between, within and overall). Default is cindex=0.
Nboot	Number of bootstrap resamplings to compute standard-error of the concordances measures, as well as a percentile 95% confidence interval. Nboot=0 indicates no bootstrap procedure. Maximum admitted is 1000. Minimum admitted is 2. Default is 0. Resampling is done at the group level. If Cox model is used, resampling is done at individual level.
tau	Time used to limit the interval on which the concordance is estimated. Note that the survival function for the underlying censoring time distribution needs to be positive at tau. If tau=0, the maximum of the observed event times is used. Default is tau=0.
data.val	A dataframe. It is possible to specify a different dataset than the one used in the model input in the argument 'fitc'. This new dataset will be a validation population and the function will compute new concordance measures from the parameters estimated on the development population. In this case for conditional measures, the frailties are a posteriori predicted. The two datasets must have the same covariates with the same coding without missing data.

Value

call	The shared frailty model evaluated.
Frailty	Logical value. Was model with frailties fitted.
frequencies	Numbers of patients, events and groups used to fit the model.
Npairs	Number of pairs of subjects, between-groups, within-groups and over all the population. If cindex=1, number of comparable (useable) pairs also available.
Nboot	Number of bootstrap resamplings required.
ties	A binary, indicating if the tied pairs on prediction were used to compute the concordance.
CPEcond	Values of Gonen & Heller's measure (conditional). If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).
Cunocond	Values of Uno's measure (conditional). If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).
marginal	A binary, indicating if the marginal values were computed.
CPEmarg	Values of Gonen & Heller's measure (marginal), if marginal=1. If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).
Cunomarg	Values of Uno's measure (marginal), if marginal=1. If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).
cindex	A binary, indicating if the c-indexes were computed.
cindexcond	Values of the C-index of Harrell (conditional). If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).
cindexmarg	Values of the C-index of Harrell (marginal), if marginal=1. If Nboot>0, give SE, the standard-error of the parameters evaluated by bootstrap, IC.low and IC.high, the lower and upper bounds of the percentile confidence interval evaluated by bootstrap (2.5% and 97.5% percentiles).

References

Mauguen, A., Collette, S., Pignon, J. P. and Rondeau, V. (2013). Concordance measures in shared frailty models: application to clustered data in cancer prognosis. Statistics in Medicine 32, 27, 4803-4820

Harrell, F.E. et al. (1996). Tutorial in biostatistics: multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15, 361-387.

Gonen, M., Heller, G. (2005). Concordance probability and discriminatory power in proportional hazards regression. Biometrika 92, 965-970.

See Also

print.Cmeasures,frailtyPenal

Examples

#-- load data
data(readmission)

#-- a frailtypenal fit
fit <- frailtyPenal(Surv(time,event)~cluster(id)+dukes+
charlson+chemo,data=readmission,cross.validation=FALSE,
n.knots=10,kappa=1,hazard="Splines")

#-- a Cmeasures call
fit.Cmeasures <- Cmeasures(fit)
fit.Cmeasures.noties <- Cmeasures(fit, ties=0)
fit.Cmeasures.marginal <- Cmeasures(fit, marginal=1)
fit.Cmeasures.cindex <- Cmeasures(fit, cindex=1)

#-- a short summary
fit.Cmeasures
fit.Cmeasures.noties
fit.Cmeasures.marginal
fit.Cmeasures.cindex

---

Follow-up of metastatic colorectal cancer patients: times of new lesions appearance and death

Description

Randomly chosen 150 patients from the follow-up of the FFCD 2000-05 multicenter phase III clinical trial originally including 410 patients with metastatic colorectal cancer randomized into two therapeutic strategies: combination and sequential. The dataset contains times of observed appearances of new lesions censored by a terminal event (death or right-censoring) with baseline characteristics (treatment arm, age, WHO performance status and previous resection).

Usage

data(colorectal)

Format

This data frame contains the following columns:

id

identification of each subject. Repeated for each recurrence

time0

start of interval (0 or previous recurrence time)

time1

recurrence or censoring time

new.lesions

Appearance of new lesions status. 0: censsored or no event, 1: new lesions

treatment

To which treatment arm a patient was allocated? 1: sequential (S); 2: combination (C)

age

Age at baseline: 1: <50 years, 2: 50-69 years, 3: >69 years

who.PS

WHO performance status at baseline: 1: status 0, 2: status 1, 3: status 2

prev.resection

Previous resection of the primate tumor? 0: No, 1: Yes

state

death indicator. 0: alive, 1: dead

gap.time

interocurrence time or censoring time

Note

We thank the Federation Francophone de Cancerologie Digestive and Gustave Roussy for sharing the data of the FFCD 2000-05 trial supported by an unrestricted Grant from Sanofi.

References

M. Ducreux, D. Malka, J. Mendiboure, P.-L. Etienne, P. Texereau, D. Auby, P. Rougier, M. Gasmi, M. Castaing, M. Abbas, P. Michel, D. Gargot, A. Azzedine, C. Lombard- Bohas, P. Geoffroy, B. Denis, J.-P., Pignon, L.,Bedenne, and O. Bouche (2011). Sequential versus combination chemotherapy for the treatment of advanced colorectal cancer (FFCD 2000-05): an open-label, randomised, phase 3 trial. The Lancet Oncology 12, 1032-44.

---

Follow-up of metastatic colorectal cancer patients : longitudinal measurements of tumor size

Description

Randomly chosen 150 patients from the follow-up of the FFCD 2000-05 multicenter phase III clinical trial originally including 410 patients with metastatic colorectal cancer randomized into two therapeutic strategies: combination and sequential. The dataset contains measurements of tumor size (left-censored sums of the longest diameters of target lesions; transformed using Box-Cox) with baseline characteristics(treatment arm, age, WHO performance status and previous resection).

Usage

data(colorectalLongi)

Format

This data frame contains the following columns:

id

identification of each subject. Repeated for each recurrence

year

time of visit counted in years from baseline

tumor.size

Individual longitudinal measurement of transformed (Box-Cox with parameter 0.3) sums of the longest diameters, left-censored due to a detection limit (threshold $s=-3.33$).

treatment

To which treatment arm a patient was allocated? 1: sequential (S); 2: combination (C)

age

Age at baseline: 1: <50 years, 2: 50-69 years, 3: >69 years

who.PS

WHO performance status at baseline: 1: status 0, 2: status 1, 3: status 2

prev.resection

Previous resection of the primate tumor? 0: No, 1: Yes

Note

We thank the Federation Francophone de Cancerologie Digestive and Gustave Roussy for sharing the data of the FFCD 2000-05 trial supported by an unrestricted Grant from Sanofi.

References

Ducreux, M., Malka, D., Mendiboure, J., Etienne, P.-L., Texereau, P., Auby, D., Rougier, P., Gasmi, M., Castaing, M., Abbas, M., Michel, P., Gargot, D., Azzedine, A., Lombard- Bohas, C., Geoffroy, P., Denis, B., Pignon, J.-P., Bedenne, L., and Bouche, O. (2011). Sequential versus combination chemotherapy for the treatment of advanced colorectal cancer (FFCD 2000-05): an open-label, randomised, phase 3 trial. The Lancet Oncology 12, 1032-44.

---

Simulated data as a gathering of clinical trials databases

Description

This contains simulated samples of 100 clusters with 100 subjects in each cluster, like a gathering of clinical trials databases. Two correlated centred gaussian random effects are generated with the same variance fixed at 0.3 and the covariance at -0.2. The regression coefficient $\beta$ is fixed at -0.11. The percentage of right-censored data is around 30 percent which are generated from a uniform distribution on [1,150]. The baseline hazard function is considered as a simple Weibull.

Usage

data(dataAdditive)

Format

This data frame contains the following columns:

group

identification variable

t1

start of interval (=0, because left-truncated data are not allowed)

t2

end of interval (death or censoring time)

event

censoring status (0:alive, 1:death, as acensoring indicator

var1

dichotomous covariate (=0 or 1,as a treatment variable)

var2

dichotomous covariate (=0 or 1,as a treatment variable)

Source

V. Rondeau, S. Michiels, B.Liquet, and J.P. Pignon (2008). Investigating trial and treatment heterogeneity in an individual patient data meta-analysis of survival data by mean of the maximum penalized likelihood approach. Statistics in Medecine, 27, 1894-1910.

---

Simulated data for two types of recurrent events and a terminal event

Description

This contains a simulated sample of of 800 subjects and 1652 observations. This dataset can be used to illustrate how to fit a joint multivariate frailty model. Two gaussian correlated random effects were generated with mean 0, variances 0.5 and a correlation coefficient equals to 0.5. The coefficients $\alpha_1$ and $\alpha_2$ were fixed to 1. The three baseline hazard functions followed a Weibull distribution and right censoring was fixed at 5.

Usage

data(dataMultiv)

Format

This data frame contains the following columns:

PATIENT

identification of patient

obs

number of observation for a patient

TIME0

start of interval

TIME1

end of interval (death or censoring time)

INDICREC

recurrent of type 1 status (0:no, 1:yes)

INDICMETA

recurrent of type 2 status (0:no, 1:yes)

INDICDEATH

censoring status (0:alive, 1:death)

v1

dichotomous covariate (0,1)

v2

dichotomous covariate (0,1)

v3

dichotomous covariate (0,1)

TIMEGAP

time to event

---

Simulated data for recurrent events and a terminal event with weigths using nested case-control design

Description

This contains a simulated sample of of 819 subjects and 1510 observations. This dataset can be used to illustrate how to fit a joint frailty model for data from nested case-control studies.

Usage

data(dataNCC)

Format

This data frame contains the following columns:

id

identification of patient

cov1

dichotomous covariate (0,1)

cov2

dichotomous covariate (0,1)

t.start

start of interval

t.stop

end of interval (death or censoring time)

gaptime

time to event

event

recurrent event status (0:no, 1:yes)

deathdays

time of terminal event (death or right-censoring)

death

censoring status (0:alive, 1:death)

ncc.wts

weights for NCC design

---

Simulated data with two levels of grouping

Description

This contains a simulated sample of 400 observations which allow establishing 20 clusters with 4 subgroups and 5 subjects in each subgroup, in order to obtain two levels of grouping. This data set is useful to illustrate how to fit a nested model. Two independent gamma frailty parameters with a variance fixed at 0.1 for the cluster effect and at 0.5 for the subcluster effect were generated. Independent survival times were generated from a simple Weibull baseline risk function. The percentage of censoring data was around 30 per cent. The right-censoring variables were generated from a uniform distribution on [1,36] and a left-truncating variable was generated with a uniform distribution on [0,10]. Observations were included only if the survival time is greater than the truncated time.

Usage

data(dataNested)

Format

This data frame contains the following columns:

group

group identification variable

subgroup

subgroup identification variable

t1

start of interval (0 or truncated time)

t2

end of interval (death or censoring time)

event

censoring status (0: alive, 1: death)

cov1

dichotomous covariate (0,1)

cov2

dichotomous covariate (0,1)

Source

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

---

Advanced Ovarian Cancer dataset

Description

This dataset combines the data that were collected in four double-blind randomized clinical trials in advanced ovarian cancer. In these trials, the objective was to examine the efficacy of cyclophosphamide plus cisplatin (CP) versus cyclophosphamide plus adriamycin plus cisplatin (CAP) to treat advanced ovarian cancer. The candidate surrogate endpoint S is progression-free survival time, defined as the time (in years) from randomization to clinical progression of the disease or death. The true endpoint T is survival time, defined as the time (in years) from randomization to death of any cause

Usage

data(dataOvarian)

Format

This data frame contains the following columns:

patientID

The identification number of a patient

trialID

The center in which a patient was treated

trt

The treatment indicator, coded as 0 = cyclophosphamide plus cisplatin (CP) and 1 = cyclophosphamide plus adriamycin plus cisplatin(CAP)

timeS

The candidate surrogate (progression-free survival)

statusS

Censoring indicator for for Progression-free survival

timeT

The true endpoint (survival time)

statusT

Censoring indicator for survival time

Source

Ovarian cancer Meta-Analysis Project (1991). Cyclophosphamide plus cisplatin plus adriamycin versus Cyclophosphamide, doxorubicin, and cisplatin chemotherapy of ovarian carcinoma: A meta-analysis. Classic Papers and Current Comments, 3, 237-234.

---

Difference of Expected Prognostic Observed Cross-Entropy (EPOCE) estimators and its 95% tracking interval between two joint models.

Description

This function computes the difference of two EPOCE estimates (CVPOL and MPOL) and its 95% tracking interval between two joint models estimated using frailtyPenal, longiPenal or trivPenal. Difference in CVPOL is computed when the EPOCE was previously estimated on the same dataset as used for estimation (using an approximated cross-validation), and difference in MPOL is computed when the EPOCE was previously estimated on an external dataset.

Usage

Diffepoce(epoce1, epoce2)

Arguments

epoce1	a first object inheriting from class epoce.
epoce2	a second object inheriting from class epoce.

Details

From the EPOCE estimates and the individual contributions to the prognostic observed log-likelihood obtained with epoce function on the same dataset from two different estimated joint models, the difference of CVPOL (or MPOL) and its 95% tracking interval is computed. The 95% tracking interval is : Delta(MPOL) +/- qnorm(0.975)*sqrt(VARIANCE) for an external dataset Delta(CVPOL) +/- qnorm(0.975)*sqrt(VARIANCE) for the dataset used in frailtyPenal, longiPenal or trivPenal where Delta(CVPOL) (or Delta(MPOL)) is the difference of CVPOL (or MPOL) of the two joint models, and VARIANCE is the empirical variance of the difference of individuals contributions to the prognostic observed log-likelihoods of the two joint models.

The estimators of EPOCE from arguments epoce1 and epoce2 must have been computed on the same dataset and with the pred.times.

Value

new.data	a boolean which is FALSE if computation is done on the same data as for estimation, and TRUE otherwise
pred.times	time or vector of times used in the function
DEPOCE	the difference between the two MPOL or CVPOL for each time
TIinf	lower confidence band for the difference
TIsup	upper confidence band for the difference

References

D. Commenges, B. Liquet, C. Proust-Lima (2012). Choice of prognostic estimators in joint models by estimating differences of expected conditional Kullback-Leibler risks. Biometrics, 68(2), 380-387.

Examples

## Not run: 

#Example for joint frailty models
data(readmission)

# first joint frailty model
joint1 <- frailtyPenal(Surv(t.start,t.stop,event)~ cluster(id) +
  dukes + charlson + sex + chemo + terminal(death),
  formula.terminalEvent = ~ dukes + charlson + sex + chemo ,
  data = readmission, n.knots = 8, kappa = c(2.11e+08,9.53e+11),
  recurrentAG=TRUE)

# second joint frailty model without dukes nor charlson as covariates
joint2 <- frailtyPenal(Surv(t.start,t.stop,event)~ cluster(id) +
  sex + chemo + terminal(death),
  formula.terminalEvent = ~ sex + chemo ,
  data = readmission, n.knots = 8, kappa = c(2.11e+08,9.53e+11),
  recurrentAG=TRUE)

temps <- c(200,500,800,1100)

# computation of estimators of EPOCE for the two models
epoce1 <- epoce(joint1,temps)
epoce2 <- epoce(joint2,temps)

# computation of the difference
diff <- Diffepoce(epoce1,epoce2)

print(diff)
plot(diff)


#Example for joint models with a biomarker
data(colorectal)
data(colorectalLongi)

# Survival data preparation - only terminal events 
colorectalSurv <- subset(colorectal, new.lesions == 0)

# first joint model for a biomarker and a terminal event
modLongi <- longiPenal(Surv(time0, time1, state) ~ age +
treatment + who.PS, tumor.size ~  year*treatment + age +
who.PS, colorectalSurv, data.Longi =colorectalLongi,
random = c("1", "year"),  id = "id", link = "Random-effects", 
left.censoring = -3.33, hazard = "Weibull", 
method.GH = "Pseudo-adaptive")

# second joint model for a biomarker, recurrent events and a terminal event
# (computation takes around 30 minutes)
modTriv <- model.weib.RE.gap <-trivPenal(Surv(gap.time, new.lesions) 
~ cluster(id) + age + treatment + who.PS + prev.resection + terminal(state),
formula.terminalEvent =~ age + treatment + who.PS + prev.resection, 
tumor.size ~ year * treatment + age + who.PS, data = colorectal,
data.Longi = colorectalLongi, random = c("1", "year"), id = "id", 
link = "Random-effects", left.censoring = -3.33, recurrentAG = FALSE,
hazard = "Weibull", method.GH="Pseudo-adaptive", n.nodes=7)

time <- c(1, 1.5, 2, 2.5)

# computation of estimators of EPOCE for the two models
epoce1 <- epoce(modLongi, time)
# (computation takes around 10 minutes)
epoce2 <- epoce(modTriv, time)


# computation of the difference
diff <- Diffepoce(epoce1, epoce2)

print(diff)
plot(diff)

## End(Not run)

---

Estimators of the Expected Prognostic Observed Cross-Entropy (EPOCE) for evaluating predictive accuracy of joint models.

Description

This function computes estimators of the Expected Prognostic Observed Cross-Entropy (EPOCE) for evaluating the predictive accuracy of joint models using frailtyPenal, longiPenal, trivPenal or trivPenalNL. On the same data as used for estimation of the joint model, this function computes both the Mean Prognosis Observed Loss (MPOL) and the Cross-Validated Prognosis Observed Loss (CVPOL), two estimators of EPOCE. The latter corrects the MPOL estimate for over-optimism by approximated cross-validation. On external, this function only computes MPOL.

Usage

epoce(fit, pred.times, newdata = NULL, newdata.Longi = NULL)

Arguments

fit	A jointPenal, longiPenal, trivPenal or trivPenalNL object.
pred.times	Time or vector of times to compute epoce.
newdata	Optional. In case of joint models obtained with frailtyPenal, trivPenal or trivPenalNL. For models inheriting from trivPenal or trivPenalNL class, if newdata is given, newdata.Longi must be given as well. When missing, the data used for estimating the fit are used, and CVPOL and MPOL are computed (internal validation). When newdata is specified, only MPOL is computed on this new dataset (external validation). The new dataset and the dataset used in the estimation must have the same covariates with the same coding without missing data.
newdata.Longi	Optional. In case of joint models obtained with longiPenal, trivPenal or trivPenalNL. For models inheriting from longiPenal, if the newdata.Longi is given, newdata must be NULL, but for models from trivPenal or trivPenalNL class, if newdata.Longi is given, newdata must be provided as well. The two datasets newdata and newdata.Longi must include the information concerning the same patients with the same characteristics and the appropriate data on follow up (recurrences for newdata and longitudinal measurements for newdata.Longi).

Value

data	name of the data used to compute epoce
new.data	a boolean which is FALSE if computation is done on the same data as for estimation, and TRUE otherwise
pred.times	time or vector of times used in the function
mpol	values of MPOL for each pred.times
cvpol	values of CVPOL for each pred.times
IndivContrib	all the contributions to the log-likelihood for each pred.times
AtRisk	number of subject still at risk for each pred.times

References

D. Commenges, B. Liquet, C. Proust-Lima (2012). Choice of prognostic estimators in joint models by estimating differences of expected conditional Kullback-Leibler risks. Biometrics, 68(2), 380-387.

Examples

## Not run: 

########################################
#### EPOCE on a joint frailty model ####
########################################

data(readmission)

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

# computation on the same dataset
temps <- c(200,500,800,1100)
epoce <- epoce(modJoint.gap,temps)

print(epoce)
plot(epoce,type = "cvpol")

# computation on a new dataset
# here a sample of readmission with the first 50 subjects
s <- readmission[1:100,]
epoce <- epoce(modJoint.gap,temps,newdata=s)

print(epoce)
plot(epoce,type = "mpol")

#################################################
#### EPOCE on a joint  model for a biomarker ####
#########   and a terminal event  ###############
#################################################

data(colorectal)
data(colorectalLongi)

# Survival data preparation - only terminal events 
colorectalSurv <- subset(colorectal, new.lesions == 0)

modLongi <- longiPenal(Surv(time0, time1, state) ~ age +
treatment + who.PS, tumor.size ~  year*treatment + age +
who.PS, colorectalSurv, data.Longi =colorectalLongi,
random = c("1", "year"),  id = "id", link = "Random-effects", 
left.censoring = -3.33, hazard = "Weibull", 
method.GH = "Pseudo-adaptive")

# computation on the same dataset
time <- c(1, 1.5, 2, 2.5)
epoce <- epoce(modLongi,time)

print(epoce)
plot(epoce, type = "cvpol")

# computation on a new dataset
# here a sample of colorectal data with the first 50 subjects
s <-  subset(colorectal, new.lesions == 0 & id%in%1:50)
s.Longi <- subset(colorectalLongi, id%in%1:50)
epoce <- epoce(modLongi, time, newdata = s, newdata.Longi = s.Longi)

print(epoce)
plot(epoce, type = "mpol")


###################################################
#### EPOCE on a joint model for a biomarker, ######
#### recurrent events and a terminal event   ######
###################################################

data(colorectal)
data(colorectalLongi)

# Linear model for the biomarker
# (computation takes around 30 minutes)
model.trivPenalNL <-trivPenal(Surv(gap.time, new.lesions) ~ cluster(id)
+ age + treatment + who.PS + prev.resection + terminal(state),
formula.terminalEvent =~ age + treatment + who.PS + prev.resection, 
tumor.size ~ year * treatment + age + who.PS, data = colorectal,
data.Longi = colorectalLongi, random = c("1", "year"), id = "id", 
link = "Random-effects", left.censoring = -3.33, recurrentAG = FALSE,
hazard = "Weibull", method.GH="Pseudo-adaptive", n.nodes=7)

# computation on the same dataset
time <- c(1, 1.5, 2, 2.5)

# (computation takes around 10 minutes)
epoce <- epoce(model.trivPenalNL,time)
print(epoce)
plot(epoce, type = "cvpol")

# computation on a new dataset
# here a sample of colorectal data with the first 100 subjects
s <-  subset(colorectal,  id%in%1:100)
s.Longi <- subset(colorectalLongi, id%in%1:100)
# (computation takes around 10 minutes)
epoce <- epoce(model.trivPenalNL, time, newdata = s, newdata.Longi = s.Longi)

print(epoce)
plot(epoce, type = "mpol")



# Non-linear model for the biomarker

# No information on dose - creation of a dummy variable 
colorectalLongi$dose <- 1

# (computation can take around 40 minutes)
model.trivPenalNL <- trivPenalNL(Surv(time0, time1, new.lesions) ~ cluster(id) + age + treatment
 + terminal(state), formula.terminalEvent =~ age + treatment, biomarker = "tumor.size",
 formula.KG ~ 1, formula.KD ~ treatment, dose = "dose", time.biomarker = "year", 
 data = colorectal, data.Longi =colorectalLongi, random = c("y0", "KG"), id = "id", 
 init.B = c(-0.22, -0.16, -0.35, -0.19, 0.04, -0.41, 0.23), init.Alpha = 1.86,
 init.Eta = c(0.5, 0.57, 0.5, 2.34), init.Biomarker = c(1.24, 0.81, 1.07, -1.53),
 recurrentAG = TRUE, n.knots = 5, kappa = c(0.01, 2), method.GH = "Pseudo-adaptive")

# computation on the same dataset
time <- c(1, 1.5, 2, 2.5)

epoce <- epoce(model.trivPenalNL, time)



## End(Not run)

---

Identify event2 indicator

Description

This is a special function used in the context of multivariate frailty model with two types of recurrent events and a terminal event (e.g., censoring variable related to both recurrent events). It contains the indicator of the recurrent event of type 2, normally 0=no event, 1=event, and is used on the right hand side of a formula of a 'multivPenal' object. Using event2() in a formula implies that a multivariate frailty model for two types of recurrent events and a terminal event is fitted.

Usage

event2(x)

Arguments

x	A numeric variable but should be a boolean which equals 1 if the subject has experienced an event of type 2 and 0 if not.

Value

x	an indicator for an event of type 2

See Also

multivPenal

---

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)

---

Advanced Gastric Cancer dataset

Description

This meta-analysis was carried out by the GASTRIC (Global Advanced/Adjuvant Stomach Tumor Research international Collaboration) group, using individual data on patients with curatively resected gastric cancer. Data from all published randomized trials, with a patient recruitment end date before 2004, and comparing adjuvant chemotherapy with surgery alone for resectable gastric cancers, were searched electronically. The candidate surrogate endpoint S was Desease-free survival time, defined as the time (in days) to relapse, second cancer or dead from any cause. The true endpoint T was the overall survival time, defined as the time (in days) from randomization to death of any cause or to the last follow-up.

Usage

data(gastadj)

Format

This data frame contains the following columns:

trialID

The trial in which the patient was treated

patientID

The identification number of a patient

trt

The treatment indicator, coded as 0 = Control and 1 = Experimental

timeS

The candidate surrogate (progression-free survival in days)

statusS

Censoring indicator for for Progression-free survival (0 = alive and progression-free, 1 = with progression or dead)

timeT

The true endpoint (overall survival time in days)

statusT

Censoring indicator for survival time (0 = alive, 1 = dead)

Source

Oba K, Paoletti X, Alberts S, Bang YJ, Benedetti J, Bleiberg H, Catalona P, Lordick F, Michiels S, Morita A, Okashi Y, Pignon JP, Rougier P, Sasako M, Sakamoto J, Sargent D, Shitara K, Van Cutsem E, Buyse M, Burzykowski T on behalf of the GASTRIC group (2013). Disease-Free Survival as a Surrogate for Overall Survival in Adjuvant Trials of Gastric Cancer: A Meta-Analysis. JNCI: Journal of the National Cancer Institute;105(21):1600-1607

---

Fit a Shared or a Joint Frailty Generalized Survival Model

Description

I. SHARED FRAILTY GENERALIZED SURVIVAL MODELS

Fit a gamma Shared Frailty Generalized Survival Model using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) are allowed. It allows to obtain a parametric or flexible semi-parametric smooth hazard and survival functions.

Each frailty term $u$_i is assumed to act multiplicatively on the hazard function, and to be drawn from a Gamma distribution with unit mean and variance $\theta$. Conditional on the frailty term, the hazard function for the $j$^th subject in the $i$^th group is then expressed by

where $x$_ij is a collection of baseline covariates, $\xi$ is a vector of parameters, and $\lambda$_ij ($t$ | $x$_ij ; $\xi$) is the hazard function for an average value of the frailty. The associated conditional survival function writes

where $S$_ij ($t$ | $x$_ij ; $\xi$) designates the survival function for an average value of the frailty.

Following Liu et al. (2017, 2018), the latter function is expressed in terms of a link function $g$(.) and a linear predictor $\eta$_ij ($t$, $x$_ij; $\xi$) such that $g$[$S$_ij ($t$ | $x$_ij ; $\xi$)] = $\eta$_ij ($t$, $x$_ij; $\xi$), i.e. $S$_ij ($t$ | $x$_ij ; $\xi$) = $h$[$\eta$_ij ($t$, $x$_ij; $\xi$)] with $h$() = $g$^-1().

The conditional survival function is finally modeled by

The table below summarizes the most commonly used (inverse) link functions and their associated conditional survival, hazard and cumulative hazard functions. PHM stands for "Proportional Hazards Model", POM for "Proportional Odds Model, PROM for "Probit Model" and AHM for "Additive Hazards Model".

I.(a) Fully parametric case

In the fully parametric case, linear predictors considered are of the form

where $\rho > 0$ is a shape parameter, $\gamma > 0$ a scale parameter, $\beta$ a vector of regression coefficients, and $\xi$ = ($\rho$ ,$\gamma$, $\beta$). With the appropriate link function, such linear parametric predictors make it possible to recover a Weibull baseline survival function for PHMs and AHMs, a log-logistic baseline survival function for POMs, and a log-normal one for PROMs.

I. (b) Flexible semi-parametric case

For PHM and AHM, a more flexible splines-based approach is proposed for modeling the baseline hazard function and time-varying regression coefficients. In this case, conditional on the frailty term $u$_i, the hazard function for the $j$^th subject in the $i$^th group is still expressed by $\lambda$_ij ($t$ | $x$_ij, $u$_i ; $\xi$) = $u$_i$\lambda$_ij ($t$ | $x$_ij ; $\xi$), but we have this time

The smoothness of baseline hazard function $\lambda$₀() is ensured by penalizing the log-likelihood by a term which has large values for rough functions.

Moreover, for parametric and flexible semi-parametric AHMs, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.

II. JOINT FRAILTY GENERALIZED SURVIVAL MODELS

Fit a gamma Joint Frailty Generalized Survival Model for recurrent and terminal events using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) for the recurrent event part are allowed. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring.

This model includes an common patient-specific frailty term $u$_i for the two survival functions which will take into account the unmeasured heterogeneity in the data, associated with unobserved covariates. The frailty term acts differently for the two survival functions ($u$_i for the recurrent survival function and $u$_i^α for the death one). The covariates could be different for the recurrent and terminal event parts.

II.(a) Fully parametric case

For the $j$^th recurrence (j=1,...,n_i) and the $i$^th patient (i=1,...,N), the gamma Joint Frailty Generalized Survival Model for recurrent event survival function $S$_Rij(.) and death survival function $S$_Di(.) is

- $\eta$_Rij (resp. $\eta$_Di) is the linear predictor for the recurrent (resp. terminal) event process. The form of these linear predictors is the same as the one presented in I.(a).

- $h$_R(.) (resp. $h$_D(.)) is the inverse link function associated with recurrent events (resp. terminal event).

- $x$_Rij and $x$_Di are two vectors of baseline covariates associated with recurrent and terminal events.

- $\xi$_R and $\xi$_D are the parameter vectors for recurrent and terminal events.

- $\alpha$ is a parameter allowing more flexibility in the association between recurrent and terminal events processes.

- The random frailties $u$_i are still assumed iid and drown from a $\Gamma$(1/$\theta$,1/$\theta$).

II.(b) Flexible semi-parametric case

If one chooses to fit a PHM or an AHM for recurrent and/or terminal events, a splines-based approach for modeling baseline hazard functions and time-varying regression coefficients is still available. In this approach, the submodel for recurrent events is expressed as $\lambda$_Rij ($t$ | $x$_Rij, $u$_i ; $\xi$_R) = $u$_i$\lambda$_Rij ($t$ | $x$_Rij ; $\xi$_R), where

The submodel for terminal event is expressed as $\lambda$_Di ($t$ | $x$_Di, $u$_i ; $\xi$_D) = $u$_i^α$\lambda$_Di ($t$ | $x$_Di ; $\xi$_D), where

Baseline hazard functions $\lambda$_R0(.) and $\lambda$_D0(.) are estimated using cubic M-splines (of order 4) with positive coefficients, and the time-varying coefficients $\beta$_R(.) and $\beta$_D(.) are estimated using B-splines of order q.

The smoothness of baseline hazard functions is ensured by penalizing the log-likelihood by two terms which has large values for rough functions.

Moreover, if one chooses an AHM for recurrent and/or terminal event submodel, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.

Usage

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

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. Interactions are possible using ' * ' or ' : '.
formula.terminalEvent	Only for joint frailty models: a formula object, only requires terms on the right to indicate which variables are used for 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.
family	Type of Generalized Survival Model to fit. "PH" for a proportional hazards model, "AH" for an additive hazards model, "PO" for a proportional odds model and "probit" for a probit model. A vector of length 2 is expected for joint models (e.g., family=c("PH","PH")).
hazard	Type of hazard functions: "Splines" for semi-parametric hazard functions using equidistant intervals, or "parametric" for parametric distribution functions. In case of family="PH" or family="AH", the "parametric" option corresponds to a Weibull distribution. In case of family="PO" and family="probit", the "parametric" option corresponds to a log-logistic and a log-normal distribution, respectively. So far, the "Splines" option is only available for PH and AH submodels. Default is "Splines".
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 (i.e. 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. The coefficient kappa tunes the intensity of the penalization (the integral of the squared second derivative of hazard function). 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. 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).
betaknots	Number of inner knots used for the B-splines time-varying coefficient estimation. Default is 1. See 'timedep' function for more details.
betaorder	Order of the B-splines used for the time-varying coefficient estimation. Default is cubic B-splines (order=3). See 'timedep' function for more details. Not implemented for Proportional Odds and Probit submodels.
RandDist	Type of random effect distribution: "Gamma" for a gamma distribution, and "LogN" for a log-normal distribution (not implemented yet). Default is "Gamma".
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 Generalized Survival and Shared Frailty Models) or 0.5 (for Generalized Joint Frailty Models).
init.Theta	Initial value for frailty variance.
init.Alpha	Only for Generalized Joint Frailty Models: initial value for parameter alpha.
Alpha	Only for Generalized Joint Frailty Models: input "None" so as to fit a joint model without the parameter alpha.
maxit	Maximum number of iterations for the Marquardt algorithm. Default is 300
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.
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.
cross.validation	Not implemented yet for the generalized settings.
jointGeneral	Not implemented yet for the generalized settings.
nb.int	Not implemented yet for the generalized settings.
initialize	Not implemented yet for the generalized settings.
init.Ksi	Not implemented yet for the generalized settings.
Ksi	Not implemented yet for the generalized settings.
init.Eta	Not implemented yet for the generalized settings.

Details

TYPICAL USES

For a Generalized Survival Model:

GenfrailtyPenal(
formula=Surv(time,event)~var1+var2, 
data, family, \dots)

For a Shared Frailty Generalized Survival Model:

GenfrailtyPenal(
formula=Surv(time,event)~cluster(group)+var1+var2, 
data, family, \dots)

For a Joint Frailty Generalized Survival Model:

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

OPTIMIZATION ALGORITHM

The estimated parameters 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 is small ($<10$^-3), the estimated coefficients are stable (consecutive values $<10$^-3, and the gradient small enough ($<10$^-3). When the frailty variance 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). For Proportional Hazards and Additive Hazards submodels, cubic M-splines of order 4 can be used to estimate the hazard function. In this case, I-splines (integrated M-splines) are used to compute the cumulative hazard function.

The inverse of the Hessian matrix is the variance estimator. 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 integrations in the full log likelihood are evaluated using Gaussian quadrature. Laguerre polynomials with 20 points are used to treat the integrations on $[0,\infty[$.

INITIAL VALUES

In case of a shared frailty model, the splines and the regression coefficients are initialized to 0.1. 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 firstly, an adjusted Cox model to have new initial values for the splines and the regression coefficients. The variance of the frailty term $\theta$ and the association parameter $\alpha$ are initialized to 1. Then, a joint frailty model is fitted.

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.
n	The number of observations used in the fit.
groups	The maximum number of groups used in the fit.
n.events	The number of events observed in the fit.
n.eventsbygrp	A vector of length the number of groups giving the number of observed events in each group.
loglik	The marginal log-likelihood in the parametric case.
loglikPenal	The marginal penalized log-likelihood in the semiparametric case.
coef	The regression coefficients.
varH	The variance matrix of the regression coefficients 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.
varHtotal	The variance matrix of all the 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 the regression coefficients
varHIHtotal	The robust estimation of the variance matrix of all parameters.
x	Matrix of times where the hazard functions are estimated.
xSu	Matrix of times where the survival functions are estimated.
lam	Array (dim=3) of baseline hazard estimates and confidence bands.
surv	Array (dim=3) of baseline survival estimates and confidence bands.
type	Character string specifying the type of censoring, see the Surv function for more details.
n.strat	Number of strata.
n.iter	Number of iterations needed to converge.
median	The value of the median survival and its confidence bands. If there are two strata or more, the first value corresponds to the value for the first strata, etc.
LCV	The approximated likelihood cross-validation criterion in the semiparametric case. With H (resp. Hpen) the hessian matrix of log-likelihood (resp. penalized log-likelihood), EDF = Hpen-1 H the effective degrees of freedom, L($\xi$,$\theta$) the log-likelihood and n the number of observations, $LCV = 1/n x (trace(EDF) - L(\xi,\theta)).$
AIC	The Akaike information Criterion for the parametric case. With p the number of parameters, n the number of observations and L($\xi$,$\theta$) the log-likelihood, $AIC = 1/n x (p - L(\xi,\theta)).$
npar	Number of parameters.
nvar	Number of explanatory variables.
typeof	Indicator of the type of hazard functions computed : 0 for "Splines", 2 for "parametric".
istop	Convergence indicator: 1 if convergence is reached, 2 if convergence is not reached, 3 if the hessian matrix is not positive definite, 4 if a numerical problem has occurred in the likelihood calculation
shape.param	Shape parameter for the parametric hazard function (a Weibull distribution is used for proportional and additive hazards models, a log-logistic distribution is used for proportional odds models, a log-normal distribution is used for probit models).
scale.param	Scale parameter for the parametric hazard function.
Names.data	Name of the dataset.
Frailty	Logical value. Was model with frailties fitted ?
linear.pred	Linear predictor: $\beta$'$X$ in the generalized survival models or $\beta$'$X$ + log($u$i) in the shared frailty generalized survival models.
BetaTpsMat	Matrix of time varying-effects and confidence bands (the first column used for abscissa of times).
nvartimedep	Number of covariates with time-varying effects.
Names.vardep	Name of the covariates with time-varying effects.
EPS	Convergence criteria concerning the parameters, the likelihood and the gradient.
family	Type of Generalized Survival Model fitted (0 for PH, 1 for PO, 2 for probit, 3 for AH).
global_chisq.test	A binary variable equals to 0 when no multivariate Wald is given, 1 otherwise.
beta_p.value	p-values of the Wald test for the estimated 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 smoothing parameter kappa.
kappa	A vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components.
n.knots	Number of knots for estimating the baseline functions in the penalized likelihood estimation.
n.knots.temp	Initial value for the number of knots.
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.
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.

The following components are specific to shared models.

equidistant	Indicator for the intervals used in the spline estimation of baseline hazard functions : 1 for equidistant intervals ; 0 for intervals using percentile (note: equidistant = 2 in case of parametric estimation).
Names.cluster	Cluster names.
theta	Variance of the gamma frailty parameter, i.e. Var($u$i).
varTheta	Variance of parameter theta.
theta_p.value	p-value of the Wald test for the estimated variance of the gamma frailty.

The following components are specific to joint models.

formula	The formula part of the code used for the recurrent events.
formula.terminalEvent	The formula part of the code used for the terminal model.
n.deaths	Number of observed deaths.
n.censored	Number of censored individuals.
theta	Variance of the gamma frailty parameter, i.e. Var($u$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.
nvar	A vector with the number of covariates of each type of hazard function as components.
nvarnotdep	A vector with the number of constant effect 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.
Names.vardep	Name of the covariates with time-varying effects for the recurrent events.
Names.vardepdc	Name of the covariates with time-varying effects for the terminal event.
xR	Matrix of times where both survival and hazard function are estimated for the recurrent event.
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.
nb.gh	Number of nodes for the Gaussian-Hermite quadrature.
nb.gl	Number of nodes for the Gaussian-Laguerre quadrature.
medianR	The value of the median survival for the recurrent events and its confidence bands.
medianD	The value of the median survival for the terminal event and its confidence bands.
names.factor	Names of the "as.factor" variables for the recurrent events.
names.factordc	Names of the "as.factor" variables for the terminal event.
Xlevels	Vector of the values that factor might have taken for the recurrent events.
Xlevels2	Vector of the values that factor might have taken for the terminal event.
linear.pred	Linear predictor for the recurrent part: $\beta$'$X$ + log($u$i).
lineardeath.pred	Linear predictor for the terminal part: $\beta$'$X$ + $\alpha$ x log($u$i).
Xlevels	Vector of the values that factor might have taken for the recurrent part.
Xlevels2	vector of the values that factor might have taken 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$.

Note

In the flexible semiparametric case, smoothing parameters kappa and number of knots 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

J. Chauvet and V. Rondeau (2021). A flexible class of generalized joint frailty models for the analysis of survival endpoints. In revision.

Liu XR, Pawitan Y, Clements M. (2018) Parametric and penalized generalized survival models. Statistical Methods in Medical Research 27(5), 1531-1546.

Liu XR, Pawitan Y, Clements MS. (2017) Generalized survival models for correlated time-to-event data. Statistics in Medicine 36(29), 4743-4762.

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.

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, 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

Surv, terminal, timedep

Examples

#############################################################################
# -----        GENERALIZED SURVIVAL MODELS (without frailties)       -----  #
#############################################################################

adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50

### ---  Parametric PH, AH, PO and probit models  --- ###

GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino, 
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino, 
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino, 
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino, 
hazard="parametric", family="probit")

### ---  Semi-parametric PH and AH models  --- ###

GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino, 
family="PH", hazard="Splines", n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino, 
family="AH", hazard="Splines", n.knots=8, kappa=10^10, betaknots=1, betaorder=2)



#############################################################################
# -----          SHARED FRAILTY GENERALIZED SURVIVAL MODELS          -----  #
#############################################################################

adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50

### ---  Parametric PH, AH, PO and probit models  --- ###

GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino, 
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino, 
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino, 
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino, 
hazard="parametric", family="probit")

### ---  Semi-parametric PH and AH models  --- ###

GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt), 
data=adult.retino, family="PH", hazard="Splines", 
n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt), 
data=adult.retino, family="AH", hazard="Splines", 
n.knots=8, kappa=10^10, betaknots=1, betaorder=2)



#############################################################################
# -----         JOINT FRAILTY GENERALIZED SURVIVAL MODELS            -----  #
#############################################################################

data("readmission") 
readmission[, 3:5] = readmission[, 3:5]/365.25

### ---  Parametric dual-PH, AH, PO and probit models  --- ###

GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE, 
hazard="parametric", family=c("PH","PH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE, 
hazard="parametric", family=c("AH","AH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE, 
hazard="parametric", family=c("PO","PO"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE, 
hazard="parametric", family=c("probit","probit"))

### ---  Semi-parametric dual-PH and AH models  --- ###

GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE, 
hazard="Splines", family=c("PH","PH"), 
n.knots=5, kappa=c(100,100), betaknots=1, betaorder=3)
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE, 
hazard="Splines", family=c("AH","AH"), 
n.knots=5, kappa=c(600,600), betaknots=1, betaorder=3)

---