frailtyPenal: Fit a Shared, Joint or Nested Frailty model
In frailtypack: Shared, Joint (Generalized) Frailty Models; Surrogate Endpoints
frailtyPenal R Documentation
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 \omegai, of a
shared gamma frailty model for the jth subject in the ith
group:
where \lambda0(t) is the baseline hazard function, \beta
the vector of the regression coefficient associated to the covariate vector
Zij for the jth individual in the ith
group.
Otherwise, in case of a shared log-normal frailty model, we have for the
jth subject in the ith 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 \omegai 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 (\omegai for the
recurrent rate and \omegai\eqn{\alpha} for the death rate). The
covariates could be different for the recurrent rate and death rate.
For the jth recurrence (j=1,...,ni) and the
ith subject (i=1,...,G), the joint gamma frailty model
for recurrent event hazard function rij(.) and death rate
\lambdai is :
where r0(t) (resp. \lambda0(t)) is the recurrent (resp.
terminal) event baseline hazard function, \beta1 (resp.
\beta2) the regression coefficient vector, Zi(t)
the covariate vector. The random effects of frailties \omegai ~ \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 ui represents the unobserved association between recurrences and
death. The frailty term vi is specific to the recurrent event rate.
Thus, the general joint frailty model is:
where the iid random effects
ui ~ \Gamma(1/\theta,1/\theta) and the
iid random effects
vi ~ \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
Zi(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 vi and
wij for the kth individual of the jth subgroup of
the ith group is :
where \lambda0(t) is the baseline hazard function, Xijk
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 (ufi) and one for the group (wf) 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 (ufi, wf\eqn{\xi} for the recurrent rate and
ufi\eqn{\alpha}, wi for the terminal event rate). The covariates
could be different for the recurrent rate and death rate.
For the jth recurrence (j = 1, ..., ni) of the ith
individual (i = 1, ..., mf) of the fth group (f = 1,...,
n), the joint nested gamma frailty model for recurrent event hazard function
rfij(.) and for terminal event hazard function \lambdafi
is :
where r0(resp. \lambda0) is the recurrent (resp.
terminal) event baseline hazard function, \beta (resp. \gamma)
the regression coefficient vector, Xfij(t) the covariates
vector. The random effects are \omegaf ~ \Gamma(1/\eta,1/\eta)
and ufi ~ \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))
###--- 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")
frailtypack documentation built on Nov. 25, 2023, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| frailtyPenal | R Documentation |
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 \omegai, of a
shared gamma frailty model for the jth subject in the ith
group:
where \lambda0(t) is the baseline hazard function, \beta
the vector of the regression coefficient associated to the covariate vector
Zij for the jth individual in the ith
group.
Otherwise, in case of a shared log-normal frailty model, we have for the jth subject in the ith 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 \omegai 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 (\omegai for the
recurrent rate and \omegai\eqn{\alpha} for the death rate). The
covariates could be different for the recurrent rate and death rate.
For the jth recurrence (j=1,...,ni) and the
ith subject (i=1,...,G), the joint gamma frailty model
for recurrent event hazard function rij(.) and death rate
\lambdai is :
where r0(t) (resp. \lambda0(t)) is the recurrent (resp.
terminal) event baseline hazard function, \beta1 (resp.
\beta2) the regression coefficient vector, Zi(t)
the covariate vector. The random effects of frailties \omegai ~ \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 ui represents the unobserved association between recurrences and
death. The frailty term vi is specific to the recurrent event rate.
Thus, the general joint frailty model is:
where the iid random effects
ui ~ \Gamma(1/\theta,1/\theta) and the
iid random effects
vi ~ \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
Zi(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 vi and
wij for the kth individual of the jth subgroup of
the ith group is :
where \lambda0(t) is the baseline hazard function, Xijk
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 (ufi) and one for the group (wf) 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 (ufi, wf\eqn{\xi} for the recurrent rate and
ufi\eqn{\alpha}, wi for the terminal event rate). The covariates
could be different for the recurrent rate and death rate.
For the jth recurrence (j = 1, ..., ni) of the ith
individual (i = 1, ..., mf) of the fth group (f = 1,...,
n), the joint nested gamma frailty model for recurrent event hazard function
rfij(.) and for terminal event hazard function \lambdafi
is :
where r0(resp. \lambda0) is the recurrent (resp.
terminal) event baseline hazard function, \beta (resp. \gamma)
the regression coefficient vector, Xfij(t) the covariates
vector. The random effects are \omegaf ~ \Gamma(1/\eta,1/\eta)
and ufi ~ \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
|
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 |
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 |
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 |
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 |
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 |
init.B |
A vector of initial values for regression coefficients. This vector should be of the same size as the whole vector of covariates with the first elements for the covariates related to the recurrent events and then to the terminal event (interactions in the end of each component). Default is 0.1 for each (for Cox and shared model) or 0.5 (for joint and joint nested frailty models). |
init.Theta |
Initial value for variance of the frailties. |
init.Alpha |
Only for joint and joint nested frailty models : initial value for parameter alpha. |
Alpha |
Only for joint and joint nested frailty model : input "None" so as to fit a joint model without the parameter alpha. |
init.Ksi |
Only for joint nested frailty model : initial value for
parameter |
Ksi |
Only for joint nested frailty model : input |
init.Eta |
Only for general joint and joint nested frailty models :
initial value for the variance |
LIMparam |
Convergence threshold of the Marquardt algorithm for the
parameters (see Details), |
LIMlogl |
Convergence threshold of the Marquardt algorithm for the
log-likelihood (see Details), |
LIMderiv |
Convergence threshold of the Marquardt algorithm for the
gradient (see Details), |
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).
|
AIC |
the Akaike information Criterion for the parametric case.
|
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:
|
intcens |
Logical value. Indicator if a joint frailty model with interval-censored data was fitted) |
theta |
variance of the
gamma frailty parameter |
sigma2 |
variance
of the log-normal frailty parameter |
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 |
sigma2 |
variance of the log-normal frailty parameter
|
eta |
variance
of the second gamma frailty parameter in general joint frailty models
|
indic_alpha |
indicator if a joint frailty
model with |
alpha |
the coefficient
|
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 |
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 |
eta |
variance of the subcluster effect |
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 |
eta |
variance of the cluster effect |
alpha |
the power coefficient |
ksi |
the power coefficient |
indic_alpha |
indicator if a joint frailty model with |
indic_ksi |
indicator if a joint frailty
model with |
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 |
ksi_p.value |
p-value of the Wald test for the estimated power
coefficient |
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))
###--- 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")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.