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("eha-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('eha') Loading required package: survival Loading required package: splines > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "age.window" > > ### * age.window > > flush(stderr()); flush(stdout()) > > ### Name: age.window > ### Title: Age cut of survival data > ### Aliases: age.window > ### Keywords: survival > > ### ** Examples > > dat <- data.frame(enter = 0, exit = 5.731, event = 1, x = 2) > window <- c(2, 5.3) > dat.trim <- age.window(dat, window) > > > > cleanEx(); ..nameEx <- "cal.window" > > ### * cal.window > > flush(stderr()); flush(stdout()) > > ### Name: cal.window > ### Title: Calendar time cut of survival data > ### Aliases: cal.window > ### Keywords: survival > > ### ** Examples > > dat <- data.frame(enter = 0, exit = 5.731, event = 1, + birthdate = 1962.505, x = 2) > window <- c(1963, 1965) > dat.trim <- cal.window(dat, window) > > > > cleanEx(); ..nameEx <- "check.surv" > > ### * check.surv > > flush(stderr()); flush(stdout()) > > ### Name: check.surv > ### Title: Check the integrity of survival data. > ### Aliases: check.surv > ### Keywords: manip survival > > ### ** Examples > > xx <- data.frame(enter = c(0, 1), exit = c(1.5, 3), event = c(0, 1), id = + c(1,1)) > check.surv(xx$enter, xx$exit, xx$event, xx$id) [1] "1" > > > > cleanEx(); ..nameEx <- "coxreg" > > ### * coxreg > > flush(stderr()); flush(stdout()) > > ### Name: coxreg > ### Title: Cox regression > ### Aliases: coxreg > ### Keywords: survival regression > > ### ** Examples > > > dat <- data.frame(time= c(4, 3,1,1,2,2,3), + status=c(1,1,1,0,1,1,0), + x= c(0, 2,1,1,1,0,0), + sex= c(0, 0,0,0,1,1,1)) > coxreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model Call: coxreg(formula = Surv(time, status) ~ x + strata(sex), data = dat) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.625 0.802 2.231 0.329 Events 5 Total time at risk 16 Max. log. likelihood -3.3277 LR test statistic 1.09 Degrees of freedom 1 Overall p-value 0.297117 > # Same as: > rs <- risksets(Surv(dat$time, dat$status), strata = dat$sex) > coxreg( Surv(time, status) ~ x, data = dat, rs = rs) #stratified model Call: coxreg(formula = Surv(time, status) ~ x, data = dat, rs = rs) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.625 0.802 2.231 0.329 Events 5 Total time at risk 16 Max. log. likelihood -3.3277 LR test statistic 1.09 Degrees of freedom 1 Overall p-value 0.297117 > > > > cleanEx(); ..nameEx <- "coxreg.fit" > > ### * coxreg.fit > > flush(stderr()); flush(stdout()) > > ### Name: coxreg.fit > ### Title: Cox regression > ### Aliases: coxreg.fit > ### Keywords: survival regression > > ### ** Examples > > X <- as.matrix(data.frame( + x= c(0, 2,1,4,1,0,3), + sex= c(1, 0,0,0,1,1,1))) > time <- c(1,2,3,4,5,6,7) > status <- c(1,1,1,0,1,1,0) > stratum <- rep(1, length(time)) > > coxreg.fit(X, Surv(time, status), strats = stratum, max.survs = 6, + control = list(eps=1.e-4, maxiter = 10, trace = FALSE)) $coefficients [1] -0.7716698 -1.7543183 $var [,1] [,2] [1,] 0.3057497 0.3037635 [2,] 0.3037635 1.8157100 $loglik [1] -7.138867 -5.273533 $score [1] 3.274765 $linear.predictors [,1] 1 -1.7543183 2 -1.5433395 3 -0.7716698 4 -3.0866790 5 -2.5259880 6 -1.7543183 7 -4.0693275 $residuals [1] 0.8514396 0.6010801 -0.4571630 -0.1439171 0.4517570 -1.0961680 -0.2070286 $noOfRisksets [1] 5 $risktimes [1] 1 2 3 5 6 $hazard [1] 0.8586059 1.0084167 1.2853682 3.7024245 5.2600009 $means x sex 1.5714286 0.5714286 $bootstrap NULL $boot.sd NULL $conver [1] 1 $fail [1] 0 $iter [1] 5 > > > > cleanEx(); ..nameEx <- "cro" > > ### * cro > > flush(stderr()); flush(stdout()) > > ### Name: cro > ### Title: Creates a minimal representation of a data frame. > ### Aliases: cro > ### Keywords: manip > > ### ** Examples > > dat <- data.frame(y = c(1.1, 2.3, 0.7), x1 = c(1, 0, 1), x2 = c(0, 1, 0)) > cro(dat) $y [1] 1.1 2.3 0.7 $covar x1 x2 1 1 0 2 0 1 $keys [1] 1 2 1 > > > > cleanEx(); ..nameEx <- "frail.fit" > > ### * frail.fit > > flush(stderr()); flush(stdout()) > > ### Name: frail.fit > ### Title: Fits a frailty proportional hazards model > ### Aliases: frail.fit > ### Keywords: survival regression > > ### ** Examples > > fam.size <- rep(20, 5) > n <- sum(fam.size) > fam <- rep(1:length(fam.size), fam.size) > beta <- 1 > sigma <- 1 > frail <- rep(rnorm(length(fam.size), 0, sigma), fam.size) > x <- runif(n) > exit <- rexp(n, exp(x * beta + frail)) > event <- rep(1, n) > enter <- numeric(n) > mlreg(Surv(enter, exit, event) ~ x, frailty = frail) Call: mlreg(formula = Surv(enter, exit, event) ~ x, frailty = frail) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.451 0.966 2.627 0.017 Frailty standard deviation = 1.22 S.E. = 0.175 Events 100 Total time at risk 75.625 Max. log. likelihood -441.63 LR test statistic 5.64 Degrees of freedom 1 Overall p-value 0.0175679 > > > > cleanEx(); ..nameEx <- "geome.fit" > > ### * geome.fit > > flush(stderr()); flush(stdout()) > > ### Name: geome.fit > ### Title: Constant intensity discrete time proportional hazards > ### Aliases: geome.fit > ### Keywords: survival > > ### ** Examples > > > > > > cleanEx(); ..nameEx <- "make.communal" > > ### * make.communal > > flush(stderr()); flush(stdout()) > > ### Name: make.communal > ### Title: Put communals in "fixed" data frame > ### Aliases: make.communal > ### Keywords: survival > > ### ** Examples > > dat <- data.frame(enter = 0, exit = 5.731, event = 1, + birthdate = 1962.505, x = 2) > ## Birth date: July 2, 1962 (approximately). > com.dat <- data.frame(price = c(12, 3, -5, 6, -8, -9, 1, 7)) > dat.com <- make.communal(dat, com.dat, start = 1962.000) beg.per = 1962 end.per = 1970 > > > > cleanEx(); ..nameEx <- "mlreg" > > ### * mlreg > > flush(stderr()); flush(stdout()) > > ### Name: mlreg > ### Title: ML proportional hazards regression > ### Aliases: mlreg > ### Keywords: survival regression > > ### ** Examples > > > dat <- data.frame(time= c(4, 3,1,1,2,2,3), + status=c(1,1,1,0,1,1,0), + x= c(0, 2,1,1,1,0,0), + sex= c(0, 0,0,0,1,1,1)) > mlreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model Call: mlreg(formula = Surv(time, status) ~ x + strata(sex), data = dat) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.625 1.240 3.455 0.180 Events 5 Total time at risk 16 Max. log. likelihood -4.4118 LR test statistic 2.27 Degrees of freedom 1 Overall p-value 0.132173 > # Same as: > rs <- risksets(Surv(dat$time, dat$status), strata = dat$sex) > mlreg( Surv(time, status) ~ x, data = dat, rs = rs) #stratified model Call: mlreg(formula = Surv(time, status) ~ x, data = dat, rs = rs) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.625 1.240 3.455 0.180 Events 5 Total time at risk 16 Max. log. likelihood -4.4118 LR test statistic 2.27 Degrees of freedom 1 Overall p-value 0.132173 > > > > cleanEx(); ..nameEx <- "mlreg.fit" > > ### * mlreg.fit > > flush(stderr()); flush(stdout()) > > ### Name: mlreg.fit > ### Title: ML proportional hazards regression > ### Aliases: mlreg.fit > ### Keywords: survival regression > > ### ** Examples > > X <- as.matrix(data.frame( + x= c(0, 2,1,4,1,0,3), + sex= c(1, 0,0,0,1,1,1))) > time <- c(1,2,3,4,5,6,7) > status <- c(1,1,1,0,1,1,0) > stratum <- rep(1, length(time)) > > mlreg.fit(X, Surv(time, status), strats = stratum, max.survs = 6, + control = list(eps=1.e-4, maxiter = 10, trace = TRUE)) *** Iteration [1] 0 L2 = [1] 2.069264 loglik = [1] -11.37203 *** Iteration [1] 1 L2 = [1] 0.247308 loglik = [1] -8.807664 *** Iteration [1] 2 L2 = [1] 0.04405284 loglik = [1] -8.71684 *** Iteration [1] 3 L2 = [1] 0.001357857 loglik = [1] -8.71475 *** Iteration [1] 4 L2 = [1] 1.234757e-06 loglik = [1] -8.714748 *** Iteration [1] 5 Convergence loglik = [1] -8.714748 $coefficients [1] -0.9033924 -2.1480920 $var [,1] [,2] [1,] 0.2909465 0.3022695 [2,] 0.3022695 1.9457905 $loglik [1] -11.372030 -8.714748 $score [1] 4.826028 $linear.predictors [,1] 1 -2.1480920 2 -1.8067848 3 -0.9033924 4 -3.6135695 5 -3.0514843 6 -2.1480920 7 -4.8582691 $residuals [1] 0.990395574 0.970303199 0.855130680 -0.009637434 0.964992026 [6] 0.818547527 -0.012071129 $noOfRisksets [1] 5 $risktimes [1] 1 2 3 5 6 $hazard [1] 0.08229547 0.09858277 0.17665362 0.38277161 0.81447091 $means x sex 1.5714286 0.5714286 $bootstrap NULL $boot.sd NULL $conver [1] 1 $fail [1] 0 $iter [1] 5 > > > > cleanEx(); ..nameEx <- "plot.Surv" > > ### * plot.Surv > > flush(stderr()); flush(stdout()) > > ### Name: plot.Surv > ### Title: Plots of survivor functions. > ### Aliases: plot.Surv > ### Keywords: survival > > ### ** Examples > > time0 <- numeric(50) > group <- c(rep(0, 25), rep(1, 25)) > time1 <- rexp( 50, exp(group) ) > event <- rep(1, 50) > plot.Surv(Surv(time0, time1, event), strata = group) > > > > cleanEx(); ..nameEx <- "plot.cum" > > ### * plot.cum > > flush(stderr()); flush(stdout()) > > ### Name: plot.cum > ### Title: Plots of Cumulative Hazards Functions. > ### Aliases: plot.cum > ### Keywords: survival > > ### ** Examples > > time0 <- numeric(50) > group <- c(rep(0, 25), rep(1, 25)) > time1 <- rexp( 50, exp(group) ) > event <- rep(1, 50) > plot.cum(Surv(time0, time1, event), group) Warning in plot.cum(Surv(time0, time1, event), group) : 'plot.cum' is deprecated (use 'plot.Surv') > plot.cum(Surv(time0, time1, event), group, log.scale = TRUE) Warning in plot.cum(Surv(time0, time1, event), group, log.scale = TRUE) : 'plot.cum' is deprecated (use 'plot.Surv') > > > > cleanEx(); ..nameEx <- "plot.weibreg" > > ### * plot.weibreg > > flush(stderr()); flush(stdout()) > > ### Name: plot.weibreg > ### Title: Plots output from a Weibull regression > ### Aliases: plot.weibreg > ### Keywords: dplot survival > > ### ** Examples > > ## The function is currently defined as > function (x, new.data = rep(0, length(x$means))) + { + if (!inherits(x, "weibreg")) + stop("Works only with 'weibreg' objects.") + ncov <- length(x$means) + ns <- x$n.strata + lambda <- exp(x$coefficients[ncov + (1:ns) * 2 - 1]) + p <- exp(x$coefficients[ncov + (1:ns) * 2]) + xlim <- c(min(x$y[, 1]), max(x$y[, 2])) + npts <- 199 + xx <- seq(xlim[1], xlim[2], length = npts) + if (xx[1] <= 0) + xx[1] <- 0.001 + haz <- matrix(ncol = npts, nrow = ns) + for (i in 1:ns) { + tl <- xx/lambda[i] + haz[i, ] <- (p[i]/lambda[i]) * tl^(p[i] - 1) + } + ylim <- c(0, max(haz)) + cat("length(xx) = ", length(xx), "\n") + cat("dim(haz) = ", dim(haz), "\n") + plot(xx, haz[1, ], type = "l", xlim = xlim, ylim = ylim) + if (ns > 1) { + for (i in 2:ns) { + lines(xx, haz[i, ], type = "l", lty = i) + } + } + } function (x, new.data = rep(0, length(x$means))) { if (!inherits(x, "weibreg")) stop("Works only with 'weibreg' objects.") ncov <- length(x$means) ns <- x$n.strata lambda <- exp(x$coefficients[ncov + (1:ns) * 2 - 1]) p <- exp(x$coefficients[ncov + (1:ns) * 2]) xlim <- c(min(x$y[, 1]), max(x$y[, 2])) npts <- 199 xx <- seq(xlim[1], xlim[2], length = npts) if (xx[1] <= 0) xx[1] <- 0.001 haz <- matrix(ncol = npts, nrow = ns) for (i in 1:ns) { tl <- xx/lambda[i] haz[i, ] <- (p[i]/lambda[i]) * tl^(p[i] - 1) } ylim <- c(0, max(haz)) cat("length(xx) = ", length(xx), "\n") cat("dim(haz) = ", dim(haz), "\n") plot(xx, haz[1, ], type = "l", xlim = xlim, ylim = ylim) if (ns > 1) { for (i in 2:ns) { lines(xx, haz[i, ], type = "l", lty = i) } } } > > > > cleanEx(); ..nameEx <- "risksets" > > ### * risksets > > flush(stderr()); flush(stdout()) > > ### Name: risksets > ### Title: Finds the compositions and sizes of risk sets > ### Aliases: risksets > ### Keywords: survival > > ### ** Examples > > enter <- c(0, 1, 0, 0) > exit <- c(1, 2, 3, 4) > event <- c(1, 1, 1, 0) > risksets(Surv(enter, exit, event)) $ns [1] 1 $antrs [1] 3 $risktimes [1] 1 2 3 $n.events [1] 1 1 1 $size [1] 3 3 2 $eventset [1] 1 2 3 $riskset [1] 1 3 4 2 3 4 3 4 > > > > cleanEx(); ..nameEx <- "summary.coxreg" > > ### * summary.coxreg > > flush(stderr()); flush(stdout()) > > ### Name: summary.coxreg > ### Title: Prints coxreg objects > ### Aliases: summary.coxreg > ### Keywords: survival print > > ### ** Examples > > ## The function is currently defined as > function (object, ...) + print(object) function (object, ...) print(object) > > > > cleanEx(); ..nameEx <- "summary.weibreg" > > ### * summary.weibreg > > flush(stderr()); flush(stdout()) > > ### Name: summary.weibreg > ### Title: Prints a weibreg object > ### Aliases: summary.weibreg > ### Keywords: survival print > > ### ** Examples > > ## The function is currently defined as > function (object, ...) + print(object) function (object, ...) print(object) > > > > cleanEx(); ..nameEx <- "table.events" > > ### * table.events > > flush(stderr()); flush(stdout()) > > ### Name: table.events > ### Title: Calculating failure times, risk set sizes and No. of events in > ### each risk set > ### Aliases: table.events > ### Keywords: survival > > ### ** Examples > > exit = c(1,2,3,4,5) > event = c(1,1,0,1,1) > table.events(exit = exit, event = event) $times [1] 1 2 4 $events [1] 1 1 1 $riskset.sizes [1] 5 4 2 > > > > cleanEx(); ..nameEx <- "weibreg" > > ### * weibreg > > flush(stderr()); flush(stdout()) > > ### Name: weibreg > ### Title: Weibull regression > ### Aliases: weibreg > ### Keywords: survival regression > > ### ** Examples > > dat <- data.frame(time= c(4, 3,1,1,2,2,3), + status=c(1,1,1,0,1,1,0), + x= c(0, 2,1,1,1,0,0), + sex= c(0, 0,0,0,1,1,1)) > weibreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model fit$fail = 0 Call: weibreg(formula = Surv(time, status) ~ x + strata(sex), data = dat) Covariate Mean Coef Rel.Risk L-R p Wald p x 0.625 0.590 1.804 0.320 log(scale):1 0.000 1.175 3.238 0.000 log(shape):1 0.000 1.040 2.829 0.036 log(scale):2 0.000 0.930 2.533 0.000 log(shape):2 0.000 1.359 3.893 0.017 Events 5 Total time at risk 16 Max. log. likelihood -7.5616 LR test statistic 0.99 Degrees of freedom 5 Overall p-value 0.962976 > > > > ### *