R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("frailtypack-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('frailtypack')
Loading required package: survival
Loading required package: splines
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "frailtyPenal"
> 
> ### * frailtyPenal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: frailtyPenal
> ### Title: Fit Shared Gamma Frailty model using penalized likelihood
> ###   estimation
> ### Aliases: frailtyPenal
> ### Keywords: survival
> 
> ### ** Examples
> 
> data(kidney)
> #Shared frailty model
> frailtyPenal(Surv(time,status)~sex+age+cluster(id),
+              n.knots=12,kappa1=1000,data=kidney)
Call:
frailtyPenal(formula = Surv(time, status) ~ sex + age + cluster(id), 
    data = kidney, n.knots = 12, kappa1 = 1000)


  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
sex -1.58858     0.204      0.4488        0.4243 -3.539 0.0004
age  0.00739     1.007      0.0118        0.0118  0.624 0.5300

    Frailty parameter, Theta: 0.431 (SE (H): 0.233 ) (SE (HIH): 0.228 ) 
 
    penalized marginal log-likelihood = -319.88
    n= 76
    n events= 58  n groups= 38
    number of iterations:  34
    Exact number of knots used:  12 
    Value of the smoothing parameter:  1000, DoF:  13.98
> 
> #model without frailties (e.g., Cox proportional hazards 
> #                          estimated via penalized likelihood)
> frailtyPenal(Surv(time,status)~sex+age+cluster(id),
+              n.knots=12,kappa1=1000,data=kidney,Frailty=FALSE)
Call:
frailtyPenal(formula = Surv(time, status) ~ sex + age + cluster(id), 
    data = kidney, Frailty = FALSE, n.knots = 12, kappa1 = 1000)


  Cox proportional hazards model parameter estimates  
  using a Penalized Likelihood on the hazard function 

        coef exp(coef) SE coef (H) SE coef (HIH)     z      p
sex -0.92873     0.395     0.30734       0.30698 -3.02 0.0025
age  0.00358     1.004     0.00967       0.00966  0.37 0.7100

 
    penalized marginal log-likelihood = -323.25
    n= 76
    n events= 58  n groups= 38
    number of iterations:  46
    Exact number of knots used:  12 
    Value of the smoothing parameter:  1000, DoF:  13.99
> 
> # truncated data
> 
> # first, we create 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=1000,data=kidney)
Call:
frailtyPenal(formula = Surv(tt0, time, status) ~ sex + age + 
    cluster(id), data = kidney, n.knots = 12, kappa1 = 1000)

      left truncated structure used

  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
sex -1.59913     0.202      0.4641        0.4442 -3.446 0.00057
age  0.00718     1.007      0.0119        0.0118  0.603 0.55000

    Frailty parameter, Theta: 0.432 (SE (H): 0.235 ) (SE (HIH): 0.231 ) 
 
    penalized marginal log-likelihood = -319.65
    n= 76
    n events= 58  n groups= 38
    number of iterations:  36
    Exact number of knots used:  12 
    Value of the smoothing parameter:  1000, DoF:  13.98
> 
> #stratified data. Let's use another dataset
> data(readmission)
> frailtyPenal(Surv(time,event)~as.factor(dukes)+cluster(id)+strata(sex),
+              n.knots=10,kappa1=10000,kappa2=10000,data=readmission)
Call:
frailtyPenal(formula = Surv(time, event) ~ as.factor(dukes) + 
    cluster(id) + strata(sex), data = readmission, n.knots = 10, 
    kappa1 = 10000, kappa2 = 10000)


  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
as.factor(dukes)2 0.444      1.56       0.149         0.149 2.97 2.9e-03
as.factor(dukes)3 1.287      3.62       0.182         0.182 7.06 1.7e-12

    Frailty parameter, Theta: 0.693 (SE (H): 0.145 ) (SE (HIH): 0.145 ) 
 
    penalized marginal log-likelihood = -3244.22
    n= 861
    n events= 458  n groups= 403
    number of iterations:  19
    Exact number of knots used:  10 
    Value of the smoothing parameter:  10000  10000, DoF:  11.01
> 
> #Andersen-Gill counting-process approach with time-dependent covariate
> frailtyPenal(Surv(t.start,t.stop,event)~as.factor(sex)+as.factor(dukes)+
+              as.factor(charlson)+cluster(id),data=readmission, Frail=TRUE,
+              n.knots=6,kappa1=100000,recurrentAG=TRUE)
Call:
frailtyPenal(formula = Surv(t.start, t.stop, event) ~ as.factor(sex) + 
    as.factor(dukes) + as.factor(charlson) + cluster(id), data = readmission, 
    Frailty = TRUE, recurrentAG = TRUE, n.knots = 6, kappa1 = 1e+05)

      left truncated structure used

  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
as.factor(sex)2      -0.630     0.532       0.162         0.162 -3.88 1.0e-04
as.factor(dukes)2     0.372     1.451       0.177         0.177  2.10 3.6e-02
as.factor(dukes)3     1.334     3.797       0.225         0.225  5.92 3.3e-09
as.factor(charlson)1  0.456     1.577       0.298         0.298  1.53 1.3e-01
as.factor(charlson)3  0.608     1.837       0.151         0.151  4.03 5.6e-05

    Frailty parameter, Theta: 1.28 (SE (H): 0.19 ) (SE (HIH): 0.19 ) 
 
    penalized marginal log-likelihood = -3315.09
    n= 861
    n events= 458  n groups= 403
    number of iterations:  6
    Exact number of knots used:  6 
    Value of the smoothing parameter:  1e+05, DoF:  8.00
> 
> # with the use of the cross validation approach, to find the smoothing parameter 
> frailtyPenal(Surv(t.start,t.stop,event)~as.factor(sex)+as.factor(dukes)+
+              as.factor(charlson)+cluster(id),data=readmission, Frail=TRUE,
+              n.knots=6,kappa1=5000,recurrentAG=TRUE,cross.validation=TRUE)
Call:
frailtyPenal(formula = Surv(t.start, t.stop, event) ~ as.factor(sex) + 
    as.factor(dukes) + as.factor(charlson) + cluster(id), data = readmission, 
    Frailty = TRUE, recurrentAG = TRUE, cross.validation = TRUE, 
    n.knots = 6, kappa1 = 5000)

      left truncated structure used

  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
as.factor(sex)2      -0.607     0.545       0.163         0.164 -3.73 1.9e-04
as.factor(dukes)2     0.397     1.488       0.178         0.179  2.24 2.5e-02
as.factor(dukes)3     1.361     3.901       0.225         0.225  6.05 1.4e-09
as.factor(charlson)1  0.473     1.605       0.298         0.298  1.59 1.1e-01
as.factor(charlson)3  0.607     1.835       0.151         0.151  4.02 5.7e-05

    Frailty parameter, Theta: 1.27 (SE (H): 0.190 ) (SE (HIH): 0.19 ) 
 
    penalized marginal log-likelihood = -3317
    n= 861
    n events= 458  n groups= 403
    number of iterations:  10
    Exact number of knots used:  6 
    Best smoothing parameter estimated by
       an approximated Cross validation:  1.43e+12, DoF:  5.90
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "plotfrailtyPenal"
> 
> ### * plotfrailtyPenal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.frailtyPenal
> ### Title: Plot Method for an object of class `frailtyPenal'.
> ### Aliases: plot.frailtyPenal lines.frailtyPenal
> ### Keywords: survival
> 
> ### ** Examples
> 
> 
> data(readmission)
> 
> # Let's compare shared frailty model with Cox proportional hazards model 
> mod.sha<-frailtyPenal(Surv(time,event)~as.factor(dukes)+cluster(id),
+              n.knots=10,kappa1=10000,data=readmission)
> plot(mod.sha,type="surv",conf=FALSE)
> mod.cox<-frailtyPenal(Surv(time,event)~as.factor(dukes)+cluster(id),
+              n.knots=10,kappa1=10000,data=readmission,Frailty=FALSE)
> lines(mod.cox,type="surv",conf=FALSE,col=2)
> 
> 
> 
> # Stratified model
> mod<-frailtyPenal(Surv(time,event)~as.factor(dukes)+cluster(id)+strata(sex),
+              n.knots=10,kappa1=10000,kappa2=10000,data=readmission)
> plot(mod)
> 
> # no confidence bands
> plot(mod,conf.bands=FALSE)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "printfrailtyPenal"
> 
> ### * printfrailtyPenal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: print.frailtyPenal
> ### Title: Print a Short Summary of parameter estimates of a shared gamma
> ###   frailty model
> ### Aliases: print.frailtyPenal
> ### Keywords: survival
> 
> ### ** Examples
> 
> data(kidney)
> mod<-frailtyPenal(Surv(time,status)~sex+age+cluster(id),
+      n.knots=8,kappa1=10000,data=kidney)
> print(mod)
Call:
frailtyPenal(formula = Surv(time, status) ~ sex + age + cluster(id), 
    data = kidney, n.knots = 8, kappa1 = 10000)


  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
sex -1.73262     0.177      0.5004        0.4715 -3.462 0.00054
age  0.00812     1.008      0.0124        0.0123  0.655 0.51000

    Frailty parameter, Theta: 0.499 (SE (H): 0.256 ) (SE (HIH): 0.253 ) 
 
    penalized marginal log-likelihood = -325.64
    n= 76
    n events= 58  n groups= 38
    number of iterations:  25
    Exact number of knots used:  8 
    Value of the smoothing parameter:  10000, DoF:  9.98
> 
> 
> 
> cleanEx(); ..nameEx <- "summaryfrailtyPenal"
> 
> ### * summaryfrailtyPenal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.frailtyPenal
> ### Title: summary of 'frailtyPenal'
> ### Aliases: summary.frailtyPenal print.summary.frailtyPenal
> ###   intervals.frailtyPenal
> ### Keywords: survival
> 
> ### ** Examples
> 
> 
> data(kidney)
> mod<-frailtyPenal(Surv(time,status)~age+sex+cluster(id),
+      data=kidney,n.knots=8,kappa1=1000)
> summary(mod)
         hr     95%     C.I. 
 age   1.01 (   0.98 -   1.03 ) 
 sex   0.18 (   0.07 -   0.47 ) 
> 
> # confidence interval at 99
> 
> summary(mod,level=0.99)
         hr     99%     C.I. 
 age   1.01 (   0.98 -   1.04 ) 
 sex   0.18 (   0.05 -   0.64 ) 
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  11.98 0.08 12.08 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')