Learn R Programming

frailtypack (version 2.7.1)

frailtyPenal: 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, conditionnal on the frailty term $\omega_i$, of a shared gamma frailty model for the $j^{th}$ subject in the $i^{th}$ group: $$\lambda_{ij}(t|\omega_i)=\lambda_0(t)\omega_i\exp(\bold{\beta^{'}Z_{ij}})$$ $$\omega_i\sim\Gamma\left(\frac{1}{\theta},\frac{1}{\theta}\right) \hspace{0.5cm} \bold{E}(\omega_i)=1 \hspace{0.5cm}\bold{Var}(\omega_i)=\theta$$ where $\lambda_0(t)$ is the baseline hazard function, $\bold{\beta}$ the vector of the regression coefficient associated to the covariate vector $\bold{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: $$\lambda_{ij}(t|\eta_i)=\lambda_0(t)\exp(\eta_i+\bold{\beta^{'}Z_{ij}})$$ $$\eta_i\sim N(0,\sigma^2)$$ 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^{\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. 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.$$ 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. 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.

Usage

frailtyPenal(formula, formula.terminalEvent, data,
             recurrentAG = FALSE, cross.validation = FALSE,
             n.knots, kappa, maxit = 350, hazard = "Splines",
             nb.int, RandDist = "Gamma", betaknots = 1,
             betaorder = 3, init.B, init.Theta, init.Alpha,
             Alpha, 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 $\texttildelow$ 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 inte
formula.terminalEvent
only for joint model: a formula object, only requires terms on the right to indicate which variables are modelling the terminal event.
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 cova
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 se
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 estim
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 f
maxit
maximum number of iterations for the Marquardt algorithm. Default is 350
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 perce
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 o
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.
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 model.
betaorder
Order of the B-splines. Defaut is cubic B-splines (order = 3). See 'timedep' function for more details. Not implemented for nested model.
init.B
A vector of inital values for regression coefficients. This vector should be of the same size as the whole vector of covariates. Default is 0.1 for each (for Cox and shared model) or 0.5 (for joint model).
init.Theta
Inital value for variance of the frailties.
init.Alpha
Only for joint model : initial value for parameter alpha.
Alpha
Only for joint model : input "none" so as to fit a joint model without the parameter alpha.
LIMparam
Convergence threshold of the Marquard algorithm for the parameters (see Details), $10^{-3}$ by default.
LIMlogl
Convergence threshold of the Marquard algorithm for the log-likelihood (see Details), $10^{-3}$ by default.
LIMderiv
Convergence threshold of the Marquard algorithm for the gradient (see Details), $10^{-3}$ by default.
print.times
a logical parameter to print iteration process. Default is TRUE.

Value

  • The following components are included in a 'frailtyPenal' object for each model.
  • bsequence of the corresponding estimation of the coefficients for the hazard functions (parametric or semiparametric), the random effects variances and the regression coefficients.
  • callThe code used for the model.
  • formulathe formula part of the code used for the model.
  • coefthe regression coefficients.
  • cross.ValLogical value. Is cross validation procedure used for estimating the smoothing parameters in the penalized likelihood estimation?
  • DoFDegrees of freedom associated with the "kappa".
  • formulathe formula part of the code used for the model.
  • groupsthe maximum number of groups used in the fit.
  • kappaA vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components.
  • loglikPenalthe complete marginal penalized log-likelihood in the semiparametric case.
  • loglikthe marginal log-likelihood in the parametric case.
  • nthe number of observations used in the fit.
  • n.eventsthe number of events observed in the fit.
  • n.iternumber of iterations needed to converge.
  • n.knotsnumber of knots for estimating the baseline functions in the penalized likelihood estimation.
  • n.stratnumber of stratum.
  • varHthe variance matrix of all parameters before positivity constraint transformation. Thenafter, 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.
  • varHIHthe robust estimation of the variance matrix of all parameters.
  • xmatrix 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.
  • lamarray (dim=3) of hazard estimates and confidence bands.
  • survarray (dim=3) of baseline survival estimates and confidence bands.
  • type.of.hazardType of hazard functions (0:"Splines", "1:Piecewise", "2:Weibull").
  • type.of.PiecewiseType of Piecewise hazard functions (1:"percentile", 0:"equidistant").
  • nbintervRNumber of intervals (between 1 and 20) for the parametric hazard functions ("Piecewise-per", "Piecewise-equi").
  • nparnumber of parameters.
  • nvarnumber of explanatory variables.
  • noVarindicator of explanatory variables.
  • LCVthe approximated likelihood cross-validation criterion in the semiparametric case (with H minus the converged hessien matrix, and l(.) the full log-likelihood).$$LCV=\frac{1}{n}(trace(H^{-1}_{pl}H) - l(.))$$
  • AICthe Akaike information Criterion for the parametric case.$$AIC=\frac{1}{n}(np - l(.))$$
  • n.knots.tempinitial value for the number of knots.
  • shape.weibshape parameter for the Weibull hazard function.
  • scale.weibscale parameter for the Weibull hazard function.
  • martingale.resmartingale residuals for each cluster.
  • martingaleCoxmartingale residuals for observation in the Cox model.
  • FrailtyLogical value. Was model with frailties fitted ?
  • frailty.predempirical Bayes prediction of the frailty term (ie, using conditional posterior distributions).
  • frailty.varvariance of the empirical Bayes prediction of the frailty term (only for gamma frailty models).
  • frailty.sdstandard error of the frailty empirical Bayes prediction (only for gamma frailty models).
  • global_chisqa vector with the values of each multivariate Wald test.
  • dof_chisqa vector with the degree of freedom for each multivariate Wald test.
  • global_chisq.testa binary variable equals to 0 when no multivariate Wald is given, 1 otherwise.
  • p.global_chisqa vector with the p_values for each global multivariate Wald test.
  • names.factorNames of the "as.factor" variables.
  • Xlevelsvector of the values that factor might have taken.
  • contraststype of contrast for factor variable.
  • The following components are specific to shared models.
  • thetavariance of the gamma frailty parameter $(\bold{Var}(\omega_i))$
  • sigma2variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$
  • linear.predlinear 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.
  • BetaTpsMatmatrix of time varying-effects and confidence bands (the first column used for absciss of times)
  • The following components are specific to joint models.
  • thetavariance of the gamma frailty parameter $(\bold{Var}(\omega_i))$ or $(\bold{Var}(u_i))$
  • sigma2variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$ or $(\bold{Var}(v_i))$
  • indic_alphaindicator if a joint frailty model with $\alpha$ parameter was fitted
  • alphathe coefficient $\alpha$ associated with the frailty parameter in the terminal hazard function.
  • nbintervRNumber of intervals (between 1 and 20) for the recurrent parametric hazard functions ("Piecewise-per", "Piecewise-equi").
  • nbintervDCNumber of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewise-per", "Piecewise-equi").
  • nvarA vector with the number of covariates of each type of hazard function as components.
  • nvarRecnumber of recurrent explanatory variables.
  • nvarEndnumber of death explanatory variables.
  • noVar1indicator of recurrent explanatory variables.
  • noVar2indicator of death explanatory variables.
  • xRmatrix 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.
  • xDmatrix of times for the terminal event.
  • lamRarray (dim=3) of hazard estimates and confidence bands for recurrent event.
  • lamDthe same value as lamR for the terminal event.
  • survRarray (dim=3) of baseline survival estimates and confidence bands for recurrent event.
  • survDthe same value as survR for the terminal event.
  • martingale.resmartingale residuals for each cluster (recurrent).
  • martingaledeath.resmartingale residuals for each cluster (death).
  • linear.predlinear predictor: uses "Beta'X + log w_i" in the gamma frailty model, otherwise uses "Beta'X + eta_i" for log-normal frailty distribution
  • lineardeath.predlinear predictor for the terminal part : "Beta'X + alpha.log w_i" for gamma, "Beta'X + alpha.eta_i" for log-normal frailty distribution
  • Xlevelsvector of the values that factor might have taken for the recurrent part.
  • contraststype of contrast for factor variable for the recurrent part.
  • Xlevels2vector of the values that factor might have taken for the death part.
  • contrasts2type of contrast for factor variable for the death part.
  • BetaTpsMatmatrix of time varying-effects and confidence bands for recurrent event (the first column used for absciss of times of recurrence)
  • BetaTpsMatDcmatrix of time varying-effects and confidence bands for terminal event (the first column used for absciss of times of death)
  • The following components are specific to nested models.
  • alphavariance of the cluster effect $(\bold{Var}(v_{i}))$
  • etavariance of the subcluster effect $(\bold{Var}(w_{ij}))$
  • subgroupsthe maximum number of subgroups used in the fit.
  • frailty.pred.groupempirical Bayes prediction of the frailty term by group.
  • frailty.pred.subgroupempirical Bayes prediction of the frailty term by subgroup.
  • linear.predlinear predictor: uses "Beta'X + log v_i.w_ij".
  • subgbygsubgroup by group.
  • n.stratA vector with the number of covariates of each type of hazard function as components.

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 nested model frailtyPenal(Surv(time,event)~cluster(group)+subcluster(sbgroup)+ var1+var2, 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 loglikelihoods was small $(<10^{-3})$, 4="" 20="" 2003="" the="" estimated="" coefficients="" were="" stable="" (consecutive="" values="" $(<10^{-3})$,="" and="" gradient="" small="" enough="" $(<10^{-3})$.="" when="" frailty="" parameter="" is="" small,="" numerical="" problems="" may="" arise.="" to="" solve="" this="" problem,="" an="" alternative="" formula="" of="" penalized="" log-likelihood="" used="" (see="" rondeau,="" for="" further="" details).="" cubic="" m-splines="" order="" are="" hazard="" function,="" i-splines="" (integrated="" m-splines)="" cumulative="" function.="" inverse="" hessian="" matrix="" variance="" estimator="" deal="" with="" positivity="" constraint="" component="" spline="" coefficients,="" a="" squared="" transformation="" standard="" errors="" computed="" by="" $\delta$-method="" (knight="" &="" xekalaki,="" 2000).="" smooth="" can="" be="" chosen="" maximizing="" likelihood="" cross="" validation="" criterion="" (joly="" other,="" 1998).="" integrations="" in="" full="" log="" evaluated="" using="" gaussian="" quadrature.="" laguerre="" polynomials="" points="" treat="" 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 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.

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, 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 Medecine, 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. McGilchrist CA and Aisbett CW (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, subcluster, terminal, num.id, timedep

Examples

Run this code
###---  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 a generated cluster (31 clusters of 13 subjects)
readmission <- transform(readmission,group=id%%31+1)

joi.clus <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(group)+
num.id(id)+dukes+charlson+sex+chemo+terminal(death),
formula.terminalEvent=~dukes+charlson+sex+chemo,
data=readmission,recurrentAG=TRUE,n.knots=10,
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