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^{\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 :
$$\left{
\begin{array}{ll}
r_{ij}(t|\omega_i)=\omega_ir_0(t)\exp(\bold{\beta_1^{'}Z_i(t)}) & \mbox{(Recurrent)} \
\lambda_i(t|\omega_i)=\omega_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_i(t)}) & \mbox{(Death)} \
\end{array}
\right.$$
where $r_0(t)$ (resp. $\lambda_0(t)$) is the recurrent (resp. terminal) event baseline hazard function, $\bold{\beta_1}$ (resp. $\bold{\beta_2}$) the regression coefficient vector, $\bold{Z_i(t)}$ the covariate vector. The random effects of frailties $\omega_i\sim\bold{\Gamma}(\frac{1}{\theta},\frac{1}{\theta})$ and are iid.
The joint log-normal frailty model will be :
$$\left{
\begin{array}{ll}
r_{ij}(t|\eta_i)=r_0(t)\exp(\eta_i+\bold{\beta_1^{'}Z_i(t)}) & \mbox{(Recurrent)} \
\lambda_i(t|\eta_i)=\lambda_0(t)\exp(\alpha \eta_i+\bold{\beta_2^{'}Z_i(t)}) & \mbox{(Death)} \
\end{array}
\right.$$
where $$\eta_i\sim N(0,\sigma^2)$$
- 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 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.
$$\left{
\begin{array}{ll}
r_{ij}(t|u_i)=u_ir_0(t)\exp(\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to event)} \
\lambda_{ij}(t|u_i)=u_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_{ij}(t)}) & \mbox{(Death)} \
\end{array}
\right.$$
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 :
$$\left{
\begin{array}{ll}
r_{ij}(t|v_i)=r_0(t)\exp(v_i+\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to event)} \
\lambda_{ij}(t|v_i)=\lambda_0(t)\exp(\alpha v_i+\bold{\beta_2^{'}Z_{ij}(t)}) & \mbox{(Death)} \
\end{array}
\right.$$
where $$v_i\sim N(0,\sigma^2)$$
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.
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:
$\left{
\begin{array}{ll}
r_{ij}(t|u_i,v_i)=\textcolor{red}{u_iv_i}r_0(t)\exp(\bold{\beta_1^{'}Z_{ij}(t)}) =\textcolor{red}{u_iv_i}r_{ij}(t) & \mbox{(Reccurent)} \
\lambda_{i}(t|u_i)=\textcolor{red}{u_i}\lambda_0(t)\exp(\bold{\beta_1^{'}Z_{i}(t)}) = \textcolor{red}{u_i}\lambda_{i}(t) & \mbox{(Death)} \
\end{array}
\right.$
where the $iid$ random effects $\bold{u_i}\sim\Gamma(\frac{1}{\theta},\frac{1}{\theta})$ and the $iid$ random effects $\bold{v_i}\sim\Gamma(\frac{1}{\eta},\frac{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 $\bold{Z}_i(t)$ are allowed.
Nested Frailty model
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 :
$$\left{
\begin{array}{ll}
\lambda_{ijk}(t|v_i,w_{ij})=v_iw_{ij}\lambda_0(t)exp(\bold{\beta^{'}X_{ijk}}) \
v_i\sim\Gamma\left(\frac{1}{\alpha},\frac{1}{\alpha}\right) \hspace{0.05cm}i.i.d. \hspace{0.2cm} \bold{E}(v_i)=1 \hspace{0.2cm}\bold{Var}(v_i)=\alpha \
w_{ij}\sim\Gamma\left(\frac{1}{\eta},\frac{1}{\eta}\right)\hspace{0.05cm}i.i.d. \hspace{0.2cm} \bold{E}(w_{ij})=1 \hspace{0.2cm} \bold{Var}(w_{ij})=\eta
\end{array}
\right.$$
where $\lambda_0(t)$ is the baseline hazard function, $X_{ijk}$
denotes the covariate vector and $\beta$ the corresponding vector of regression parameters.frailtyPenal(formula, formula.terminalEvent, data,
recurrentAG = FALSE, cross.validation = FALSE,
jointGeneral,n.knots, kappa, maxit = 350,
hazard = "Splines", nb.int, RandDist = "Gamma",
betaknots = 1, betaorder = 3, init.B, init.Theta,
init.Alpha, Alpha, init.Eta, LIMparam = 1e-3, LIMlogl = 1e-3,
LIMderiv = 1e-3, print.times = TRUE, ...)jointGeneral = TRUE , the log-normal distribution cannot be chosen.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.SurvIC, subcluster, terminal, num.id, timedep###--- 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 clustered data ---###
#-- here is generated cluster (5 clusters)
readmission <- transform(readmission,group=id
#-- 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 ---###
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))Run the code above in your browser using DataLab