frailtyPenal for Shared frailty models: Fit a Shared Frailty model using a semiparametric penalized likelihood estimation or parametric estimation

Description

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 (max=2) 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)$$

Usage

frailtyPenal(formula, formula.terminalEvent, data, Frailty = FALSE,
            joint = FALSE, recurrentAG = FALSE, cross.validation =
            FALSE, n.knots, kappa1, kappa2, maxit = 350,
            hazard = "Splines", nb.int1, nb.int2, RandDist = "Gamma")

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 int

formula.terminalEvent

Not required.

data

a 'data.frame' in which to interpret the variables named in the 'formula'.

Frailty

Logical value. Is model with frailties fitted? If so, variance of frailty parameter is estimated. If not, Cox proportional hazard model is estimated using Penalized Likelihood on the hazard function. The default is TRUE

joint

Not required.

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 cov

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 kappa1 as the

n.knots

integer giving the number of knots to use. Value required in the penalized likelihood estimation. It corresponds to the (n.knots+2) splines functions for the approximation of the hazard or the survival functions. Number

kappa1

positive smoothing parameter in the penalized likelihood estimation. 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

kappa2

positive smoothing parameter in the penalized likelihood estimation for the second stratum, when data are stratified. See kappa1.

maxit

maximum number of iterations for the Marquardt algorithm. Default is 350

hazard

Type of hazard functions: "Splines" for semiparametric hazard functions with the penalized likelihood estimation, "Piecewise-per" for piecewise constant hazard function using percentile (not available for interval-censored data), "Piecewise-

nb.int1

Number of intervals (between 1 and 20) for the parametric hazard functions ("Piecewise-per", "Piecewise-equi").

nb.int2

Argument used only for the Joint frailty model.

RandDist

Type of random effect distribution: "Gamma" for a gamma distribution, "LogN" for a log-normal distribution. Default is "Gamma".

Value

The following components are included in a 'frailtyPenal' object for shared frailty models.
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.
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.
lammatrix of hazard estimates and confidence bands.
lam2the same value as lam for the second stratum.
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.
survmatrix of baseline survival estimates and confidence bands.
surv2the same value as surv for the second stratum.
thetavariance of the gamma frailty parameter $(\bold{Var}(\omega_i))$
sigma2variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$
varHthe variance matrix of all parameters before positivity constraint transformation (theta, the regression coefficients and the spline coefficients). Thenafter, the delta method is needed to obtain the estimated variance parameters.
varHIHthe robust estimation of the variance matrix of all parameters (theta, the regression coefficients and the spline coefficients).
x1vector 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.
x2vector of times for the second stratum (see x1 value).
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).
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" in log-normal.
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.

Details

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. 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. INITIAL VALUES 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.

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

Examples

Run this code

###  Shared model  ###

data(kidney)
frailtyPenal(Surv(time,status)~cluster(id)+sex+age,
             n.knots=12,kappa1=10000,data=kidney,Frailty=TRUE)

###  COX proportional hazard model (SHARED without frailties) ###
###  estimated with penalized likelihood ###

frailtyPenal(Surv(time,status)~sex+age,
             n.knots=12,kappa1=10000,data=kidney,Frailty=FALSE)

###  Shared model with truncated data ###

# Here is created a hypothetical truncated data

kidney$tt0<-rep(0,nrow(kidney))
kidney$tt0[1:3]<-c(2,9,13)

# then, we fit the model
frailtyPenal(Surv(tt0,time,status)~sex+age+cluster(id),
             n.knots=12,kappa1=10000,data=kidney,Frailty=TRUE)

###  Shared model - stratified analysis ###

data(readmission)
frailtyPenal(Surv(time,event)~cluster(id)+dukes+strata(sex),
             n.knots=10,kappa1=10000,kappa2=10000,data=readmission,
             Frailty=TRUE)

###  Shared model - recurrentAG=TRUE ###

frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex+dukes+
             charlson,data=readmission,Frailty=TRUE,
             n.knots=6,kappa1=100000,recurrentAG=TRUE)

###  Shared model - cross.validation=TRUE ###

frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex+dukes+
             charlson,data=readmission, Frailty=TRUE,n.knots=6,
             kappa1=5000,recurrentAG=TRUE,cross.validation=TRUE)


### Shared model - log-normal distribution ###

frailtyPenal(Surv(t.start,t.stop,event)~sex+dukes+charlson+cluster(id),
             data=readmission, Frailty=TRUE,n.knots=6,kappa1=5000,
             recurrentAG=TRUE,RandDist="LogN")

Run the code above in your browser using DataLab