frailtySurv is an R package for simulating and fitting semiparametricshared frailty models.
For the latest stable version, install from CRAN:
install.packages("frailtySurv")The development version can be installed from github usingdevtools.
devtools::install_github("vmonaco/frailtySurv")The following code shows how to generate data and fit the model.
set.seed(1234)library(frailtySurv)#> Loading required package: survival#> Welcome to frailtySurv v1.3.3dat<-genfrail(N=200,K=2,beta=c(log(2),log(3)),frailty="gamma",theta=2,censor.rate=0.35,Lambda_0=function(t,tau=4.6,C=0.01) (C*t)^tau)fit<-fitfrail(Surv(time, status)~ Z1+ Z2+cluster(family),dat=dat,frailty="gamma")fit#> Call: fitfrail(formula = Surv(time, status) ~ Z1 + Z2 + cluster(family),#> dat = dat, frailty = "gamma")#>#> Covariate Coefficient#> Z1 0.654#> Z2 1.006#>#> Frailty distribution gamma(1.802), VAR of frailty variates = 1.802#> Log-likelihood -1575.752#> Converged (method) 10 iterations, 3.24 secs (maximized log-likelihood)Parameter traces are given by
plot(fit,type="trace")
The estimated cumulative baseline hazard is given by
plot(fit,type="cumhaz")
The results can be compared to other estimation techniques.
library(survival)library(frailtypack)coxph(Surv(time, status)~ Z1+ Z2+frailty.gamma(family),data=dat)#> Call:#> coxph(formula = Surv(time, status) ~ Z1 + Z2 + frailty.gamma(family),#> data = dat)#>#> coef se(coef) se2 Chisq DF p#> Z1 0.6651 0.0941 0.0852 49.9190 1 1.6e-12#> Z2 1.0282 0.1046 0.0941 96.5685 1 < 2e-16#> frailty.gamma(family) 475.6283 150 < 2e-16#>#> Iterations: 7 outer, 58 Newton-Raphson#> Variance of random effect= 1.88 I-likelihood = -1319.6#> Degrees of freedom for terms= 0.8 0.8 149.7#> Likelihood ratio test=576 on 151 df, p=<2e-16 n= 400frailtyPenal(Surv(time, status)~ Z1+ Z2+cluster(family),data=dat,n.knots=10,kappa=2)#>#> Be patient. The program is computing ...#> The program took 0.22 seconds#> Call:#> frailtyPenal(formula = Surv(time, status) ~ Z1 + Z2 + cluster(family),#> data = dat, n.knots = 10, kappa = 2)#>#>#> Shared Gamma Frailty model parameter estimates#> using a Penalized Likelihood on the hazard function#>#> coef exp(coef) SE coef (H) SE coef (HIH) z p#> Z1 0.67236 1.95885 0.0994732 0.0994732 6.75921 1.3875e-11#> Z2 1.04884 2.85433 0.1191210 0.1191210 8.80481 0.0000e+00#>#> Frailty parameter, Theta: 1.92378 (SE (H): 0.319435 ) p = 8.5895e-10#>#> penalized marginal log-likelihood = -1378.15#> Convergence criteria:#> parameters = 0.000362 likelihood = 0.000201 gradient = 2.44e-08#>#> LCV = the approximate likelihood cross-validation criterion#> in the semi parametrical case = 3.48287#>#> n= 400#> n events= 256 n groups= 200#> number of iterations: 13#>#> Exact number of knots used: 10#> Value of the smoothing parameter: 2, DoF: 11.96You can clone the repository and build the project from source usingRStudio. To create a project in RStudio from this repository, you musthave both RStudio and git installed.
https://github.com/vmonaco/frailtySurvSee theR packages book
To load the package (simulate installing and loading withlibrary(frailty)):
devtools::load_all()or shortuct:Ctrl+Shift+L
Use this command to build a source package (binary = FALSE bydefault)
devtools::build()Other usefull commands and shortcuts:
devtools::check():Ctrl+Shift+Edevtools::test():Ctrl+Shift+TBuild the vignettes with the command:
devtools::build_vignettes()Compiled vignettes will reside ininst/doc. Inparticular, see frailtySurv.pdf
Use RStudio to pull the latest code if you’ve already cloned therepository. In the “git” pane, look for the Pull button. You can alsoview the commit history from there and see the diff between modifiedfiles, if any.