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. > ### *
> ### > 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("gcmrec-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('gcmrec') Loading required package: survrec Loading required package: boot > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "addCenTime" > > ### * addCenTime > > flush(stderr()); flush(stdout()) > > ### Name: addCenTime > ### Title: Add censored time equal to 0 > ### Aliases: addCenTime > ### Keywords: survival > > ### ** Examples > > library(survival) Loading required package: splines Attaching package: 'survival' The following object(s) are masked from package:boot : aml > data(bladder2) Warning in data(bladder2) : data set 'bladder2' not found > > # we compute the interocurrence time > bladder2$time<-bladder2$stop-bladder2$start > > # If we execute: > # gcmrec(Survr(id,time,event)~rx+size+number,data=bladder2,s=2060) > > # We will obtain the following error message: > # Error in Survr(id, time, event) : Data doesn't match... > > # This means that we have some patients without right-censoring time. So, > # we understand that the last event coincides with the end of study. > # Consequently,we need to add a line with time 0 and status value equal > # to 0, too. To do so, we can use the function "addCenTime" as follows: > > # for example: > # bladder2[bladder2$id==12,] > > # id rx number size start stop event enum time > # 45 12 1 1 1 0 3 1 1 3 > # 46 12 1 1 1 3 16 1 2 13 > # 47 12 1 1 1 16 23 1 3 7 > > # there is no censored time for 12th patient. So, if we execute > > bladderOK<-addCenTime(bladder2) > > # we get > > # id rx number size start stop event enum time > # 45 12 1 1 1 0 3 1 1 3 > # 46 12 1 1 1 3 16 1 2 13 > # 47 12 1 1 1 16 23 1 3 7 > # 471 12 1 1 1 16 23 0 3 0 > > > > > cleanEx(); ..nameEx <- "gcmrec" > > ### * gcmrec > > flush(stderr()); flush(stdout()) > > ### Name: gcmrec > ### Title: General Class of Models for recurrent event data > ### Aliases: gcmrec > ### Keywords: survival > > ### ** Examples > > > > ################################### > ## Models using different data formats > ################################### > > # > # Data input as a data frame > # > > # We use the well-known bladder cancer data set from survival package > > library(survival) Loading required package: splines Attaching package: 'survival' The following object(s) are masked from package:boot : aml > data(bladder2) Warning in data(bladder2) : data set 'bladder2' not found > > # we compute the interocurrence time > bladder2$time<-bladder2$stop-bladder2$start > > # If we execute: > # gcmrec(Survr(id,time,event)~rx+size+number,data=bladder2,s=2060) > > # We will obtain the following error message: > # Error in Survr(id, time, event) : Data doesn't match... > > # This means that we have some patients without right-censoring time. So, > # we understand that the last event coincides with the end of study. > # Consequently,we need to add a line with time 0 and status value equal > # to 0, too. To do so, we can use the function "addCenTime" as follows: > > bladderOK<-addCenTime(bladder2) > > # Now, we can fit the model using this new data set: > > gcmrec(Survr(id,time,event)~rx+size+number,data=bladderOK,s=2060) Call: gcmrec(formula = Survr(id, time, event) ~ rx + size + number, data = bladderOK, s = 2060) coef exp(coef) se(coef) z p rx -0.3189 0.727 0.2051 -1.555 0.1200 size -0.0154 0.985 0.0695 -0.222 0.8200 number 0.1354 1.145 0.0511 2.649 0.0081 General class model parameter estimates rho function: Alpha to k alpha (s.e.): 0.983 (0.0737) log-likelihood= -521.42 n= 85 n times= 197 number of iterations: 5 Newton-Raphson > > > # > # Data as a list. See either GeneratedData or hydraulic data > # sets as an example. > # > > # > # We can fit the model by transforming our data in a data frame > # using "List.to.Dataframe" function: > # > > data(hydraulic) > hydraulicOK<-List.to.Dataframe(hydraulic) > gcmrec(Survr(id,time,event)~covar.1+covar.2,data=hydraulicOK,s=4753) Call: gcmrec(formula = Survr(id, time, event) ~ covar.1 + covar.2, data = hydraulicOK, s = 4753) coef exp(coef) se(coef) z p covar.1 -0.0764 0.926 0.201 -0.381 0.70 covar.2 -0.0537 0.948 0.206 -0.261 0.79 General class model parameter estimates rho function: Alpha to k alpha (s.e.): 1.03 (0.0107) log-likelihood= -612.33 n= 6 n times= 158 number of iterations: 4 Newton-Raphson > > > # > # Our model allows us to incorporate effective age information > # > # To illustrate this example, we will use a simulated data set generated > # under the minimal repair model with probability of perfect repair equal to 0.6 > # > # As we have the data in a list, first we need to obtain a data frame containing > # the time, event, and covariates information: > # > > data(GeneratedData) > temp<-List.to.Dataframe(GeneratedData) > > # then, we can fit the model incorporating the information about the effective > # age in the effageData argument: > > gcmrec(Survr(id,time,event)~covar.1+covar.2, data=temp, + effageData=GeneratedData, s=100) Call: gcmrec(formula = Survr(id, time, event) ~ covar.1 + covar.2, data = temp, effageData = GeneratedData, s = 100) coef exp(coef) se(coef) z p covar.1 0.498 1.65 0.227 2.20 2.8e-02 covar.2 1.018 2.77 0.155 6.57 5.0e-11 General class model parameter estimates rho function: Alpha to k alpha (s.e.): 0.995 (0.0163) log-likelihood= -393.61 n= 10 n times= 135 number of iterations: 7 Newton-Raphson > > > > ##################################################################### > ## How to fit minimal or perfect repair models, with and without frailties > ##################################################################### > > # Model with frailties > > mod.Fra<-gcmrec(Survr(id,time,event)~rx+size+number,data=bladderOK,s=2060,Frailty=TRUE) > print(mod.Fra) Call: gcmrec(formula = Survr(id, time, event) ~ rx + size + number, data = bladderOK, s = 2060, Frailty = TRUE) coef exp(coef) se(coef) z p rx -0.3189 0.727 NA NA NA size -0.0154 0.985 NA NA NA number 0.1354 1.145 NA NA NA General class model parameter estimates rho function: Alpha to k alpha (s.e.): 0.983 (NA) Frailty parameter, Xi (s.e. Jacknife): 3.17e+13 ( NA ) Marginal log-likelihood= -521.42 n= 85 n times= 197 number of iterations: 13 EM steps > > # effective age function: perfect repair and minimal repair models > # (models without frailties) > > data(readmission) > > # perfect > mod.per<-gcmrec(Survr(id,time,event)~as.factor(dukes),data=readmission, + s=3000,typeEffage="per") > print(mod.per) Call: gcmrec(formula = Survr(id, time, event) ~ as.factor(dukes), data = readmission, s = 3000, typeEffage = "per") coef exp(coef) se(coef) z p as.factor(dukes)2 0.365 1.44 0.112 3.26 1.1e-03 as.factor(dukes)3 0.965 2.63 0.133 7.24 4.4e-13 General class model parameter estimates rho function: Alpha to k alpha (s.e.): 1.13 (0.0143) log-likelihood= -2726.23 n= 403 n times= 861 number of iterations: 7 Newton-Raphson > > # minimal > mod.min<-gcmrec(Survr(id,time,event)~as.factor(dukes),data=readmission, + s=3000,typeEffage="min") > print(mod.min) Call: gcmrec(formula = Survr(id, time, event) ~ as.factor(dukes), data = readmission, s = 3000, typeEffage = "min") coef exp(coef) se(coef) z p as.factor(dukes)2 0.402 1.50 0.112 3.60 0.00032 as.factor(dukes)3 1.214 3.37 0.134 9.03 0.00000 General class model parameter estimates rho function: Alpha to k alpha (s.e.): 1.21 (0.0150) log-likelihood= -2448.15 n= 403 n times= 861 number of iterations: 8 Newton-Raphson > > ##################################################################### > ## How to fit models with \rho function equal to identity > ##################################################################### > > data(lymphoma) > > gcmrec(Survr(id, time, event) ~ as.factor(distrib), + data = lymphoma, s = 1000, Frailty = TRUE, rhoFunc = "Ident") Error in gcmrec(Survr(id, time, event) ~ as.factor(distrib), data = lymphoma, : NA/NaN/Inf in foreign function call (arg 19) Execution halted