GenfrailtyPenal: Fit a Shared or a Joint Frailty Generalized Survival Model

Description

I. SHARED FRAILTY GENERALIZED SURVIVAL MODELSFit a gamma Shared Frailty Generalized Survival Model using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) are allowed. It allows to obtain a parametric or flexible semi-parametric smooth hazard and survival functions.Each frailty term $u$_i is assumed to act multiplicatively on the hazard function, and to be drawn from a Gamma distribution with unit mean and variance $\theta$. Conditional on the frailty term, the hazard function for the $j$^th subject in the $i$^th group is then expressed bywhere $x$_ij is a collection of baseline covariates, $\xi$ is a vector of parameters, and $\lambda$_ij ($t$ | $x$_ij ; $\xi$) is the hazard function for an average value of the frailty. The associated conditional survival function writeswhere $S$_ij ($t$ | $x$_ij ; $\xi$) designates the survival function for an average value of the frailty.Following Liu et al. (2017, 2018), the latter function is expressed in terms of a link function $g$(.) and a linear predictor $\eta$_ij ($t$, $x$_ij; $\xi$) such that $g$[$S$_ij ($t$ | $x$_ij ; $\xi$)] = $\eta$_ij ($t$, $x$_ij; $\xi$), i.e. $S$_ij ($t$ | $x$_ij ; $\xi$) = $h$[$\eta$_ij ($t$, $x$_ij; $\xi$)] with $h$() = $g$^-1().The conditional survival function is finally modeled byThe table below summarizes the most commonly used (inverse) link functions and their associated conditional survival, hazard and cumulative hazard functions. PHM stands for "Proportional Hazards Model", POM for "Proportional Odds Model, PROM for "Probit Model" and AHM for "Additive Hazards Model".I.(a) Fully parametric caseIn the fully parametric case, linear predictors considered are of the formwhere $\rho > 0$ is a shape parameter, $\gamma > 0$ a scale parameter, $\beta$ a vector of regression coefficients, and $\xi$ = ($\rho$ ,$\gamma$, $\beta$). With the appropriate link function, such linear parametric predictors make it possible to recover a Weibull baseline survival function for PHMs and AHMs, a log-logistic baseline survival function for POMs, and a log-normal one for PROMs. I. (b) Flexible semi-parametric caseFor PHM and AHM, a more flexible splines-based approach is proposed for modeling the baseline hazard function and time-varying regression coefficients. In this case, conditional on the frailty term $u$_i, the hazard function for the $j$^th subject in the $i$^th group is still expressed by $\lambda$_ij ($t$ | $x$_ij, $u$_i ; $\xi$) = $u$_i $\lambda$_ij ($t$ | $x$_ij ; $\xi$), but we have this time The smoothness of baseline hazard function $\lambda$₀() is ensured by penalizing the log-likelihood by a term which has large values for rough functions. Moreover, for parametric and flexible semi-parametric AHMs, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.II. JOINT FRAILTY GENERALIZED SURVIVAL MODELSFit a gamma Joint Frailty Generalized Survival Model for recurrent and terminal events using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) for the recurrent event part are allowed. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring.This model includes an common patient-specific frailty term $u$_i for the two survival functions which will take into account the unmeasured heterogeneity in the data, associated with unobserved covariates. The frailty term acts differently for the two survival functions ($u$_i for the recurrent survival function and $u$_i^α for the death one). The covariates could be different for the recurrent and terminal event parts.II.(a) Fully parametric caseFor the $j$^th recurrence (j=1,...,n_i) and the $i$^th patient (i=1,...,N), the gamma Joint Frailty Generalized Survival Model for recurrent event survival function $S$_Rij(.) and death survival function $S$_Di(.) is- $\eta$_Rij (resp. $\eta$_Di) is the linear predictor for the recurrent (resp. terminal) event process. The form of these linear predictors is the same as the one presented in I.(a).- $h$_R(.) (resp. $h$_D(.)) is the inverse link function associated with recurrent events (resp. terminal event).- $x$_Rij and $x$_Di are two vectors of baseline covariates associated with recurrent and terminal events.- $\xi$_R and $\xi$_D are the parameter vectors for recurrent and terminal events.- $\alpha$ is a parameter allowing more flexibility in the association between recurrent and terminal events processes.- The random frailties $u$_i are still assumed iid and drown from a $\Gamma$(1/$\theta$,1/$\theta$).II.(b) Flexible semi-parametric caseIf one chooses to fit a PHM or an AHM for recurrent and/or terminal events, a splines-based approach for modeling baseline hazard functions and time-varying regression coefficients is still available. In this approach, the submodel for recurrent events is expressed as $\lambda$_Rij ($t$ | $x$_Rij, $u$_i ; $\xi$_R) = $u$_i $\lambda$_Rij ($t$ | $x$_Rij ; $\xi$_R), where The submodel for terminal event is expressed as $\lambda$_Di ($t$ | $x$_Di, $u$_i ; $\xi$_D) = $u$_i^α $\lambda$_Di ($t$ | $x$_Di ; $\xi$_D), where Baseline hazard functions $\lambda$_R0(.) and $\lambda$_D0(.) are estimated using cubic M-splines (of order 4) with positive coefficients, and the time-varying coefficients $\beta$_R(.) and $\beta$_D(.) are estimated using B-splines of order q.The smoothness of baseline hazard functions is ensured by penalizing the log-likelihood by two terms which has large values for rough functions. Moreover, if one chooses an AHM for recurrent and/or terminal event submodel, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.

Usage

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

Value

The following components are included in a 'frailtyPenal' object for each model.

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

The following components are specific to shared models.

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

The following components are specific to joint models.

formula: The formula part of the code used for the recurrent events.
formula.terminalEvent: The formula part of the code used for the terminal model.
n.deaths: Number of observed deaths.
n.censored: Number of censored individuals.
theta: Variance of the gamma frailty parameter, i.e. Var($u$_i).
indic_alpha: Indicator if a joint frailty model with $\alpha$ parameter was fitted.
alpha: The coefficient $\alpha$ associated with the frailty parameter in the terminal hazard function.
nvar: A vector with the number of covariates of each type of hazard function as components.
nvarnotdep: A vector with the number of constant effect covariates of each type of hazard function as components.
nvarRec: Number of recurrent explanatory variables.
nvarEnd: Number of death explanatory variables.
noVar1: Indicator of recurrent explanatory variables.
noVar2: Indicator of death explanatory variables.
Names.vardep: Name of the covariates with time-varying effects for the recurrent events.
Names.vardepdc: Name of the covariates with time-varying effects for the terminal event.
xR: Matrix of times where both survival and hazard function are estimated for the recurrent event.
xD: Matrix of times for the terminal event.
lamR: Array (dim=3) of hazard estimates and confidence bands for recurrent event.
lamD: The same value as lamR for the terminal event.
survR: Array (dim=3) of baseline survival estimates and confidence bands for recurrent event.
survD: The same value as survR for the terminal event.
nb.gh: Number of nodes for the Gaussian-Hermite quadrature.
nb.gl: Number of nodes for the Gaussian-Laguerre quadrature.
medianR: The value of the median survival for the recurrent events and its confidence bands.
medianD: The value of the median survival for the terminal event and its confidence bands.
names.factor: Names of the "as.factor" variables for the recurrent events.
names.factordc: Names of the "as.factor" variables for the terminal event.
Xlevels: Vector of the values that factor might have taken for the recurrent events.
Xlevels2: Vector of the values that factor might have taken for the terminal event.
linear.pred: Linear predictor for the recurrent part: $\beta$'$X$ + log($u$_i).
lineardeath.pred: Linear predictor for the terminal part: $\beta$'$X$ + $\alpha$ x log($u$_i).
Xlevels: Vector of the values that factor might have taken for the recurrent part.
Xlevels2: vector of the values that factor might have taken for the death part.
BetaTpsMat: Matrix of time varying-effects and confidence bands for recurrent event (the first column used for abscissa of times of recurrence).
BetaTpsMatDc: Matrix of time varying-effects and confidence bands for terminal event (the first column used for abscissa of times of death).
alpha_p.value: p-value of the Wald test for the estimated $\alpha$.

Arguments

formula: A formula object, with the response on the left of a $\sim$ operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function like in survival package. Interactions are possible using ' * ' or ' : '.
formula.terminalEvent: Only for joint frailty models: a formula object, only requires terms on the right to indicate which variables are used for the terminal event. Interactions are possible using ' * ' or ' : '.
data: A 'data.frame' with the variables used in 'formula'.
recurrentAG: Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent covariates. The default is FALSE.
family: Type of Generalized Survival Model to fit. "PH" for a proportional hazards model, "AH" for an additive hazards model, "PO" for a proportional odds model and "probit" for a probit model. A vector of length 2 is expected for joint models (e.g., family=c("PH","PH")).
hazard: Type of hazard functions: "Splines" for semi-parametric hazard functions using equidistant intervals, or "parametric" for parametric distribution functions. In case of family="PH" or family="AH", the "parametric" option corresponds to a Weibull distribution. In case of family="PO" and family="probit", the "parametric" option corresponds to a log-logistic and a log-normal distribution, respectively. So far, the "Splines" option is only available for PH and AH submodels. Default is "Splines".
n.knots: Integer giving the number of knots to use. Value required in the penalized likelihood estimation. It corresponds to the n.knots+2 splines functions for the approximation of the hazard or the survival functions. We estimate I- or M-splines of order 4. When the user set a number of knots equals to k (i.e. n.knots=k), then the number of interior knots is k-2 and the number of splines is (k-2)+order. Number of knots must be between 4 and 20. (See Note)
kappa: Positive smoothing parameter in the penalized likelihood estimation. The coefficient kappa tunes the intensity of the penalization (the integral of the squared second derivative of hazard function). In a stratified shared model, this argument must be a vector with kappas for both strata. In a stratified joint model, this argument must be a vector with kappas for both strata for recurrent events plus one kappa for terminal event. We advise the user to identify several possible tuning parameters, note their defaults and look at the sensitivity of the results to varying them. Value required. (See Note).
betaknots: Number of inner knots used for the B-splines time-varying coefficient estimation. Default is 1. See 'timedep' function for more details.
betaorder: Order of the B-splines used for the time-varying coefficient estimation. Default is cubic B-splines (order=3). See 'timedep' function for more details. Not implemented for Proportional Odds and Probit submodels.
RandDist: Type of random effect distribution: "Gamma" for a gamma distribution, and "LogN" for a log-normal distribution (not implemented yet). Default is "Gamma".
init.B: A vector of initial values for regression coefficients. This vector should be of the same size as the whole vector of covariates with the first elements for the covariates related to the recurrent events and then to the terminal event (interactions in the end of each component). Default is 0.1 for each (for Generalized Survival and Shared Frailty Models) or 0.5 (for Generalized Joint Frailty Models).
init.Theta: Initial value for frailty variance.
init.Alpha: Only for Generalized Joint Frailty Models: initial value for parameter alpha.
Alpha: Only for Generalized Joint Frailty Models: input "None" so as to fit a joint model without the parameter alpha.
maxit: Maximum number of iterations for the Marquardt algorithm. Default is 300
nb.gh: Number of nodes for the Gaussian-Hermite quadrature. It can be chosen among 5, 7, 9, 12, 15, 20 and 32. The default is 20 if hazard="Splines", 32 otherwise.
nb.gl: Number of nodes for the Gaussian-Laguerre quadrature. It can be chosen between 20 and 32. The default is 20 if hazard="Splines", 32 otherwise.
LIMparam: Convergence threshold of the Marquardt algorithm for the parameters (see Details), $10$^-3 by default.
LIMlogl: Convergence threshold of the Marquardt algorithm for the log-likelihood (see Details), $10$^-3 by default.
LIMderiv: Convergence threshold of the Marquardt algorithm for the gradient (see Details), $10$^-3 by default.
print.times: A logical parameter to print iteration process. Default is TRUE.
cross.validation: Not implemented yet for the generalized settings.
jointGeneral: Not implemented yet for the generalized settings.
nb.int: Not implemented yet for the generalized settings.
initialize: Not implemented yet for the generalized settings.
init.Ksi: Not implemented yet for the generalized settings.
Ksi: Not implemented yet for the generalized settings.
init.Eta: Not implemented yet for the generalized settings.

Details

TYPICAL USESFor a Generalized Survival Model:

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

For a Shared Frailty Generalized Survival Model:

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

For a Joint Frailty Generalized Survival Model:

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

OPTIMIZATION ALGORITHMThe estimated parameters are obtained using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm. The iterations are stopped when the difference between two consecutive log-likelihoods is small ($<10$^-3), the estimated coefficients are stable (consecutive values $<10$^-3, and the gradient small enough ($<10$^-3). When the frailty variance is small, numerical problems may arise. To solve this problem, an alternative formula of the penalized log-likelihood is used (see Rondeau, 2003 for further details). For Proportional Hazards and Additive Hazards submodels, cubic M-splines of order 4 can be used to estimate the hazard function. In this case, I-splines (integrated M-splines) are used to compute the cumulative hazard function.The inverse of the Hessian matrix is the variance estimator. To deal with the positivity constraint of the variance component and the spline coefficients, a squared transformation is used and the standard errors are computed by the $\Delta$-method (Knight & Xekalaki, 2000). The integrations in the full log likelihood are evaluated using Gaussian quadrature. Laguerre polynomials with 20 points are used to treat the integrations on $[0,\infty[$.INITIAL VALUESIn case of a shared frailty model, the splines and the regression coefficients are initialized to 0.1. The program fits, firstly, an adjusted Cox model to give new initial values for the splines and the regression coefficients. The variance of the frailty term $\theta$ is initialized to 0.1. Then, a shared frailty model is fitted.In case of a joint frailty model, the splines and the regression coefficients are initialized to 0.5. The program fits firstly, an adjusted Cox model to have new initial values for the splines and the regression coefficients. The variance of the frailty term $\theta$ and the association parameter $\alpha$ are initialized to 1. Then, a joint frailty model is fitted.

References

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

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

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

A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels and V. Rondeau (2017). Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event. Journal of Statistical Software 81(3), 1-52.

V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametric estimation. Journal of Statistical Software 47, 1-28.

V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the dependance between clustered times to tumour progression and deaths: A meta-analysis of chemotherapy in head and neck cancer. Statistical methods in medical research 897, 1-19.

V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P. Soubeyran (2007). Joint frailty models for recurring events and death using maximum penalized likelihood estimation:application on cancer events. Biostatistics 8,4, 708-721.

V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.

C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in survival analysis. Biometrics 47, 461-466.

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

Examples

Run this code

# \donttest{

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

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

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

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

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

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



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

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

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

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

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

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



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

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

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

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

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

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

# }

Run the code above in your browser using DataLab