File: C:/Program Files/R/R-4.4.0/library/survival/tests/Examples/survival-Ex.Rout.save
R Under development (unstable) (2024-02-07 r85873) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu
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.
Natural language support but running in an English locale
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 an HTML browser interface to help.
Type 'q()' to quit R.
> pkgname <- "survival"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('survival')
>
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("Surv")
> ### * Surv
>
> flush(stderr()); flush(stdout())
>
> ### Name: Surv
> ### Title: Create a Survival Object
> ### Aliases: Surv is.Surv [.Surv
> ### Keywords: survival
>
> ### ** Examples
>
> with(aml, Surv(time, status))
[1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+ 5 5 8 8
[16] 12 16+ 23 27 30 33 43 45
> survfit(Surv(time, status) ~ ph.ecog, data=lung)
Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung)
1 observation deleted due to missingness
n events median 0.95LCL 0.95UCL
ph.ecog=0 63 37 394 348 574
ph.ecog=1 113 82 306 268 429
ph.ecog=2 50 44 199 156 288
ph.ecog=3 1 1 118 NA NA
> Surv(heart$start, heart$stop, heart$event)
[1] ( 0.0, 50.0] ( 0.0, 6.0] ( 0.0, 1.0+] ( 1.0, 16.0]
[5] ( 0.0, 36.0+] ( 36.0, 39.0] ( 0.0, 18.0] ( 0.0, 3.0]
[9] ( 0.0, 51.0+] ( 51.0, 675.0] ( 0.0, 40.0] ( 0.0, 85.0]
[13] ( 0.0, 12.0+] ( 12.0, 58.0] ( 0.0, 26.0+] ( 26.0, 153.0]
[17] ( 0.0, 8.0] ( 0.0, 17.0+] ( 17.0, 81.0] ( 0.0, 37.0+]
[21] ( 37.0,1387.0] ( 0.0, 1.0] ( 0.0, 28.0+] ( 28.0, 308.0]
[25] ( 0.0, 36.0] ( 0.0, 20.0+] ( 20.0, 43.0] ( 0.0, 37.0]
[29] ( 0.0, 18.0+] ( 18.0, 28.0] ( 0.0, 8.0+] ( 8.0,1032.0]
[33] ( 0.0, 12.0+] ( 12.0, 51.0] ( 0.0, 3.0+] ( 3.0, 733.0]
[37] ( 0.0, 83.0+] ( 83.0, 219.0] ( 0.0, 25.0+] ( 25.0,1800.0+]
[41] ( 0.0,1401.0+] ( 0.0, 263.0] ( 0.0, 71.0+] ( 71.0, 72.0]
[45] ( 0.0, 35.0] ( 0.0, 16.0+] ( 16.0, 852.0] ( 0.0, 16.0]
[49] ( 0.0, 17.0+] ( 17.0, 77.0] ( 0.0, 51.0+] ( 51.0,1587.0+]
[53] ( 0.0, 23.0+] ( 23.0,1572.0+] ( 0.0, 12.0] ( 0.0, 46.0+]
[57] ( 46.0, 100.0] ( 0.0, 19.0+] ( 19.0, 66.0] ( 0.0, 4.5+]
[61] ( 4.5, 5.0] ( 0.0, 2.0+] ( 2.0, 53.0] ( 0.0, 41.0+]
[65] ( 41.0,1408.0+] ( 0.0, 58.0+] ( 58.0,1322.0+] ( 0.0, 3.0]
[69] ( 0.0, 2.0] ( 0.0, 40.0] ( 0.0, 1.0+] ( 1.0, 45.0]
[73] ( 0.0, 2.0+] ( 2.0, 996.0] ( 0.0, 21.0+] ( 21.0, 72.0]
[77] ( 0.0, 9.0] ( 0.0, 36.0+] ( 36.0,1142.0+] ( 0.0, 83.0+]
[81] ( 83.0, 980.0] ( 0.0, 32.0+] ( 32.0, 285.0] ( 0.0, 102.0]
[85] ( 0.0, 41.0+] ( 41.0, 188.0] ( 0.0, 3.0] ( 0.0, 10.0+]
[89] ( 10.0, 61.0] ( 0.0, 67.0+] ( 67.0, 942.0+] ( 0.0, 149.0]
[93] ( 0.0, 21.0+] ( 21.0, 343.0] ( 0.0, 78.0+] ( 78.0, 916.0+]
[97] ( 0.0, 3.0+] ( 3.0, 68.0] ( 0.0, 2.0] ( 0.0, 69.0]
[101] ( 0.0, 27.0+] ( 27.0, 842.0+] ( 0.0, 33.0+] ( 33.0, 584.0]
[105] ( 0.0, 12.0+] ( 12.0, 78.0] ( 0.0, 32.0] ( 0.0, 57.0+]
[109] ( 57.0, 285.0] ( 0.0, 3.0+] ( 3.0, 68.0] ( 0.0, 10.0+]
[113] ( 10.0, 670.0+] ( 0.0, 5.0+] ( 5.0, 30.0] ( 0.0, 31.0+]
[117] ( 31.0, 620.0+] ( 0.0, 4.0+] ( 4.0, 596.0+] ( 0.0, 27.0+]
[121] ( 27.0, 90.0] ( 0.0, 5.0+] ( 5.0, 17.0] ( 0.0, 2.0]
[125] ( 0.0, 46.0+] ( 46.0, 545.0+] ( 0.0, 21.0] ( 0.0, 210.0+]
[129] (210.0, 515.0+] ( 0.0, 67.0+] ( 67.0, 96.0] ( 0.0, 26.0+]
[133] ( 26.0, 482.0+] ( 0.0, 6.0+] ( 6.0, 445.0+] ( 0.0, 428.0+]
[137] ( 0.0, 32.0+] ( 32.0, 80.0] ( 0.0, 37.0+] ( 37.0, 334.0]
[141] ( 0.0, 5.0] ( 0.0, 8.0+] ( 8.0, 397.0+] ( 0.0, 60.0+]
[145] ( 60.0, 110.0] ( 0.0, 31.0+] ( 31.0, 370.0+] ( 0.0, 139.0+]
[149] (139.0, 207.0] ( 0.0, 160.0+] (160.0, 186.0] ( 0.0, 340.0]
[153] ( 0.0, 310.0+] (310.0, 340.0+] ( 0.0, 28.0+] ( 28.0, 265.0+]
[157] ( 0.0, 4.0+] ( 4.0, 165.0] ( 0.0, 2.0+] ( 2.0, 16.0]
[161] ( 0.0, 13.0+] ( 13.0, 180.0+] ( 0.0, 21.0+] ( 21.0, 131.0+]
[165] ( 0.0, 96.0+] ( 96.0, 109.0+] ( 0.0, 21.0] ( 0.0, 38.0+]
[169] ( 38.0, 39.0+] ( 0.0, 31.0+] ( 0.0, 11.0+] ( 0.0, 6.0]
>
>
>
> cleanEx()
> nameEx("aareg")
> ### * aareg
>
> flush(stderr()); flush(stdout())
>
> ### Name: aareg
> ### Title: Aalen's additive regression model for censored data
> ### Aliases: aareg
> ### Keywords: survival
>
> ### ** Examples
>
> # Fit a model to the lung cancer data set
> lfit <- aareg(Surv(time, status) ~ age + sex + ph.ecog, data=lung,
+ nmin=1)
> ## Not run:
> ##D lfit
> ##D Call:
> ##D aareg(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung, nmin = 1
> ##D )
> ##D
> ##D n=227 (1 observations deleted due to missing values)
> ##D 138 out of 138 unique event times used
> ##D
> ##D slope coef se(coef) z p
> ##D Intercept 5.26e-03 5.99e-03 4.74e-03 1.26 0.207000
> ##D age 4.26e-05 7.02e-05 7.23e-05 0.97 0.332000
> ##D sex -3.29e-03 -4.02e-03 1.22e-03 -3.30 0.000976
> ##D ph.ecog 3.14e-03 3.80e-03 1.03e-03 3.70 0.000214
> ##D
> ##D Chisq=26.73 on 3 df, p=6.7e-06; test weights=aalen
> ##D
> ##D plot(lfit[4], ylim=c(-4,4)) # Draw a plot of the function for ph.ecog
> ## End(Not run)
> lfit2 <- aareg(Surv(time, status) ~ age + sex + ph.ecog, data=lung,
+ nmin=1, taper=1:10)
> ## Not run: lines(lfit2[4], col=2) # Nearly the same, until the last point
>
> # A fit to the mulitple-infection data set of children with
> # Chronic Granuomatous Disease. See section 8.5 of Therneau and Grambsch.
> fita2 <- aareg(Surv(tstart, tstop, status) ~ treat + age + inherit +
+ steroids + cluster(id), data=cgd)
> ## Not run:
> ##D n= 203
> ##D 69 out of 70 unique event times used
> ##D
> ##D slope coef se(coef) robust se z p
> ##D Intercept 0.004670 0.017800 0.002780 0.003910 4.55 5.30e-06
> ##D treatrIFN-g -0.002520 -0.010100 0.002290 0.003020 -3.36 7.87e-04
> ##D age -0.000101 -0.000317 0.000115 0.000117 -2.70 6.84e-03
> ##D inheritautosomal 0.001330 0.003830 0.002800 0.002420 1.58 1.14e-01
> ##D steroids 0.004620 0.013200 0.010600 0.009700 1.36 1.73e-01
> ##D
> ##D Chisq=16.74 on 4 df, p=0.0022; test weights=aalen
> ## End(Not run)
>
>
>
> cleanEx()
> nameEx("aggregate.survfit")
> ### * aggregate.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: aggregate.survfit
> ### Title: Average survival curves
> ### Aliases: aggregate.survfit
> ### Keywords: survival
>
> ### ** Examples
>
> cfit <- coxph(Surv(futime, death) ~ sex + age*hgb, data=mgus2)
> # marginal effect of sex, after adjusting for the others
> dummy <- rbind(mgus2, mgus2)
> dummy$sex <- rep(c("F", "M"), each=nrow(mgus2)) # population data set
> dummy <- na.omit(dummy) # don't count missing hgb in our "population
> csurv <- survfit(cfit, newdata=dummy)
> dim(csurv) # 2 * 1384 survival curves
data
2676
> csurv2 <- aggregate(csurv, dummy$sex)
>
>
>
> cleanEx()
> nameEx("anova.coxph")
> ### * anova.coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: anova.coxph
> ### Title: Analysis of Deviance for a Cox model.
> ### Aliases: anova.coxph anova.coxphlist
> ### Keywords: models regression survival
>
> ### ** Examples
>
> fit <- coxph(Surv(futime, fustat) ~ resid.ds *rx + ecog.ps, data = ovarian)
> anova(fit)
Analysis of Deviance Table
Cox model: response is Surv(futime, fustat)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -34.985
resid.ds -33.105 3.7594 1 0.05251 .
rx -32.269 1.6733 1 0.19582
ecog.ps -31.970 0.5980 1 0.43934
resid.ds:rx -30.946 2.0469 1 0.15251
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> fit2 <- coxph(Surv(futime, fustat) ~ resid.ds +rx + ecog.ps, data=ovarian)
> anova(fit2,fit)
Analysis of Deviance Table
Cox model: response is Surv(futime, fustat)
Model 1: ~ resid.ds + rx + ecog.ps
Model 2: ~ resid.ds * rx + ecog.ps
loglik Chisq Df Pr(>|Chi|)
1 -31.970
2 -30.946 2.0469 1 0.1525
>
>
>
> cleanEx()
> nameEx("attrassign")
> ### * attrassign
>
> flush(stderr()); flush(stdout())
>
> ### Name: attrassign
> ### Title: Create new-style "assign" attribute
> ### Aliases: attrassign.default attrassign attrassign.lm
> ### Keywords: models
>
> ### ** Examples
>
> formula <- Surv(time,status)~factor(ph.ecog)
> tt <- terms(formula)
> mf <- model.frame(tt,data=lung)
> mm <- model.matrix(tt,mf)
> ## a few rows of data
> mm[1:3,]
(Intercept) factor(ph.ecog)1 factor(ph.ecog)2 factor(ph.ecog)3
1 1 1 0 0
2 1 0 0 0
3 1 0 0 0
> ## old-style assign attribute
> attr(mm,"assign")
[1] 0 1 1 1
> ## alternate style assign attribute
> attrassign(mm,tt)
$`(Intercept)`
[1] 1
$`factor(ph.ecog)`
[1] 2 3 4
>
>
>
> cleanEx()
> nameEx("blogit")
> ### * blogit
>
> flush(stderr()); flush(stdout())
>
> ### Name: blogit
> ### Title: Bounded link functions
> ### Aliases: blogit bcloglog bprobit blog
> ### Keywords: survival
>
> ### ** Examples
>
> py <- pseudo(survfit(Surv(time, status) ~1, lung), time=730) #2 year survival
> range(py)
[1] -0.335248 1.693831
> pfit <- glm(py ~ ph.ecog, data=lung, family=gaussian(link=blogit()))
> # For each +1 change in performance score, the odds of 2 year survival
> # are multiplied by 1/2 = exp of the coefficient.
>
>
>
> cleanEx()
> nameEx("brier")
> ### * brier
>
> flush(stderr()); flush(stdout())
>
> ### Name: brier
> ### Title: Compute the Brier score for a Cox model
> ### Aliases: brier
> ### Keywords: survival
>
> ### ** Examples
>
> cfit <- coxph(Surv(rtime, recur) ~ age + meno + size + pmin(nodes,11),
+ data= rotterdam)
> round(cfit$concordance["concordance"], 3) # some predictive power
concordance
0.675
> brier(cfit, times=c(4,6)*365.25) # values at 4 and 6 years
$rsquared
[1] 0.1501931 0.1514488
$brier
[1] 0.1952731 0.2097106
$times
[1] 1461.0 2191.5
>
>
>
> cleanEx()
> nameEx("cch")
> ### * cch
>
> flush(stderr()); flush(stdout())
>
> ### Name: cch
> ### Title: Fits proportional hazards regression model to case-cohort data
> ### Aliases: cch
> ### Keywords: survival
>
> ### ** Examples
>
> ## The complete Wilms Tumor Data
> ## (Breslow and Chatterjee, Applied Statistics, 1999)
> ## subcohort selected by simple random sampling.
> ##
>
> subcoh <- nwtco$in.subcohort
> selccoh <- with(nwtco, rel==1|subcoh==1)
> ccoh.data <- nwtco[selccoh,]
> ccoh.data$subcohort <- subcoh[selccoh]
> ## central-lab histology
> ccoh.data$histol <- factor(ccoh.data$histol,labels=c("FH","UH"))
> ## tumour stage
> ccoh.data$stage <- factor(ccoh.data$stage,labels=c("I","II","III","IV"))
> ccoh.data$age <- ccoh.data$age/12 # Age in years
>
> ##
> ## Standard case-cohort analysis: simple random subcohort
> ##
>
> fit.ccP <- cch(Surv(edrel, rel) ~ stage + histol + age, data =ccoh.data,
+ subcoh = ~subcohort, id=~seqno, cohort.size=4028)
>
>
> fit.ccP
Case-cohort analysis,x$method, Prentice
with subcohort of 668 from cohort of 4028
Call: cch(formula = Surv(edrel, rel) ~ stage + histol + age, data = ccoh.data,
subcoh = ~subcohort, id = ~seqno, cohort.size = 4028)
Coefficients:
Value SE Z p
stageII 0.73457084 0.16849620 4.359569 1.303187e-05
stageIII 0.59708356 0.17345094 3.442377 5.766257e-04
stageIV 1.38413197 0.20481982 6.757803 1.400990e-11
histolUH 1.49806307 0.15970515 9.380180 0.000000e+00
age 0.04326787 0.02373086 1.823274 6.826184e-02
>
> fit.ccSP <- cch(Surv(edrel, rel) ~ stage + histol + age, data =ccoh.data,
+ subcoh = ~subcohort, id=~seqno, cohort.size=4028, method="SelfPren")
>
> summary(fit.ccSP)
Case-cohort analysis,x$method, SelfPrentice
with subcohort of 668 from cohort of 4028
Call: cch(formula = Surv(edrel, rel) ~ stage + histol + age, data = ccoh.data,
subcoh = ~subcohort, id = ~seqno, cohort.size = 4028, method = "SelfPren")
Coefficients:
Coef HR (95% CI) p
stageII 0.736 2.088 1.501 2.905 0.000
stageIII 0.597 1.818 1.294 2.553 0.001
stageIV 1.392 4.021 2.692 6.008 0.000
histolUH 1.506 4.507 3.295 6.163 0.000
age 0.043 1.044 0.997 1.094 0.069
>
> ##
> ## (post-)stratified on instit
> ##
> stratsizes<-table(nwtco$instit)
> fit.BI<- cch(Surv(edrel, rel) ~ stage + histol + age, data =ccoh.data,
+ subcoh = ~subcohort, id=~seqno, stratum=~instit, cohort.size=stratsizes,
+ method="I.Borgan")
>
> summary(fit.BI)
Exposure-stratified case-cohort analysis, I.Borgan method.
1 2
subcohort 952 202
cohort 3622 406
Call: cch(formula = Surv(edrel, rel) ~ stage + histol + age, data = ccoh.data,
subcoh = ~subcohort, id = ~seqno, stratum = ~instit, cohort.size = stratsizes,
method = "I.Borgan")
Coefficients:
Coef HR (95% CI) p
stageII 0.737 2.090 1.501 2.909 0.000
stageIII 0.602 1.825 1.301 2.561 0.000
stageIV 1.395 4.036 2.702 6.029 0.000
histolUH 1.522 4.580 3.450 6.080 0.000
age 0.043 1.044 0.996 1.093 0.072
>
>
>
> cleanEx()
> nameEx("cipoisson")
> ### * cipoisson
>
> flush(stderr()); flush(stdout())
>
> ### Name: cipoisson
> ### Title: Confidence limits for the Poisson
> ### Aliases: cipoisson
>
> ### ** Examples
>
> cipoisson(4) # 95% confidence limit
lower upper
1.089865 10.241589
> # lower upper
> # 1.089865 10.24153
> ppois(4, 10.24153) #chance of seeing 4 or fewer events with large rate
[1] 0.02500096
> # [1] 0.02500096
> 1-ppois(3, 1.08986) #chance of seeing 4 or more, with a small rate
[1] 0.02499961
> # [1] 0.02499961
>
>
>
>
> cleanEx()
> nameEx("clogit")
> ### * clogit
>
> flush(stderr()); flush(stdout())
>
> ### Name: clogit
> ### Title: Conditional logistic regression
> ### Aliases: clogit
> ### Keywords: survival models
>
> ### ** Examples
>
> ## Not run: clogit(case ~ spontaneous + induced + strata(stratum), data=infert)
>
> # A multinomial response recoded to use clogit
> # The revised data set has one copy per possible outcome level, with new
> # variable tocc = target occupation for this copy, and case = whether
> # that is the actual outcome for each subject.
> # See the reference below for the data.
> resp <- levels(logan$occupation)
> n <- nrow(logan)
> indx <- rep(1:n, length(resp))
> logan2 <- data.frame(logan[indx,],
+ id = indx,
+ tocc = factor(rep(resp, each=n)))
> logan2$case <- (logan2$occupation == logan2$tocc)
> clogit(case ~ tocc + tocc:education + strata(id), logan2)
Call:
clogit(case ~ tocc + tocc:education + strata(id), logan2)
coef exp(coef) se(coef) z p
toccfarm -1.8964629 0.1500986 1.3807822 -1.373 0.16961
toccoperatives 1.1667502 3.2115388 0.5656465 2.063 0.03914
toccprofessional -8.1005492 0.0003034 0.6987244 -11.593 < 2e-16
toccsales -5.0292297 0.0065438 0.7700862 -6.531 6.54e-11
tocccraftsmen:education -0.3322842 0.7172835 0.0568682 -5.843 5.13e-09
toccfarm:education -0.3702858 0.6905370 0.1164100 -3.181 0.00147
toccoperatives:education -0.4222188 0.6555906 0.0584328 -7.226 4.98e-13
toccprofessional:education 0.2782469 1.3208122 0.0510212 5.454 4.94e-08
toccsales:education NA NA 0.0000000 NA NA
Likelihood ratio test=665.5 on 8 df, p=< 2.2e-16
n= 4190, number of events= 838
>
>
>
> cleanEx()
> nameEx("cluster")
> ### * cluster
>
> flush(stderr()); flush(stdout())
>
> ### Name: cluster
> ### Title: Identify clusters.
> ### Aliases: cluster
> ### Keywords: survival
>
> ### ** Examples
>
> marginal.model <- coxph(Surv(time, status) ~ rx, data= rats, cluster=litter,
+ subset=(sex=='f'))
> frailty.model <- coxph(Surv(time, status) ~ rx + frailty(litter), rats,
+ subset=(sex=='f'))
>
>
>
> cleanEx()
> nameEx("concordance")
> ### * concordance
>
> flush(stderr()); flush(stdout())
>
> ### Name: concordance
> ### Title: Compute the concordance statistic for data or a model
> ### Aliases: concordance concordance.coxph concordance.formula
> ### concordance.lm concordance.survreg
> ### Keywords: survival
>
> ### ** Examples
>
> fit1 <- coxph(Surv(ptime, pstat) ~ age + sex + mspike, mgus2)
> concordance(fit1, timewt="n/G2") # Uno's weighting
Call:
concordance.coxph(object = fit1, timewt = "n/G2")
n= 1373
Concordance= 0.6132 se= 0.1026
concordant discordant tied.x tied.y tied.xy
461425.07 290956.09 265.66 120.39 0.00
>
> # logistic regression
> fit2 <- glm(I(sex=='M') ~ age + log(creatinine), binomial, data= flchain)
> concordance(fit2) # equal to the AUC
Call:
concordance.lm(object = fit2)
n= 6524
Concordance= 0.8151 se= 0.005304
concordant discordant tied.x tied.y tied.xy
8568768 1931502 31474 10689870 56412
>
> # compare multiple models
> options(na.action = na.exclude) # predict all 1384 obs, including missing
> fit3 <- glm(pstat ~ age + sex + mspike + offset(log(ptime)),
+ poisson, data= mgus2)
> fit4 <- coxph(Surv(ptime, pstat) ~ age + sex + mspike, mgus2)
> fit5 <- coxph(Surv(ptime, pstat) ~ age + sex + hgb + creat, mgus2)
>
> tdata <- mgus2; tdata$ptime <- 60 # prediction at 60 months
> p3 <- -predict(fit3, newdata=tdata)
> p4 <- -predict(fit4) # high risk scores predict shorter survival
> p5 <- -predict(fit5)
> options(na.action = na.omit) # return to the R default
>
> cfit <- concordance(Surv(ptime, pstat) ~p3 + p4 + p5, mgus2)
> cfit
Call:
concordance.formula(object = Surv(ptime, pstat) ~ p3 + p4 + p5,
data = mgus2)
n=1338 (46 observations deleted due to missingness)
concordance se
p3 0.6598 0.0313
p4 0.6618 0.0310
p5 0.6000 0.0293
concordant discordant tied.x tied.y tied.xy
p3 51105 26333 74 28 0
p4 51258 26180 74 28 0
p5 46507 31003 2 28 0
> round(coef(cfit), 3)
p3 p4 p5
0.660 0.662 0.600
> round(cov2cor(vcov(cfit)), 3) # high correlation
[,1] [,2] [,3]
[1,] 1.000 0.994 0.236
[2,] 0.994 1.000 0.258
[3,] 0.236 0.258 1.000
>
> test <- c(1, -1, 0) # contrast vector for model 1 - model 2
> round(c(difference = test %*% coef(cfit),
+ sd= sqrt(test %*% vcov(cfit) %*% test)), 3)
difference sd
-0.002 0.003
>
>
>
> cleanEx()
> nameEx("cox.zph")
> ### * cox.zph
>
> flush(stderr()); flush(stdout())
>
> ### Name: cox.zph
> ### Title: Test the Proportional Hazards Assumption of a Cox Regression
> ### Aliases: cox.zph [.cox.zph print.cox.zph
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps,
+ data=ovarian)
> temp <- cox.zph(fit)
> print(temp) # display the results
chisq df p
age 0.698 1 0.40
ecog.ps 2.371 1 0.12
GLOBAL 3.633 2 0.16
> plot(temp) # plot curves
>
>
>
> cleanEx()
> nameEx("coxph")
> ### * coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: coxph
> ### Title: Fit Proportional Hazards Regression Model
> ### Aliases: coxph print.coxph.null print.coxph.penal coxph.penalty
> ### coxph.getdata summary.coxph.penal
> ### Keywords: survival
>
> ### ** Examples
>
> # Create the simplest test data set
> test1 <- list(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))
> # Fit a stratified model
> coxph(Surv(time, status) ~ x + strata(sex), test1)
Call:
coxph(formula = Surv(time, status) ~ x + strata(sex), data = test1)
coef exp(coef) se(coef) z p
x 0.8023 2.2307 0.8224 0.976 0.329
Likelihood ratio test=1.09 on 1 df, p=0.2971
n= 7, number of events= 5
> # Create a simple data set for a time-dependent model
> test2 <- list(start=c(1,2,5,2,1,7,3,4,8,8),
+ stop=c(2,3,6,7,8,9,9,9,14,17),
+ event=c(1,1,1,1,1,1,1,0,0,0),
+ x=c(1,0,0,1,0,1,1,1,0,0))
> summary(coxph(Surv(start, stop, event) ~ x, test2))
Call:
coxph(formula = Surv(start, stop, event) ~ x, data = test2)
n= 10, number of events= 7
coef exp(coef) se(coef) z Pr(>|z|)
x -0.02111 0.97912 0.79518 -0.027 0.979
exp(coef) exp(-coef) lower .95 upper .95
x 0.9791 1.021 0.2061 4.653
Concordance= 0.526 (se = 0.129 )
Likelihood ratio test= 0 on 1 df, p=1
Wald test = 0 on 1 df, p=1
Score (logrank) test = 0 on 1 df, p=1
>
> #
> # Create a simple data set for a time-dependent model
> #
> test2 <- list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
+ stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
+ event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
+ x =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )
>
>
> summary( coxph( Surv(start, stop, event) ~ x, test2))
Call:
coxph(formula = Surv(start, stop, event) ~ x, data = test2)
n= 10, number of events= 7
coef exp(coef) se(coef) z Pr(>|z|)
x -0.02111 0.97912 0.79518 -0.027 0.979
exp(coef) exp(-coef) lower .95 upper .95
x 0.9791 1.021 0.2061 4.653
Concordance= 0.526 (se = 0.129 )
Likelihood ratio test= 0 on 1 df, p=1
Wald test = 0 on 1 df, p=1
Score (logrank) test = 0 on 1 df, p=1
>
> # Fit a stratified model, clustered on patients
>
> bladder1 <- bladder[bladder$enum < 5, ]
> coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum),
+ cluster = id, bladder1)
Call:
coxph(formula = Surv(stop, event) ~ (rx + size + number) * strata(enum),
data = bladder1, cluster = id)
coef exp(coef) se(coef) robust se z p
rx -0.52598 0.59097 0.31583 0.31524 -1.669 0.09521
size 0.06961 1.07209 0.10156 0.08863 0.785 0.43220
number 0.23818 1.26894 0.07588 0.07459 3.193 0.00141
rx:strata(enum)enum=2 -0.10633 0.89913 0.50424 0.33396 -0.318 0.75019
rx:strata(enum)enum=3 -0.17251 0.84155 0.55780 0.39868 -0.433 0.66523
rx:strata(enum)enum=4 -0.10945 0.89632 0.65730 0.50636 -0.216 0.82886
size:strata(enum)enum=2 -0.14737 0.86298 0.16803 0.11409 -1.292 0.19646
size:strata(enum)enum=3 -0.28345 0.75318 0.20894 0.15220 -1.862 0.06255
size:strata(enum)enum=4 -0.27607 0.75876 0.25222 0.18904 -1.460 0.14418
number:strata(enum)enum=2 -0.10125 0.90370 0.11904 0.11759 -0.861 0.38920
number:strata(enum)enum=3 -0.06467 0.93738 0.12925 0.12035 -0.537 0.59101
number:strata(enum)enum=4 0.09429 1.09888 0.14594 0.11973 0.788 0.43097
Likelihood ratio test=30.09 on 12 df, p=0.002708
n= 340, number of events= 112
>
> # Fit a time transform model using current age
> coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,
+ tt=function(x,t,...) pspline(x + t/365.25))
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + tt(age), data = lung,
tt = function(x, t, ...) pspline(x + t/365.25))
coef se(coef) se2 Chisq DF p
ph.ecog 0.4528 0.1178 0.1174 14.7704 1.00 0.00012
tt(age), linear 0.0112 0.0093 0.0093 1.4414 1.00 0.22991
tt(age), nonlin 2.6992 3.08 0.45431
Iterations: 4 outer, 10 Newton-Raphson
Theta= 0.796
Degrees of freedom for terms= 1.0 4.1
Likelihood ratio test=22.5 on 5.07 df, p=5e-04
n= 227, number of events= 164
(1 observation deleted due to missingness)
>
>
>
> cleanEx()
> nameEx("coxph.detail")
> ### * coxph.detail
>
> flush(stderr()); flush(stdout())
>
> ### Name: coxph.detail
> ### Title: Details of a Cox Model Fit
> ### Aliases: coxph.detail
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- coxph(Surv(futime,fustat) ~ age + rx + ecog.ps, ovarian, x=TRUE)
> fitd <- coxph.detail(fit)
> # There is one Schoenfeld residual for each unique death. It is a
> # vector (covariates for the subject who died) - (weighted mean covariate
> # vector at that time). The weighted mean is defined over the subjects
> # still at risk, with exp(X beta) as the weight.
>
> events <- fit$y[,2]==1
> etime <- fit$y[events,1] #the event times --- may have duplicates
> indx <- match(etime, fitd$time)
> schoen <- fit$x[events,] - fitd$means[indx,]
>
>
>
> cleanEx()
> nameEx("diabetic")
> ### * diabetic
>
> flush(stderr()); flush(stdout())
>
> ### Name: diabetic
> ### Title: Ddiabetic retinopathy
> ### Aliases: diabetic
> ### Keywords: datasets survival
>
> ### ** Examples
>
> # juvenile diabetes is defined as and age less than 20
> juvenile <- 1*(diabetic$age < 20)
> coxph(Surv(time, status) ~ trt + juvenile, cluster= id,
+ data= diabetic)
Call:
coxph(formula = Surv(time, status) ~ trt + juvenile, data = diabetic,
cluster = id)
coef exp(coef) se(coef) robust se z p
trt -0.77893 0.45890 0.16893 0.14851 -5.245 1.56e-07
juvenile -0.05388 0.94754 0.16211 0.17864 -0.302 0.763
Likelihood ratio test=22.48 on 2 df, p=1.312e-05
n= 394, number of events= 155
>
>
>
> cleanEx()
> nameEx("dsurvreg")
> ### * dsurvreg
>
> flush(stderr()); flush(stdout())
>
> ### Name: dsurvreg
> ### Title: Distributions available in survreg.
> ### Aliases: dsurvreg psurvreg qsurvreg rsurvreg
> ### Keywords: distribution
>
> ### ** Examples
>
> # List of distributions available
> names(survreg.distributions)
[1] "extreme" "logistic" "gaussian" "weibull" "exponential"
[6] "rayleigh" "loggaussian" "lognormal" "loglogistic" "t"
> ## Not run:
> ##D [1] "extreme" "logistic" "gaussian" "weibull" "exponential"
> ##D [6] "rayleigh" "loggaussian" "lognormal" "loglogistic" "t"
> ## End(Not run)
> # Compare results
> all.equal(dsurvreg(1:10, 2, 5, dist='lognormal'), dlnorm(1:10, 2, 5))
[1] TRUE
>
> # Hazard function for a Weibull distribution
> x <- seq(.1, 3, length=30)
> haz <- dsurvreg(x, 2, 3)/ (1-psurvreg(x, 2, 3))
> ## Not run:
> ##D plot(x, haz, log='xy', ylab="Hazard") #line with slope (1/scale -1)
> ## End(Not run)
>
> # Estimated CDF of a simple Weibull
> fit <- survreg(Surv(time, status) ~ 1, data=lung)
> pp <- 1:99/100
> q1 <- qsurvreg(pp, coef(fit), fit$scale)
> q2 <- qweibull(pp, shape= 1/fit$scale, scale= exp(coef(fit)))
> all.equal(q1, q2)
[1] TRUE
> ## Not run:
> ##D plot(q1, pp, type='l', xlab="Months", ylab="CDF")
> ## End(Not run)
> # per the help page for dweibull, the mean is scale * gamma(1 + 1/shape)
> c(mean = exp(coef(fit))* gamma(1 + fit$scale))
mean.(Intercept)
384.8529
>
>
>
>
> cleanEx()
> nameEx("finegray")
> ### * finegray
>
> flush(stderr()); flush(stdout())
>
> ### Name: finegray
> ### Title: Create data for a Fine-Gray model
> ### Aliases: finegray
> ### Keywords: survival
>
> ### ** Examples
>
> # Treat time to death and plasma cell malignancy as competing risks
> etime <- with(mgus2, ifelse(pstat==0, futime, ptime))
> event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
> event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))
>
> # FG model for PCM
> pdata <- finegray(Surv(etime, event) ~ ., data=mgus2)
> fgfit <- coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex,
+ weight=fgwt, data=pdata)
>
> # Compute the weights separately by sex
> adata <- finegray(Surv(etime, event) ~ . + strata(sex),
+ data=mgus2, na.action=na.pass)
>
>
>
> cleanEx()
> nameEx("flchain")
> ### * flchain
>
> flush(stderr()); flush(stdout())
>
> ### Name: flchain
> ### Title: Assay of serum free light chain for 7874 subjects.
> ### Aliases: flchain
> ### Keywords: datasets
>
> ### ** Examples
>
> data(flchain)
> age.grp <- cut(flchain$age, c(49,54, 59,64, 69,74,79, 89, 110),
+ labels= paste(c(50,55,60,65,70,75,80,90),
+ c(54,59,64,69,74,79,89,109), sep='-'))
> table(flchain$sex, age.grp)
age.grp
50-54 55-59 60-64 65-69 70-74 75-79 80-89 90-109
F 881 766 625 589 541 408 459 81
M 796 714 591 524 405 269 202 23
>
>
>
> cleanEx()
> nameEx("frailty")
> ### * frailty
>
> flush(stderr()); flush(stdout())
>
> ### Name: frailty
> ### Title: Random effects terms
> ### Aliases: frailty frailty.gamma frailty.gaussian frailty.t
> ### Keywords: survival
>
> ### ** Examples
>
> # Random institutional effect
> coxph(Surv(time, status) ~ age + frailty(inst, df=4), lung)
Call:
coxph(formula = Surv(time, status) ~ age + frailty(inst, df = 4),
data = lung)
coef se(coef) se2 Chisq DF p
age 0.01937 0.00933 0.00925 4.31149 1.00 0.038
frailty(inst, df = 4) 3.33459 3.99 0.501
Iterations: 3 outer, 10 Newton-Raphson
Variance of random effect= 0.038 I-likelihood = -743.6
Degrees of freedom for terms= 1 4
Likelihood ratio test=9.96 on 4.97 df, p=0.08
n= 227, number of events= 164
(1 observation deleted due to missingness)
>
> # Litter effects for the rats data
> rfit2a <- coxph(Surv(time, status) ~ rx +
+ frailty.gaussian(litter, df=13, sparse=FALSE), rats,
+ subset= (sex=='f'))
> rfit2b <- coxph(Surv(time, status) ~ rx +
+ frailty.gaussian(litter, df=13, sparse=TRUE), rats,
+ subset= (sex=='f'))
>
>
>
> cleanEx()
> nameEx("hoel")
> ### * hoel
>
> flush(stderr()); flush(stdout())
>
> ### Name: hoel
> ### Title: Mouse cancer data
> ### Aliases: hoel
> ### Keywords: datasets
>
> ### ** Examples
>
> hsurv <- survfit(Surv(days, outcome) ~ trt, data = hoel, id= id)
> plot(hsurv, lty=1:2, col=rep(1:3, each=2), lwd=2, xscale=30.5,
+ xlab="Months", ylab= "Death")
> legend("topleft", c("Lymphoma control", "Lymphoma germ free",
+ "Sarcoma control", "Sarcoma germ free",
+ "Other control", "Other germ free"),
+ col=rep(1:3, each=2), lty=1:2, lwd=2, bty='n')
> hfit <- coxph(Surv(days, outcome) ~ trt, data= hoel, id = id)
>
>
>
> cleanEx()
> nameEx("is.ratetable")
> ### * is.ratetable
>
> flush(stderr()); flush(stdout())
>
> ### Name: is.ratetable
> ### Title: Verify that an object is of class ratetable.
> ### Aliases: is.ratetable Math.ratetable Ops.ratetable
> ### Keywords: survival
>
> ### ** Examples
>
> is.ratetable(survexp.us) # True
[1] TRUE
> is.ratetable(lung) # False
[1] FALSE
>
>
>
> cleanEx()
> nameEx("kidney")
> ### * kidney
>
> flush(stderr()); flush(stdout())
>
> ### Name: kidney
> ### Title: Kidney catheter data
> ### Aliases: kidney
> ### Keywords: survival
>
> ### ** Examples
>
> kfit <- coxph(Surv(time, status)~ age + sex + disease + frailty(id), kidney)
> kfit0 <- coxph(Surv(time, status)~ age + sex + disease, kidney)
> kfitm1 <- coxph(Surv(time,status) ~ age + sex + disease +
+ frailty(id, dist='gauss'), kidney)
>
>
>
> cleanEx()
> nameEx("levels.Surv")
> ### * levels.Surv
>
> flush(stderr()); flush(stdout())
>
> ### Name: levels.Surv
> ### Title: Return the states of a multi-state Surv object
> ### Aliases: levels.Surv
> ### Keywords: survival
>
> ### ** Examples
>
> y1 <- Surv(c(1,5, 9, 17,21, 30),
+ factor(c(0, 1, 2,1,0,2), 0:2, c("censored", "progression", "death")))
> levels(y1)
[1] "progression" "death"
>
> y2 <- Surv(1:6, rep(0:1, 3))
> y2
[1] 1+ 2 3+ 4 5+ 6
> levels(y2)
NULL
>
>
>
> cleanEx()
> nameEx("lines.survfit")
> ### * lines.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: lines.survfit
> ### Title: Add Lines or Points to a Survival Plot
> ### Aliases: lines.survfit points.survfit lines.survexp
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- survfit(Surv(time, status==2) ~ sex, pbc,subset=1:312)
> plot(fit, mark.time=FALSE, xscale=365.25,
+ xlab='Years', ylab='Survival')
> lines(fit[1], lwd=2) #darken the first curve and add marks
>
>
> # Add expected survival curves for the two groups,
> # based on the US census data
> # The data set does not have entry date, use the midpoint of the study
> efit <- survexp(~sex, data=pbc, times= (0:24)*182, ratetable=survexp.us,
+ rmap=list(sex=sex, age=age*365.35, year=as.Date('1979/01/01')))
> temp <- lines(efit, lty=2, lwd=2:1)
> text(temp, c("Male", "Female"), adj= -.1) #labels just past the ends
> title(main="Primary Biliary Cirrhosis, Observed and Expected")
>
>
>
>
> cleanEx()
> nameEx("mgus")
> ### * mgus
>
> flush(stderr()); flush(stdout())
>
> ### Name: mgus
> ### Title: Monoclonal gammopathy data
> ### Aliases: mgus mgus1
> ### Keywords: datasets survival
>
> ### ** Examples
>
> # Create the competing risk curves for time to first of death or PCM
> sfit <- survfit(Surv(start, stop, event) ~ sex, mgus1, id=id,
+ subset=(enum==1))
> print(sfit) # the order of printout is the order in which they plot
Call: survfit(formula = Surv(start, stop, event) ~ sex, data = mgus1,
subset = (enum == 1), id = id)
n nevent rmean*
sex=female, (s0) 104 0 5762.379
sex=male, (s0) 137 0 4543.293
sex=female, pcm 104 33 2881.500
sex=male, pcm 137 31 2478.026
sex=female, death 104 63 5681.121
sex=male, death 137 100 7303.682
*restricted mean time in state (max time = 14325 )
>
> plot(sfit, xscale=365.25, lty=c(2,2,1,1), col=c(1,2,1,2),
+ xlab="Years after MGUS detection", ylab="Proportion")
> legend(0, .8, c("Death/male", "Death/female", "PCM/male", "PCM/female"),
+ lty=c(1,1,2,2), col=c(2,1,2,1), bty='n')
>
> title("Curves for the first of plasma cell malignancy or death")
> # The plot shows that males have a higher death rate than females (no
> # surprise) but their rates of conversion to PCM are essentially the same.
>
>
>
> cleanEx()
> nameEx("model.matrix.coxph")
> ### * model.matrix.coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: model.matrix.coxph
> ### Title: Model.matrix method for coxph models
> ### Aliases: model.matrix.coxph
> ### Keywords: survival
>
> ### ** Examples
>
> fit1 <- coxph(Surv(time, status) ~ age + factor(ph.ecog), data=lung)
> xfit <- model.matrix(fit1)
>
> fit2 <- coxph(Surv(time, status) ~ age + factor(ph.ecog), data=lung,
+ x=TRUE)
> all.equal(model.matrix(fit1), fit2$x)
[1] TRUE
>
>
>
> cleanEx()
> nameEx("myeloid")
> ### * myeloid
>
> flush(stderr()); flush(stdout())
>
> ### Name: myeloid
> ### Title: Acute myeloid leukemia
> ### Aliases: myeloid
> ### Keywords: datasets
>
> ### ** Examples
>
> coxph(Surv(futime, death) ~ trt + flt3, data=myeloid)
Call:
coxph(formula = Surv(futime, death) ~ trt + flt3, data = myeloid)
coef exp(coef) se(coef) z p
trtB -0.3534 0.7023 0.1122 -3.149 0.00164
flt3B 0.4114 1.5089 0.1587 2.593 0.00952
flt3C 0.7878 2.1985 0.1656 4.758 1.96e-06
Likelihood ratio test=34.04 on 3 df, p=1.94e-07
n= 646, number of events= 320
> # See the mstate vignette for a more complete analysis
>
>
>
> cleanEx()
> nameEx("myeloma")
> ### * myeloma
>
> flush(stderr()); flush(stdout())
>
> ### Name: myeloma
> ### Title: Survival times of patients with multiple myeloma
> ### Aliases: myeloma
> ### Keywords: datasets
>
> ### ** Examples
>
> # Incorrect survival curve, which ignores left truncation
> fit1 <- survfit(Surv(futime, death) ~ 1, myeloma)
> # Correct curve
> fit2 <- survfit(Surv(entry, futime, death) ~1, myeloma)
>
>
>
> cleanEx()
> nameEx("neardate")
> ### * neardate
>
> flush(stderr()); flush(stdout())
>
> ### Name: neardate
> ### Title: Find the index of the closest value in data set 2, for each
> ### entry in data set one.
> ### Aliases: neardate
> ### Keywords: manip utilities
>
> ### ** Examples
>
> data1 <- data.frame(id = 1:10,
+ entry.dt = as.Date(paste("2011", 1:10, "5", sep='-')))
> temp1 <- c(1,4,5,1,3,6,9, 2,7,8,12,4,6,7,10,12,3)
> data2 <- data.frame(id = c(1,1,1,2,2,4,4,5,5,5,6,8,8,9,10,10,12),
+ lab.dt = as.Date(paste("2011", temp1, "1", sep='-')),
+ chol = round(runif(17, 130, 280)))
>
> #first cholesterol on or after enrollment
> indx1 <- neardate(data1$id, data2$id, data1$entry.dt, data2$lab.dt)
> data2[indx1, "chol"]
[1] 186 160 NA 265 224 161 NA NA NA 205
>
> # Closest one, either before or after.
> #
> indx2 <- neardate(data1$id, data2$id, data1$entry.dt, data2$lab.dt,
+ best="prior")
> ifelse(is.na(indx1), indx2, # none after, take before
+ ifelse(is.na(indx2), indx1, #none before
+ ifelse(abs(data2$lab.dt[indx2]- data1$entry.dt) <
+ abs(data2$lab.dt[indx1]- data1$entry.dt), indx2, indx1)))
[1] 1 5 NA 6 9 11 NA 13 14 15
>
> # closest date before or after, but no more than 21 days prior to index
> indx2 <- ifelse((data1$entry.dt - data2$lab.dt[indx2]) >21, NA, indx2)
> ifelse(is.na(indx1), indx2, # none after, take before
+ ifelse(is.na(indx2), indx1, #none before
+ ifelse(abs(data2$lab.dt[indx2]- data1$entry.dt) <
+ abs(data2$lab.dt[indx1]- data1$entry.dt), indx2, indx1)))
[1] 1 5 NA 6 9 11 NA NA NA 15
>
>
>
> cleanEx()
> nameEx("nsk")
> ### * nsk
>
> flush(stderr()); flush(stdout())
>
> ### Name: nsk
> ### Title: Natural splines with knot heights as the basis.
> ### Aliases: nsk
> ### Keywords: smooth
>
> ### ** Examples
>
> # make some dummy data
> tdata <- data.frame(x= lung$age, y = 10*log(lung$age-35) + rnorm(228, 0, 2))
> fit1 <- lm(y ~ -1 + nsk(x, df=4, intercept=TRUE) , data=tdata)
> fit2 <- lm(y ~ nsk(x, df=3), data=tdata)
>
> # the knots (same for both fits)
> knots <- unlist(attributes(fit1$model[[2]])[c('Boundary.knots', 'knots')])
> sort(unname(knots))
[1] 45.35 59.00 67.00 75.00
> unname(coef(fit1)) # predictions at the knot points
[1] 22.23565 32.18784 34.82994 36.93020
>
> unname(coef(fit1)[-1] - coef(fit1)[1]) # differences: yhat[2:4] - yhat[1]
[1] 9.952189 12.594290 14.694549
> unname(coef(fit2))[-1] # ditto
[1] 9.952189 12.594290 14.694549
>
> ## Not run:
> ##D plot(y ~ x, data=tdata)
> ##D points(sort(knots), coef(fit1), col=2, pch=19)
> ##D coef(fit)[1] + c(0, coef(fit)[-1])
> ## End(Not run)
>
>
>
> cleanEx()
> nameEx("nwtco")
> ### * nwtco
>
> flush(stderr()); flush(stdout())
>
> ### Name: nwtco
> ### Title: Data from the National Wilm's Tumor Study
> ### Aliases: nwtco
> ### Keywords: datasets
>
> ### ** Examples
>
> with(nwtco, table(instit,histol))
histol
instit 1 2
1 3493 129
2 76 330
> anova(coxph(Surv(edrel,rel)~histol+instit,data=nwtco))
Analysis of Deviance Table
Cox model: response is Surv(edrel, rel)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -4666.3
histol -4532.5 267.6667 1 < 2e-16 ***
instit -4531.0 3.0397 1 0.08125 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(coxph(Surv(edrel,rel)~instit+histol,data=nwtco))
Analysis of Deviance Table
Cox model: response is Surv(edrel, rel)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -4666.3
instit -4577.5 177.714 1 < 2.2e-16 ***
histol -4531.0 92.992 1 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> cleanEx()
> nameEx("pbcseq")
> ### * pbcseq
>
> flush(stderr()); flush(stdout())
>
> ### Name: pbcseq
> ### Title: Mayo Clinic Primary Biliary Cirrhosis, sequential data
> ### Aliases: pbcseq
> ### Keywords: datasets
>
> ### ** Examples
>
> # Create the start-stop-event triplet needed for coxph
> first <- with(pbcseq, c(TRUE, diff(id) !=0)) #first id for each subject
> last <- c(first[-1], TRUE) #last id
>
> time1 <- with(pbcseq, ifelse(first, 0, day))
> time2 <- with(pbcseq, ifelse(last, futime, c(day[-1], 0)))
> event <- with(pbcseq, ifelse(last, status, 0))
>
> fit1 <- coxph(Surv(time1, time2, event) ~ age + sex + log(bili), pbcseq)
Warning in Surv(time1, time2, event) :
Invalid status value, converted to NA
>
>
>
> cleanEx()
> nameEx("plot.cox.zph")
> ### * plot.cox.zph
>
> flush(stderr()); flush(stdout())
>
> ### Name: plot.cox.zph
> ### Title: Graphical Test of Proportional Hazards
> ### Aliases: plot.cox.zph
> ### Keywords: survival
>
> ### ** Examples
>
> vfit <- coxph(Surv(time,status) ~ trt + factor(celltype) +
+ karno + age, data=veteran, x=TRUE)
> temp <- cox.zph(vfit)
> plot(temp, var=3) # Look at Karnofsy score, old way of doing plot
> plot(temp[3]) # New way with subscripting
> abline(0, 0, lty=3)
> # Add the linear fit as well
> abline(lm(temp$y[,3] ~ temp$x)$coefficients, lty=4, col=3)
> title(main="VA Lung Study")
>
>
>
> cleanEx()
> nameEx("plot.survfit")
> ### * plot.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: plot.survfit
> ### Title: Plot method for 'survfit' objects
> ### Aliases: plot.survfit
> ### Keywords: survival hplot
>
> ### ** Examples
>
> leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml)
> plot(leukemia.surv, lty = 2:3)
> legend(100, .9, c("Maintenance", "No Maintenance"), lty = 2:3)
> title("Kaplan-Meier Curves\nfor AML Maintenance Study")
> lsurv2 <- survfit(Surv(time, status) ~ x, aml, type='fleming')
> plot(lsurv2, lty=2:3, fun="cumhaz",
+ xlab="Months", ylab="Cumulative Hazard")
>
>
>
> cleanEx()
> nameEx("predict.coxph")
> ### * predict.coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: predict.coxph
> ### Title: Predictions for a Cox model
> ### Aliases: predict.coxph predict.coxph.penal
> ### Keywords: survival
>
> ### ** Examples
>
> options(na.action=na.exclude) # retain NA in predictions
> fit <- coxph(Surv(time, status) ~ age + ph.ecog + strata(inst), lung)
> #lung data set has status coded as 1/2
> mresid <- (lung$status-1) - predict(fit, type='expected') #Martingale resid
> predict(fit,type="lp")
1 2 3 4 5 6
0.215495605 -0.423532231 -0.559265038 0.183469551 -0.539432878 0.248095483
7 8 9 10 11 12
0.406461814 0.489169379 -0.047448917 0.327284344 0.040389888 0.550315552
13 14 15 16 17 18
-0.115925255 NA 0.055807340 0.110906025 0.050567124 0.493760215
19 20 21 22 23 24
0.557645717 -0.004245606 -0.127236322 -0.621260082 -0.319524466 -0.575882288
25 26 27 28 29 30
-0.345688084 0.202851214 -0.428371074 1.313400384 -0.021210624 0.761244928
31 32 33 34 35 36
0.191540147 0.749933860 0.180240469 0.459827013 0.672213041 0.625512121
37 38 39 40 41 42
0.565173220 0.085767683 0.761244928 0.076972823 0.330513426 0.511791514
43 44 45 46 47 48
-0.439682141 0.660901974 -0.164699618 0.496950353 -0.381077937 0.091073865
49 50 51 52 53 54
-0.354839644 -0.175654221 0.192873470 -0.447487689 -0.450985298 -0.562055013
55 56 57 58 59 60
0.063012023 -0.516810744 -0.297203343 0.474684682 0.034518529 0.076972823
61 62 63 64 65 66
0.678283893 -0.045992266 0.176731471 -0.149858457 0.158940268 0.718790633
67 68 69 70 71 72
0.539004484 -0.308514410 -0.543216443 0.153500561 -0.479261384 -0.078592144
73 74 75 76 77 78
0.946919127 -0.073531430 -0.049489875 0.214162281 -0.641232484 0.029078821
79 80 81 82 83 84
-0.276488357 -0.392389004 -0.439682141 0.001411510 -0.410013004 -0.151289480
85 86 87 88 89 90
-0.292311495 0.198744830 -0.039921414 -0.530162769 -0.123010230 0.738622793
91 92 93 94 95 96
-0.743642023 0.050567124 0.285269157 0.108857156 -0.437633273 0.796634781
97 98 99 100 101 102
0.158940268 0.214162281 -0.161169524 -0.400910096 -0.562055013 0.176122695
103 104 105 106 107 108
0.012722577 0.108256292 0.617817211 0.157606945 -0.189452466 0.110906025
109 110 111 112 113 114
-0.026867740 0.797968104 -0.411394980 -0.149248522 0.369011703 -0.344354760
115 116 117 118 119 120
0.006456686 0.783867062 0.503880355 0.693378524 0.527693417 0.244122624
121 122 123 124 125 126
-0.464038972 0.449575370 0.158940268 0.500480446 -0.426322206 0.005322855
127 128 129 130 131 132
-0.368298829 0.134984810 0.652115157 -0.617153698 0.131479291 -0.190511890
133 134 135 136 137 138
-0.643882217 0.001411510 -0.460255408 0.666972826 0.067118407 0.583884010
139 140 141 142 143 144
-0.036137850 -0.399002948 0.747892903 0.215495605 0.630552446 0.088283890
145 146 147 148 149 150
-0.240346995 -0.200763533 -0.558074111 -0.179200822 -0.232577411 -0.524505653
151 152 153 154 155 156
0.171077519 -0.633704981 -0.331136545 -0.190511890 0.477441161 NA
157 158 159 160 161 162
-0.031097524 0.736573925 0.123673743 -0.013515715 -0.585704233 -0.038186718
163 164 165 166 167 168
0.466547245 0.108256292 -0.209943887 -0.716429053 -0.206413793 -0.699828778
169 170 171 172 173 174
0.085634157 -0.424865554 0.069277914 -0.441093652 0.107445646 -0.874783994
175 176 177 178 179 180
-0.047448917 0.046655779 0.557645717 0.001411510 -0.047448917 -0.667994646
181 182 183 184 185 186
-0.513194586 -0.776965291 -0.614629447 0.019390401 -0.583220496 -0.651086900
187 188 189 190 191 192
0.859584155 -0.536642904 0.063145548 -0.712882451 0.024398388 0.369338475
193 194 195 196 197 198
-0.023370131 0.076972823 0.061878192 -0.368310218 -0.003231734 0.074931865
199 200 201 202 203 204
-0.629921417 -0.037164935 0.063145548 0.084500326 -0.574393166 -0.627131442
205 206 207 208 209 210
-0.658814293 0.302547317 -0.410314015 0.516017606 0.131487202 -0.302547317
211 212 213 214 215 216
-0.539432878 0.153500561 0.119700884 0.409991908 -0.149858457 -0.149858457
217 218 219 220 221 222
-0.156943432 0.781826105 0.477858312 -0.452404719 0.016633922 -0.081992053
223 224 225 226 227 228
0.212705630 0.224016697 -0.750726998 0.703662506 0.142189494 -0.085165683
> predict(fit,type="expected")
1 2 3 4 5 6 7
0.74602570 0.57892506 1.28411487 0.65144995 2.53474317 2.59935704 0.94925558
8 9 10 11 12 13 14
1.07812821 0.63137435 0.55866807 0.31809979 1.96068120 2.96879741 NA
15 16 17 18 19 20 21
2.14464916 0.39248100 1.01652225 2.53985878 0.23734050 0.15454932 0.41781121
22 23 24 25 26 27 28
0.03725873 1.07425239 0.73304358 0.71922541 1.96068538 0.91425760 0.50868712
29 30 31 32 33 34 35
1.07651355 0.10727131 1.64348011 0.22335391 1.34246079 0.18355514 0.25427967
36 37 38 39 40 41 42
0.57948554 3.87217595 1.42062915 0.50341133 2.84274107 1.90670187 0.39302876
43 44 45 46 47 48 49
1.67374788 0.56009982 1.95081502 0.39930277 0.62185372 1.18384892 1.08920268
50 51 52 53 54 55 56
1.36922169 2.72429090 0.31557423 0.04821232 0.41960993 3.07164840 0.12000994
57 58 59 60 61 62 63
0.07406041 0.17908976 1.74520134 1.10195998 1.47697029 0.54523697 0.51461493
64 65 66 67 68 69 70
0.14292300 0.18117365 0.20227027 0.70028855 1.00636733 0.31133532 0.64126839
71 72 73 74 75 76 77
0.96177399 0.46743320 0.53451717 0.16345589 0.86294287 1.44797843 1.06953116
78 79 80 81 82 83 84
1.19014609 0.03668315 0.33061179 1.90397464 0.08944145 0.20857044 0.28585781
85 86 87 88 89 90 91
1.15723874 0.87295638 1.19851949 0.14216346 1.37338069 0.92021616 1.05096221
92 93 94 95 96 97 98
0.27465006 0.47403241 0.26750987 1.01622540 0.08901343 0.32456045 0.93961618
99 100 101 102 103 104 105
0.85179714 0.14362313 0.89733451 1.74403467 0.70225748 0.15754565 0.36065915
106 107 108 109 110 111 112
0.41227011 0.29089093 0.02759911 2.54485283 1.57705739 0.02915789 0.51482474
113 114 115 116 117 118 119
1.51254632 0.24392791 1.95773713 0.16855572 0.69132758 2.65613080 1.04014324
120 121 122 123 124 125 126
0.89157179 0.40187641 0.23829273 1.56065440 0.17535194 1.02778525 0.18442460
127 128 129 130 131 132 133
0.08051722 0.20596405 1.70473379 0.86354367 0.72017118 0.27146814 0.48487446
134 135 136 137 138 139 140
1.10114414 0.51567846 1.46035831 0.93950468 1.54314328 1.12143879 0.60372302
141 142 143 144 145 146 147
1.46022571 0.88081136 0.66047105 0.18347489 0.51981101 0.28761918 0.50825077
148 149 150 151 152 153 154
0.15268490 0.06671446 0.32571666 0.39746179 0.39772440 0.38939509 0.20940447
155 156 157 158 159 160 161
0.62171971 NA 0.34080256 0.46159657 0.47539058 1.00662370 0.21472196
162 163 164 165 166 167 168
0.54619593 0.50111574 0.24481910 0.51248548 0.19954882 0.25566706 0.78817717
169 170 171 172 173 174 175
0.44798249 0.43113659 0.44847984 1.48341994 0.46620310 0.37028208 0.86812344
176 177 178 179 180 181 182
0.43844817 0.94494334 0.25935783 0.37625255 0.20649507 0.25048304 0.37569346
183 184 185 186 187 188 189
0.40334526 0.39324727 0.36799524 0.39552828 1.77501387 0.24422514 0.38021709
190 191 192 193 194 195 196
0.21501843 0.51818689 0.08032921 0.22774986 0.71502728 0.36774267 0.39500663
197 198 199 200 201 202 203
0.38445105 0.97727710 0.43520510 0.16869554 0.17219830 0.05878035 0.21716448
204 205 206 207 208 209 210
0.18384556 0.18192355 0.64682101 0.35975276 0.70106697 1.03414013 0.35317899
211 212 213 214 215 216 217
0.42921059 0.47944086 0.40234009 0.25017393 0.04470913 0.27054309 0.22137404
218 219 220 221 222 223 224
1.18698635 0.50681607 0.11190719 0.11327702 0.28954125 0.33611081 0.74776723
225 226 227 228
0.12225025 0.00000000 0.35218786 0.10231300
> predict(fit,type="risk",se.fit=TRUE)
$fit
1 2 3 4 5 6 7 8
1.2404765 0.6547301 0.5716290 1.2013784 0.5830788 1.2815823 1.5014958 1.6309609
9 10 11 12 13 14 15 16
0.9536592 1.3871959 1.0412167 1.7338000 0.8905418 NA 1.0573939 1.1172899
17 18 19 20 21 22 23 24
1.0518675 1.6384656 1.7465558 0.9957634 0.8805256 0.5372670 0.7264944 0.5622086
25 26 27 28 29 30 31 32
0.7077332 1.2248902 0.6515696 3.7187976 0.9790127 2.1409399 1.2111135 2.1168600
33 34 35 36 37 38 39 40
1.1975053 1.5838000 1.9585669 1.8692030 1.7597526 1.0895532 2.1409399 1.0800127
41 42 43 44 45 46 47 48
1.3916825 1.6682773 0.6442412 1.9365383 0.8481484 1.6437009 0.6831246 1.0953499
49 50 51 52 53 54 55 56
0.7012859 0.8389080 1.2127293 0.6392321 0.6370002 0.5700364 1.0650396 0.5964197
57 58 59 60 61 62 63 64
0.7428929 1.6075072 1.0351212 1.0800127 1.9704933 0.9550493 1.1933106 0.8608298
65 66 67 68 69 70 71 72
1.1722679 2.0519501 1.7142994 0.7345374 0.5808769 1.1659084 0.6192406 0.9244169
73 74 75 76 77 78 79 80
2.5777557 0.9291069 0.9517148 1.2388237 0.5266429 1.0295057 0.7584424 0.6754413
81 82 83 84 85 86 87 88
0.6442412 1.0014125 0.6636416 0.8595988 0.7465360 1.2198707 0.9608649 0.5885092
89 90 91 92 93 94 95 96
0.8842546 2.0930510 0.4753794 1.0518675 1.3301200 1.1150031 0.6455625 2.2180641
97 98 99 100 101 102 103 104
1.1722679 1.2388237 0.8511478 0.6697103 0.5700364 1.1925844 1.0128039 1.1143333
105 106 107 108 109 110 111 112
1.8548748 1.1707060 0.8274120 1.1172899 0.9734900 2.2210235 0.6627251 0.8613550
113 114 115 116 117 118 119 120
1.4463045 0.7086775 1.0064776 2.1899245 1.6551313 2.0004627 1.6950181 1.2765009
121 122 123 124 125 126 127 128
0.6287391 1.5676464 1.1722679 1.6495136 0.6529059 1.0053370 0.6919104 1.1445194
129 130 131 132 133 134 135 136
1.9195968 0.5394778 1.1405143 0.8265359 0.5252493 1.0014125 0.6311224 1.9483304
137 138 139 140 141 142 143 144
1.0694221 1.7929889 0.9645073 0.6709887 2.1125440 1.2404765 1.8786481 1.0922982
145 146 147 148 149 150 151 152
0.7863550 0.8181059 0.5723102 0.8359380 0.7924884 0.5918479 1.1865827 0.5306222
153 154 155 156 157 158 159 160
0.7181071 0.8265359 1.6119444 NA 0.9693810 2.0887670 1.1316466 0.9865752
161 162 163 164 165 166 167 168
0.5567137 0.9625332 1.5944793 1.1143333 0.8106297 0.4884935 0.8134964 0.4966703
169 170 171 172 173 174 175 176
1.0894077 0.6538577 1.0717340 0.6433325 1.1134303 0.4169521 0.9536592 1.0477613
177 178 179 180 181 182 183 184
1.7465558 1.0014125 0.9536592 0.5127358 0.5985803 0.4597993 0.5408413 1.0195796
185 186 187 188 189 190 191 192
0.5580981 0.5214787 2.3621782 0.5847079 1.0651819 0.4902291 1.0246985 1.4467772
193 194 195 196 197 198 199 200
0.9769008 1.0800127 1.0638328 0.6919025 0.9967735 1.0778107 0.5326337 0.9635172
201 202 203 204 205 206 207 208
1.0651819 1.0881732 0.5630464 0.5341218 0.5174645 1.3533017 0.6634419 1.6753425
209 210 211 212 213 214 215 216
1.1405233 0.7389335 0.5830788 1.1659084 1.1271597 1.5068056 0.8608298 0.8608298
217 218 219 220 221 222 223 224
0.8547524 2.1854595 1.6126170 0.6360967 1.0167730 0.9212793 1.2370205 1.2510919
225 226 227 228
0.4720233 2.0211416 1.1527951 0.9183601
$se.fit
1 2 3 4 5 6
0.094027169 0.096340319 0.096185061 0.110144705 0.091221886 0.124003567
7 8 9 10 11 12
0.106470052 0.135893441 0.104263809 0.115204660 0.048057506 0.157626321
13 14 15 16 17 18
0.058398830 NA 0.078593550 0.044525715 0.047523899 0.139753275
19 20 21 22 23 24
0.246130195 0.051683778 0.050651208 0.106747848 0.121191090 0.095563151
25 26 27 28 29 30
0.135232494 0.077970827 0.084316589 0.541641696 0.047411370 0.244541270
31 32 33 34 35 36
0.067316853 0.236761412 0.222247496 0.143779967 0.246770836 0.214866749
37 38 39 40 41 42
0.186808694 0.027994134 0.244541270 0.017746688 0.094899948 0.150429986
43 44 45 46 47 48
0.082038635 0.251128992 0.071539989 0.172653479 0.157627962 0.046664065
49 50 51 52 53 54
0.203630081 0.147427688 0.071868116 0.087051165 0.126710133 0.091078334
55 56 57 58 59 60
0.030346404 0.094111921 0.072518580 0.232795318 0.092391388 0.017746688
61 62 63 64 65 66
0.207337260 0.162712161 0.126511646 0.038549743 0.042876315 0.234595146
67 68 69 70 71 72
0.151669341 0.068462840 0.112880428 0.068678027 0.124246473 0.184637680
73 74 75 76 77 78
0.325442016 0.174862073 0.090441588 0.089040153 0.108376599 0.057550307
79 80 81 82 83 84
0.188633743 0.150191651 0.082038635 0.027564795 0.181878087 0.172125872
85 86 87 88 89 90
0.142365056 0.114741553 0.035859182 0.096819023 0.132484179 0.229864932
91 92 93 94 95 96
0.120689668 0.047523899 0.070339929 0.055381362 0.123547581 0.253870138
97 98 99 100 101 102
0.042876315 0.089040153 0.035190905 0.106227011 0.091078334 0.091298269
103 104 105 106 107 108
0.017787711 0.028641480 0.194430169 0.039989624 0.075782969 0.044525715
109 110 111 112 113 114
0.071209628 0.254965259 0.163546509 0.185211877 0.241649528 0.139074790
115 116 117 118 119 120
0.076796420 0.262556790 0.348185429 0.211911041 0.146845572 0.149423594
121 122 123 124 125 126
0.150969692 0.156065943 0.042876315 0.142648758 0.129688202 0.004890619
127 128 129 130 131 132
0.113985419 0.031310085 0.248637733 0.121183075 0.041502912 0.067248608
133 134 135 136 137 138
0.115359144 0.027564795 0.112511267 0.200585657 0.069255092 0.201817172
139 140 141 142 143 144
0.094786456 0.075667327 0.240338975 0.094027169 0.216098624 0.024974398
145 146 147 148 149 150
0.066191299 0.084423319 0.167625233 0.058808327 0.221289168 0.105873833
151 152 153 154 155 156
0.140449741 0.098993713 0.063583542 0.067248608 0.230942129 NA
157 158 159 160 161 162
0.067558237 0.245408761 0.032338223 0.075589234 0.101745759 0.174851413
163 164 165 166 167 168
0.125897325 0.028641480 0.048065722 0.111659253 0.045260623 0.125085448
169 170 171 172 173 174
0.020095538 0.093808006 0.037378627 0.093118562 0.031761359 0.135544076
175 176 177 178 179 180
0.104263809 0.016586035 0.246130195 0.027564795 0.104263809 0.174088607
181 182 183 184 185 186
0.109727836 0.166211707 0.139230772 0.017941579 0.106388490 0.137198131
187 188 189 190 191 192
0.304795981 0.089505183 0.043311645 0.114439474 0.131445121 0.192173147
193 194 195 196 197 198
0.144436340 0.017746688 0.058484070 0.121193159 0.002956631 0.025613128
199 200 201 202 203 204
0.104623286 0.033429233 0.043311645 0.080773833 0.103942128 0.124008736
205 206 207 208 209 210
0.118294076 0.078206752 0.080505144 0.235804861 0.079727031 0.057789591
211 212 213 214 215 216
0.091221886 0.068678027 0.029421496 0.124248857 0.038549743 0.038549743
217 218 219 220 221 222
0.158976598 0.269332667 0.130275218 0.089792820 0.015369862 0.085131550
223 224 225 226 227 228
0.148494109 0.160862263 0.138362860 0.225740927 0.057778343 0.074788433
> predict(fit,type="terms",se.fit=TRUE)
$fit
age ph.ecog
1 0.130878057 0.03032716
2 0.063011653 -0.54083428
3 -0.072721154 -0.54083428
4 -0.061410086 0.03032716
5 -0.027476885 -0.54083428
6 0.130878057 0.03032716
7 0.063011653 0.60148859
8 0.096944855 0.60148859
9 -0.106654355 0.03032716
10 -0.016165817 0.60148859
11 -0.061410086 0.03032716
12 0.063011653 0.60148859
13 0.063011653 0.03032716
14 NA NA
15 -0.061410086 0.03032716
16 0.051700586 0.03032716
17 0.085633788 0.03032716
18 0.006456317 0.60148859
19 -0.072721154 0.60148859
20 -0.061410086 0.03032716
21 0.051700586 0.03032716
22 -0.151898625 -0.54083428
23 -0.140587557 0.03032716
24 -0.050099019 -0.54083428
25 0.108255923 -0.54083428
26 0.085633788 0.03032716
27 -0.027476885 -0.54083428
28 0.085633788 1.17265002
29 -0.106654355 0.03032716
30 0.130878057 0.60148859
31 0.074322721 0.03032716
32 0.119566990 0.60148859
33 -0.163209692 0.60148859
34 -0.027476885 0.60148859
35 -0.016165817 0.60148859
36 -0.004854750 0.60148859
37 0.029078452 0.60148859
38 0.040389519 0.03032716
39 0.130878057 0.60148859
40 0.017767384 0.03032716
41 0.085633788 0.03032716
42 0.119566990 0.60148859
43 -0.038787952 -0.54083428
44 -0.027476885 0.60148859
45 0.063011653 0.03032716
46 0.153500192 0.60148859
47 0.130878057 -0.54083428
48 0.006456317 0.03032716
49 0.130878057 -0.54083428
50 -0.140587557 0.03032716
51 0.108255923 0.03032716
52 0.006456317 -0.54083428
53 0.063011653 -0.54083428
54 -0.050099019 -0.54083428
55 -0.038787952 0.03032716
56 -0.004854750 -0.54083428
57 0.029078452 -0.54083428
58 -0.061410086 0.60148859
59 -0.050099019 0.03032716
60 0.017767384 0.03032716
61 0.142189124 0.60148859
62 -0.163209692 0.03032716
63 0.119566990 0.03032716
64 0.029078452 0.03032716
65 0.074322721 0.03032716
66 0.063011653 0.60148859
67 0.051700586 0.60148859
68 0.017767384 -0.54083428
69 0.063011653 -0.54083428
70 0.051700586 0.03032716
71 0.006456317 -0.54083428
72 -0.163209692 0.03032716
73 0.130878057 0.60148859
74 -0.253698230 0.03032716
75 -0.106654355 0.03032716
76 0.096944855 0.03032716
77 -0.129276490 -0.54083428
78 -0.072721154 0.03032716
79 0.210055528 -0.54083428
80 0.119566990 -0.54083428
81 -0.038787952 -0.54083428
82 -0.084032221 0.03032716
83 -0.231076095 0.03032716
84 -0.208453961 0.03032716
85 -0.208453961 0.03032716
86 0.096944855 0.03032716
87 -0.004854750 0.03032716
88 -0.016165817 -0.54083428
89 -0.208453961 0.03032716
90 0.108255923 0.60148859
91 0.006456317 -0.54083428
92 0.085633788 0.03032716
93 0.040389519 0.03032716
94 -0.061410086 0.03032716
95 0.074322721 -0.54083428
96 0.108255923 0.60148859
97 0.074322721 0.03032716
98 0.096944855 0.03032716
99 0.017767384 0.03032716
100 0.085633788 -0.54083428
101 -0.050099019 -0.54083428
102 0.074322721 0.03032716
103 -0.072721154 0.03032716
104 0.006456317 0.03032716
105 -0.038787952 0.60148859
106 0.040389519 0.03032716
107 -0.095343288 0.03032716
108 0.051700586 0.03032716
109 -0.084032221 0.03032716
110 0.142189124 0.60148859
111 0.074322721 -0.54083428
112 -0.208453961 0.03032716
113 0.198744461 0.03032716
114 0.142189124 -0.54083428
115 -0.095343288 0.03032716
116 0.153500192 0.60148859
117 -0.151898625 0.60148859
118 0.063011653 0.60148859
119 0.040389519 0.60148859
120 0.198744461 0.03032716
121 0.142189124 -0.54083428
122 -0.027476885 0.60148859
123 0.074322721 0.03032716
124 0.108255923 0.60148859
125 0.085633788 -0.54083428
126 0.040389519 0.03032716
127 -0.140587557 0.03032716
128 0.017767384 0.03032716
129 0.164811259 0.60148859
130 -0.163209692 -0.54083428
131 -0.038787952 0.03032716
132 -0.106654355 0.03032716
133 -0.174520759 -0.54083428
134 -0.084032221 0.03032716
135 0.051700586 -0.54083428
136 0.130878057 0.60148859
137 -0.050099019 0.03032716
138 -0.072721154 0.60148859
139 -0.095343288 0.03032716
140 -0.072721154 -0.54083428
141 0.119566990 0.60148859
142 0.130878057 0.03032716
143 0.153500192 0.60148859
144 0.029078452 0.03032716
145 -0.061410086 0.03032716
146 -0.106654355 0.03032716
147 0.096944855 -0.54083428
148 -0.095343288 0.03032716
149 0.221366595 -0.54083428
150 -0.038787952 -0.54083428
151 0.085633788 0.03032716
152 -0.027476885 -0.54083428
153 -0.004854750 -0.54083428
154 -0.106654355 0.03032716
155 -0.084032221 0.60148859
156 NA NA
157 0.063011653 0.03032716
158 -0.004854750 0.60148859
159 0.006456317 0.03032716
160 -0.072721154 0.03032716
161 -0.004854750 -0.54083428
162 -0.208453961 0.03032716
163 0.074322721 0.60148859
164 0.006456317 0.03032716
165 0.017767384 0.03032716
166 -0.061410086 -0.54083428
167 -0.027476885 0.03032716
168 -0.185831826 -0.54083428
169 -0.016165817 0.03032716
170 0.029078452 -0.54083428
171 -0.016165817 0.03032716
172 -0.050099019 -0.54083428
173 -0.072721154 0.03032716
174 -0.219765028 -0.54083428
175 -0.106654355 0.03032716
176 -0.038787952 0.03032716
177 -0.072721154 0.60148859
178 -0.084032221 0.03032716
179 -0.106654355 0.03032716
180 0.130878057 -0.54083428
181 -0.027476885 -0.54083428
182 -0.265009297 -0.54083428
183 0.040389519 -0.54083428
184 0.029078452 0.03032716
185 -0.129276490 -0.54083428
186 -0.197142894 -0.54083428
187 0.108255923 0.60148859
188 -0.050099019 -0.54083428
189 0.017767384 0.03032716
190 -0.106654355 -0.54083428
191 0.108255923 0.03032716
192 -0.117965423 0.60148859
193 -0.140587557 0.03032716
194 0.017767384 0.03032716
195 0.096944855 0.03032716
196 0.085633788 -0.54083428
197 0.006456317 0.03032716
198 0.017767384 0.03032716
199 -0.117965423 -0.54083428
200 -0.027476885 0.03032716
201 0.017767384 0.03032716
202 0.119566990 0.03032716
203 0.006456317 -0.54083428
204 -0.140587557 -0.54083428
205 0.006456317 -0.54083428
206 -0.004854750 0.60148859
207 -0.084032221 -0.54083428
208 -0.140587557 0.60148859
209 0.074322721 0.03032716
210 -0.038787952 0.03032716
211 -0.027476885 -0.54083428
212 0.051700586 0.03032716
213 0.074322721 0.03032716
214 0.017767384 0.60148859
215 0.029078452 0.03032716
216 0.029078452 0.03032716
217 -0.242387163 0.03032716
218 0.153500192 0.60148859
219 0.085633788 0.60148859
220 -0.061410086 -0.54083428
221 0.051700586 0.03032716
222 0.096944855 0.03032716
223 0.153500192 0.03032716
224 0.164811259 0.03032716
225 -0.265009297 -0.54083428
226 0.142189124 0.60148859
227 0.040389519 0.03032716
228 -0.050099019 0.03032716
$se.fit
age ph.ecog
1 0.119930635 0.007395102
2 0.057740983 0.131879319
3 0.066638322 0.131879319
4 0.056273380 0.007395102
5 0.025178554 0.131879319
6 0.119930635 0.007395102
7 0.057740983 0.146669523
8 0.088835809 0.146669523
9 0.097733148 0.007395102
10 0.014813612 0.146669523
11 0.056273380 0.007395102
12 0.057740983 0.146669523
13 0.057740983 0.007395102
14 NA NA
15 0.056273380 0.007395102
16 0.047376041 0.007395102
17 0.078470867 0.007395102
18 0.005916272 0.146669523
19 0.066638322 0.146669523
20 0.056273380 0.007395102
21 0.047376041 0.007395102
22 0.139192917 0.131879319
23 0.128827975 0.007395102
24 0.045908438 0.131879319
25 0.099200751 0.131879319
26 0.078470867 0.007395102
27 0.025178554 0.131879319
28 0.078470867 0.285943945
29 0.097733148 0.007395102
30 0.119930635 0.146669523
31 0.068105925 0.007395102
32 0.109565693 0.146669523
33 0.149557859 0.146669523
34 0.025178554 0.146669523
35 0.014813612 0.146669523
36 0.004448670 0.146669523
37 0.026646156 0.146669523
38 0.037011098 0.007395102
39 0.119930635 0.146669523
40 0.016281214 0.007395102
41 0.078470867 0.007395102
42 0.109565693 0.146669523
43 0.035543496 0.131879319
44 0.025178554 0.146669523
45 0.057740983 0.007395102
46 0.140660519 0.146669523
47 0.119930635 0.131879319
48 0.005916272 0.007395102
49 0.119930635 0.131879319
50 0.128827975 0.007395102
51 0.099200751 0.007395102
52 0.005916272 0.131879319
53 0.057740983 0.131879319
54 0.045908438 0.131879319
55 0.035543496 0.007395102
56 0.004448670 0.131879319
57 0.026646156 0.131879319
58 0.056273380 0.146669523
59 0.045908438 0.007395102
60 0.016281214 0.007395102
61 0.130295577 0.146669523
62 0.149557859 0.007395102
63 0.109565693 0.007395102
64 0.026646156 0.007395102
65 0.068105925 0.007395102
66 0.057740983 0.146669523
67 0.047376041 0.146669523
68 0.016281214 0.131879319
69 0.057740983 0.131879319
70 0.047376041 0.007395102
71 0.005916272 0.131879319
72 0.149557859 0.007395102
73 0.119930635 0.146669523
74 0.232477395 0.007395102
75 0.097733148 0.007395102
76 0.088835809 0.007395102
77 0.118463033 0.131879319
78 0.066638322 0.007395102
79 0.192485229 0.131879319
80 0.109565693 0.131879319
81 0.035543496 0.131879319
82 0.077003264 0.007395102
83 0.211747511 0.007395102
84 0.191017627 0.007395102
85 0.191017627 0.007395102
86 0.088835809 0.007395102
87 0.004448670 0.007395102
88 0.014813612 0.131879319
89 0.191017627 0.007395102
90 0.099200751 0.146669523
91 0.005916272 0.131879319
92 0.078470867 0.007395102
93 0.037011098 0.007395102
94 0.056273380 0.007395102
95 0.068105925 0.131879319
96 0.099200751 0.146669523
97 0.068105925 0.007395102
98 0.088835809 0.007395102
99 0.016281214 0.007395102
100 0.078470867 0.131879319
101 0.045908438 0.131879319
102 0.068105925 0.007395102
103 0.066638322 0.007395102
104 0.005916272 0.007395102
105 0.035543496 0.146669523
106 0.037011098 0.007395102
107 0.087368206 0.007395102
108 0.047376041 0.007395102
109 0.077003264 0.007395102
110 0.130295577 0.146669523
111 0.068105925 0.131879319
112 0.191017627 0.007395102
113 0.182120287 0.007395102
114 0.130295577 0.131879319
115 0.087368206 0.007395102
116 0.140660519 0.146669523
117 0.139192917 0.146669523
118 0.057740983 0.146669523
119 0.037011098 0.146669523
120 0.182120287 0.007395102
121 0.130295577 0.131879319
122 0.025178554 0.146669523
123 0.068105925 0.007395102
124 0.099200751 0.146669523
125 0.078470867 0.131879319
126 0.037011098 0.007395102
127 0.128827975 0.007395102
128 0.016281214 0.007395102
129 0.151025461 0.146669523
130 0.149557859 0.131879319
131 0.035543496 0.007395102
132 0.097733148 0.007395102
133 0.159922801 0.131879319
134 0.077003264 0.007395102
135 0.047376041 0.131879319
136 0.119930635 0.146669523
137 0.045908438 0.007395102
138 0.066638322 0.146669523
139 0.087368206 0.007395102
140 0.066638322 0.131879319
141 0.109565693 0.146669523
142 0.119930635 0.007395102
143 0.140660519 0.146669523
144 0.026646156 0.007395102
145 0.056273380 0.007395102
146 0.097733148 0.007395102
147 0.088835809 0.131879319
148 0.087368206 0.007395102
149 0.202850171 0.131879319
150 0.035543496 0.131879319
151 0.078470867 0.007395102
152 0.025178554 0.131879319
153 0.004448670 0.131879319
154 0.097733148 0.007395102
155 0.077003264 0.146669523
156 NA NA
157 0.057740983 0.007395102
158 0.004448670 0.146669523
159 0.005916272 0.007395102
160 0.066638322 0.007395102
161 0.004448670 0.131879319
162 0.191017627 0.007395102
163 0.068105925 0.146669523
164 0.005916272 0.007395102
165 0.016281214 0.007395102
166 0.056273380 0.131879319
167 0.025178554 0.007395102
168 0.170287743 0.131879319
169 0.014813612 0.007395102
170 0.026646156 0.131879319
171 0.014813612 0.007395102
172 0.045908438 0.131879319
173 0.066638322 0.007395102
174 0.201382569 0.131879319
175 0.097733148 0.007395102
176 0.035543496 0.007395102
177 0.066638322 0.146669523
178 0.077003264 0.007395102
179 0.097733148 0.007395102
180 0.119930635 0.131879319
181 0.025178554 0.131879319
182 0.242842337 0.131879319
183 0.037011098 0.131879319
184 0.026646156 0.007395102
185 0.118463033 0.131879319
186 0.180652685 0.131879319
187 0.099200751 0.146669523
188 0.045908438 0.131879319
189 0.016281214 0.007395102
190 0.097733148 0.131879319
191 0.099200751 0.007395102
192 0.108098090 0.146669523
193 0.128827975 0.007395102
194 0.016281214 0.007395102
195 0.088835809 0.007395102
196 0.078470867 0.131879319
197 0.005916272 0.007395102
198 0.016281214 0.007395102
199 0.108098090 0.131879319
200 0.025178554 0.007395102
201 0.016281214 0.007395102
202 0.109565693 0.007395102
203 0.005916272 0.131879319
204 0.128827975 0.131879319
205 0.005916272 0.131879319
206 0.004448670 0.146669523
207 0.077003264 0.131879319
208 0.128827975 0.146669523
209 0.068105925 0.007395102
210 0.035543496 0.007395102
211 0.025178554 0.131879319
212 0.047376041 0.007395102
213 0.068105925 0.007395102
214 0.016281214 0.146669523
215 0.026646156 0.007395102
216 0.026646156 0.007395102
217 0.222112453 0.007395102
218 0.140660519 0.146669523
219 0.078470867 0.146669523
220 0.056273380 0.131879319
221 0.047376041 0.007395102
222 0.088835809 0.007395102
223 0.140660519 0.007395102
224 0.151025461 0.007395102
225 0.242842337 0.131879319
226 0.130295577 0.146669523
227 0.037011098 0.007395102
228 0.045908438 0.007395102
>
> # For someone who demands reference='zero'
> pzero <- function(fit)
+ predict(fit, reference="sample") + sum(coef(fit) * fit$means, na.rm=TRUE)
>
>
>
> cleanEx()
> nameEx("predict.survreg")
> ### * predict.survreg
>
> flush(stderr()); flush(stdout())
>
> ### Name: predict.survreg
> ### Title: Predicted Values for a 'survreg' Object
> ### Aliases: predict.survreg predict.survreg.penal
> ### Keywords: survival
>
> ### ** Examples
>
> # Draw figure 1 from Escobar and Meeker, 1992.
> fit <- survreg(Surv(time,status) ~ age + I(age^2), data=stanford2,
+ dist='lognormal')
> with(stanford2, plot(age, time, xlab='Age', ylab='Days',
+ xlim=c(0,65), ylim=c(.1, 10^5), log='y', type='n'))
> with(stanford2, points(age, time, pch=c(2,4)[status+1], cex=.7))
> pred <- predict(fit, newdata=list(age=1:65), type='quantile',
+ p=c(.1, .5, .9))
> matlines(1:65, pred, lty=c(2,1,2), col=1)
>
> # Predicted Weibull survival curve for a lung cancer subject with
> # ECOG score of 2
> lfit <- survreg(Surv(time, status) ~ ph.ecog, data=lung)
> pct <- 1:98/100 # The 100th percentile of predicted survival is at +infinity
> ptime <- predict(lfit, newdata=data.frame(ph.ecog=2), type='quantile',
+ p=pct, se=TRUE)
> matplot(cbind(ptime$fit, ptime$fit + 2*ptime$se.fit,
+ ptime$fit - 2*ptime$se.fit)/30.5, 1-pct,
+ xlab="Months", ylab="Survival", type='l', lty=c(1,2,2), col=1)
>
>
>
> cleanEx()
> nameEx("pseudo")
> ### * pseudo
>
> flush(stderr()); flush(stdout())
>
> ### Name: pseudo
> ### Title: Pseudo values for survival.
> ### Aliases: pseudo
> ### Keywords: survival
>
> ### ** Examples
>
> fit1 <- survfit(Surv(time, status) ~ 1, data=lung)
> yhat <- pseudo(fit1, times=c(365, 730))
> dim(yhat)
[1] 228 2
> lfit <- lm(yhat[,1] ~ ph.ecog + age + sex, data=lung)
>
> # Restricted Mean Time in State (RMST)
> rms <- pseudo(fit1, times= 730, type='RMST') # 2 years
> rfit <- lm(rms ~ ph.ecog + sex, data=lung)
> rhat <- predict(rfit, newdata=expand.grid(ph.ecog=0:3, sex=1:2), se.fit=TRUE)
> # print it out nicely
> temp1 <- cbind(matrix(rhat$fit, 4,2))
> temp2 <- cbind(matrix(rhat$se.fit, 4, 2))
> temp3 <- cbind(temp1[,1], temp2[,1], temp1[,2], temp2[,2])
> dimnames(temp3) <- list(paste("ph.ecog", 0:3),
+ c("Male RMST", "(se)", "Female RMST", "(se)"))
>
> round(temp3, 1)
Male RMST (se) Female RMST (se)
ph.ecog 0 393.7 28.6 510.6 31.5
ph.ecog 1 307.8 19.8 424.7 24.4
ph.ecog 2 221.9 29.7 338.8 33.4
ph.ecog 3 136.1 47.9 253.0 50.5
> # compare this to the fully non-parametric estimate
> fit2 <- survfit(Surv(time, status) ~ ph.ecog, data=lung)
> print(fit2, rmean=730)
Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung)
1 observation deleted due to missingness
n events rmean* se(rmean) median 0.95LCL 0.95UCL
ph.ecog=0 63 37 429 32.2 394 348 574
ph.ecog=1 113 82 366 22.4 306 268 429
ph.ecog=2 50 44 256 30.0 199 156 288
ph.ecog=3 1 1 118 0.0 118 NA NA
* restricted mean with upper limit = 730
> # the estimate for ph.ecog=3 is very unstable (n=1), pseudovalues smooth it.
> #
> # In all the above we should be using the robust variance, e.g., svyglm, but
> # a recommended package can't depend on external libraries.
> # See the vignette for a more complete exposition.
>
>
>
> cleanEx()
> nameEx("pspline")
> ### * pspline
>
> flush(stderr()); flush(stdout())
>
> ### Name: pspline
> ### Title: Smoothing splines using a pspline basis
> ### Aliases: pspline psplineinverse
> ### Keywords: survival
>
> ### ** Examples
>
> lfit6 <- survreg(Surv(time, status)~pspline(age, df=2), lung)
> plot(lung$age, predict(lfit6), xlab='Age', ylab="Spline prediction")
> title("Cancer Data")
> fit0 <- coxph(Surv(time, status) ~ ph.ecog + age, lung)
> fit1 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), lung)
> fit3 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,8), lung)
> fit0
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + age, data = lung)
coef exp(coef) se(coef) z p
ph.ecog 0.443485 1.558128 0.115831 3.829 0.000129
age 0.011281 1.011345 0.009319 1.211 0.226082
Likelihood ratio test=19.06 on 2 df, p=7.279e-05
n= 227, number of events= 164
(1 observation deleted due to missingness)
> fit1
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3),
data = lung)
coef se(coef) se2 Chisq DF p
ph.ecog 0.44802 0.11707 0.11678 14.64453 1.00 0.00013
pspline(age, 3), linear 0.01126 0.00928 0.00928 1.47231 1.00 0.22498
pspline(age, 3), nonlin 2.07924 2.08 0.37143
Iterations: 4 outer, 12 Newton-Raphson
Theta= 0.861
Degrees of freedom for terms= 1.0 3.1
Likelihood ratio test=21.9 on 4.08 df, p=2e-04
n= 227, number of events= 164
(1 observation deleted due to missingness)
> fit3
Call:
coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 8),
data = lung)
coef se(coef) se2 Chisq DF p
ph.ecog 0.47640 0.12024 0.11925 15.69732 1.00 7.4e-05
pspline(age, 8), linear 0.01172 0.00923 0.00923 1.61161 1.00 0.20
pspline(age, 8), nonlin 6.93188 6.99 0.43
Iterations: 5 outer, 15 Newton-Raphson
Theta= 0.691
Degrees of freedom for terms= 1 8
Likelihood ratio test=27.6 on 8.97 df, p=0.001
n= 227, number of events= 164
(1 observation deleted due to missingness)
>
>
>
> cleanEx()
> nameEx("pyears")
> ### * pyears
>
> flush(stderr()); flush(stdout())
>
> ### Name: pyears
> ### Title: Person Years
> ### Aliases: pyears
> ### Keywords: survival
>
> ### ** Examples
>
> # Look at progression rates jointly by calendar date and age
> #
> temp.yr <- tcut(mgus$dxyr, 55:92, labels=as.character(55:91))
> temp.age <- tcut(mgus$age, 34:101, labels=as.character(34:100))
> ptime <- ifelse(is.na(mgus$pctime), mgus$futime, mgus$pctime)
> pstat <- ifelse(is.na(mgus$pctime), 0, 1)
> pfit <- pyears(Surv(ptime/365.25, pstat) ~ temp.yr + temp.age + sex, mgus,
+ data.frame=TRUE)
> # Turn the factor back into numerics for regression
> tdata <- pfit$data
> tdata$age <- as.numeric(as.character(tdata$temp.age))
> tdata$year<- as.numeric(as.character(tdata$temp.yr))
> fit1 <- glm(event ~ year + age+ sex +offset(log(pyears)),
+ data=tdata, family=poisson)
> ## Not run:
> ##D # fit a gam model
> ##D gfit.m <- gam(y ~ s(age) + s(year) + offset(log(time)),
> ##D family = poisson, data = tdata)
> ## End(Not run)
>
> # Example #2 Create the hearta data frame:
> hearta <- by(heart, heart$id,
+ function(x)x[x$stop == max(x$stop),])
> hearta <- do.call("rbind", hearta)
> # Produce pyears table of death rates on the surgical arm
> # The first is by age at randomization, the second by current age
> fit1 <- pyears(Surv(stop/365.25, event) ~ cut(age + 48, c(0,50,60,70,100)) +
+ surgery, data = hearta, scale = 1)
> fit2 <- pyears(Surv(stop/365.25, event) ~ tcut(age + 48, c(0,50,60,70,100)) +
+ surgery, data = hearta, scale = 1)
> fit1$event/fit1$pyears #death rates on the surgery and non-surg arm
surgery
cut(age + 48, c(0, 50, 60, 70, 100)) 0 1
(0,50] 0.7615378 0.3036881
(50,60] 2.0068681 0.9979508
(60,70] 5.1083916 NaN
(70,100] NaN NaN
>
> fit2$event/fit2$pyears #death rates on the surgery and non-surg arm
surgery
tcut(age + 48, c(0, 50, 60, 70, 100)) 0 1
0+ thru 50 0.8013285 0.2636994
50+ thru 60 1.6119238 0.6564817
60+ thru 70 3.9701087 NaN
70+ thru 100 NaN NaN
>
>
>
> cleanEx()
> nameEx("quantile.survfit")
> ### * quantile.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: quantile.survfit
> ### Title: Quantiles from a survfit object
> ### Aliases: quantile.survfit quantile.survfitms median.survfit
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- survfit(Surv(time, status) ~ ph.ecog, data=lung)
> quantile(fit)
$quantile
25 50 75
ph.ecog=0 285 394 655
ph.ecog=1 181 306 550
ph.ecog=2 105 199 351
ph.ecog=3 118 118 118
$lower
25 50 75
ph.ecog=0 189 348 558
ph.ecog=1 156 268 460
ph.ecog=2 61 156 285
ph.ecog=3 NA NA NA
$upper
25 50 75
ph.ecog=0 350 574 NA
ph.ecog=1 223 429 689
ph.ecog=2 163 288 654
ph.ecog=3 NA NA NA
>
> cfit <- coxph(Surv(time, status) ~ age + strata(ph.ecog), data=lung)
> csurv<- survfit(cfit, newdata=data.frame(age=c(40, 60, 80)),
+ conf.type ="none")
> temp <- quantile(csurv, 1:5/10)
> temp[2,3,] # quantiles for second level of ph.ecog, age=80
10 20 30 40 50
92 144 181 218 270
> quantile(csurv[2,3], 1:5/10) # quantiles of a single curve, same result
10 20 30 40 50
92 144 181 218 270
>
>
>
> cleanEx()
> nameEx("reliability")
> ### * reliability
>
> flush(stderr()); flush(stdout())
>
> ### Name: reliability
> ### Title: Reliability data sets
> ### Aliases: reliability capacitor cracks genfan ifluid imotor turbine
> ### valveSeat
> ### Keywords: datasets
>
> ### ** Examples
>
> survreg(Surv(time, status) ~ temperature + voltage, capacitor)
Call:
survreg(formula = Surv(time, status) ~ temperature + voltage,
data = capacitor)
Coefficients:
(Intercept) temperature voltage
13.40701688 -0.02890466 -0.00591082
Scale= 0.3638092
Loglik(model)= -244.2 Loglik(intercept only)= -254.5
Chisq= 20.57 on 2 degrees of freedom, p= 3.41e-05
n= 64
>
>
>
> cleanEx()
> nameEx("residuals.coxph")
> ### * residuals.coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: residuals.coxph
> ### Title: Calculate Residuals for a 'coxph' Fit
> ### Aliases: residuals.coxph.penal residuals.coxph.null residuals.coxph
> ### residuals.coxphms
> ### Keywords: survival
>
> ### ** Examples
>
>
> fit <- coxph(Surv(start, stop, event) ~ (age + surgery)* transplant,
+ data=heart)
> mresid <- resid(fit, collapse=heart$id)
>
>
>
> cleanEx()
> nameEx("residuals.survfit")
> ### * residuals.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: residuals.survfit
> ### Title: IJ residuals from a survfit object.
> ### Aliases: residuals.survfit
>
> ### ** Examples
>
> fit <- survfit(Surv(time, status) ~ x, aml)
> resid(fit, times=c(24, 48), type="RMTS")
[,1] [,2]
[1,] -1.0836777 -2.076652893
[2,] -0.7200413 -1.713016529
[3,] 0.2004132 0.421074380
[4,] -0.3237345 -1.468414256
[5,] 0.1876291 -0.957050620
[6,] 0.2899019 0.965676653
[7,] 0.2899019 -0.359777893
[8,] 0.2899019 0.008403926
[9,] 0.2899019 1.726585744
[10,] 0.2899019 1.726585744
[11,] 0.2899019 1.726585744
[12,] -1.0057870 -1.475694444
[13,] -1.0057870 -1.475694444
[14,] -0.7557870 -1.225694444
[15,] -0.7557870 -1.225694444
[16,] -0.4224537 -0.892361111
[17,] 0.5636574 0.899305556
[18,] 0.4826389 -0.121527778
[19,] 0.5798611 0.267361111
[20,] 0.5798611 0.559027778
[21,] 0.5798611 0.850694444
[22,] 0.5798611 1.822916667
[23,] 0.5798611 2.017361111
>
>
>
> cleanEx()
> nameEx("residuals.survreg")
> ### * residuals.survreg
>
> flush(stderr()); flush(stdout())
>
> ### Name: residuals.survreg
> ### Title: Compute Residuals for 'survreg' Objects
> ### Aliases: residuals.survreg residuals.survreg.penal
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- survreg(Surv(futime, death) ~ age + sex, mgus2)
> summary(fit) # age and sex are both important
Call:
survreg(formula = Surv(futime, death) ~ age + sex, data = mgus2)
Value Std. Error z p
(Intercept) 8.85979 0.23842 37.16 < 2e-16
age -0.05360 0.00312 -17.19 < 2e-16
sexM -0.31874 0.06357 -5.01 5.3e-07
Log(scale) -0.02840 0.02787 -1.02 0.31
Scale= 0.972
Weibull distribution
Loglik(model)= -5528.3 Loglik(intercept only)= -5699
Chisq= 341.42 on 2 degrees of freedom, p= 7.3e-75
Number of Newton-Raphson Iterations: 5
n= 1384
>
> rr <- residuals(fit, type='matrix')
> sum(rr[,1]) - with(mgus2, sum(log(futime[death==1]))) # loglik
[1] -5528.267
>
> plot(mgus2$age, rr[,2], col= (1+mgus2$death)) # ldresp
>
>
>
> cleanEx()
> nameEx("retinopathy")
> ### * retinopathy
>
> flush(stderr()); flush(stdout())
>
> ### Name: retinopathy
> ### Title: Diabetic Retinopathy
> ### Aliases: retinopathy
> ### Keywords: datasets
>
> ### ** Examples
>
> coxph(Surv(futime, status) ~ type + trt, cluster= id, retinopathy)
Call:
coxph(formula = Surv(futime, status) ~ type + trt, data = retinopathy,
cluster = id)
coef exp(coef) se(coef) robust se z p
typeadult 0.05388 1.05536 0.16211 0.17864 0.302 0.763
trt -0.77893 0.45890 0.16893 0.14851 -5.245 1.56e-07
Likelihood ratio test=22.48 on 2 df, p=1.312e-05
n= 394, number of events= 155
>
>
>
> cleanEx()
> nameEx("rhDNase")
> ### * rhDNase
>
> flush(stderr()); flush(stdout())
>
> ### Name: rhDNase
> ### Title: rhDNASE data set
> ### Aliases: rhDNase
> ### Keywords: datasets
>
> ### ** Examples
>
> # Build the start-stop data set for analysis, and
> # replicate line 2 of table 8.13 in the book
> first <- subset(rhDNase, !duplicated(id)) #first row for each subject
> dnase <- tmerge(first, first, id=id, tstop=as.numeric(end.dt -entry.dt))
>
> # Subjects whose fu ended during the 6 day window are the reason for
> # this next line
> temp.end <- with(rhDNase, pmin(ivstop+6, end.dt-entry.dt))
> dnase <- tmerge(dnase, rhDNase, id=id,
+ infect=event(ivstart),
+ end= event(temp.end))
> # toss out the non-at-risk intervals, and extra variables
> # 3 subjects had an event on their last day of fu, infect=1 and end=1
> dnase <- subset(dnase, (infect==1 | end==0), c(id:trt, fev:infect))
> agfit <- coxph(Surv(tstart, tstop, infect) ~ trt + fev, cluster=id,
+ data=dnase)
>
>
>
> cleanEx()
> nameEx("ridge")
> ### * ridge
>
> flush(stderr()); flush(stdout())
>
> ### Name: ridge
> ### Title: Ridge regression
> ### Aliases: ridge
> ### Keywords: survival
>
> ### ** Examples
>
>
> coxph(Surv(futime, fustat) ~ rx + ridge(age, ecog.ps, theta=1),
+ ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ rx + ridge(age, ecog.ps,
theta = 1), data = ovarian)
coef se(coef) se2 Chisq DF p
rx -0.8564 0.6161 0.6156 1.9323 1 0.1645
ridge(age) 0.1229 0.0385 0.0354 10.2127 1 0.0014
ridge(ecog.ps) 0.1093 0.5734 0.5484 0.0363 1 0.8488
Iterations: 1 outer, 5 Newton-Raphson
Degrees of freedom for terms= 1.0 1.8
Likelihood ratio test=15.6 on 2.76 df, p=0.001
n= 26, number of events= 12
>
> lfit0 <- survreg(Surv(time, status) ~1, lung)
> lfit1 <- survreg(Surv(time, status) ~ age + ridge(ph.ecog, theta=5), lung)
> lfit2 <- survreg(Surv(time, status) ~ sex + ridge(age, ph.ecog, theta=1), lung)
> lfit3 <- survreg(Surv(time, status) ~ sex + age + ph.ecog, lung)
>
>
>
>
> cleanEx()
> nameEx("rotterdam")
> ### * rotterdam
>
> flush(stderr()); flush(stdout())
>
> ### Name: rotterdam
> ### Title: Breast cancer data set used in Royston and Altman (2013)
> ### Aliases: rotterdam
> ### Keywords: datasets survival
>
> ### ** Examples
>
> # liberal definition of rfs (count later deaths)
> rfs <- pmax(rotterdam$recur, rotterdam$death)
> rfstime <- with(rotterdam, ifelse(recur==1, rtime, dtime))
> fit1 <- coxph(Surv(rfstime, rfs) ~ pspline(age) + meno + size +
+ pspline(nodes) + er, data = rotterdam)
>
> # conservative (no deaths after last fu for recurrence)
> ignore <- with(rotterdam, recur ==0 & death==1 & rtime < dtime)
> table(ignore)
ignore
FALSE TRUE
2939 43
> rfs2 <- with(rotterdam, ifelse(recur==1 | ignore, recur, death))
> rfstime2 <- with(rotterdam, ifelse(recur==1 | ignore, rtime, dtime))
> fit2 <- coxph(Surv(rfstime2, rfs2) ~ pspline(age) + meno + size +
+ pspline(nodes) + er, data = rotterdam)
>
> # Note: Both age and nodes show non-linear effects.
> # Royston and Altman used fractional polynomials for the nonlinear terms
>
>
>
> cleanEx()
> nameEx("royston")
> ### * royston
>
> flush(stderr()); flush(stdout())
>
> ### Name: royston
> ### Title: Compute Royston's D for a Cox model
> ### Aliases: royston
> ### Keywords: survival
>
> ### ** Examples
>
> # An example used in Royston and Sauerbrei
> pbc2 <- na.omit(pbc) # no missing values
> cfit <- coxph(Surv(time, status==2) ~ age + log(bili) + edema + albumin +
+ stage + copper, data=pbc2, ties="breslow")
> royston(cfit)
D se(D) R.D R.KO R.N C.GH
2.6917766 0.2273352 0.6336693 0.5554885 0.4714442 0.7735923
>
>
>
> cleanEx()
> nameEx("rttright")
> ### * rttright
>
> flush(stderr()); flush(stdout())
>
> ### Name: rttright
> ### Title: Compute redistribute-to-the-right weights
> ### Aliases: rttright
> ### Keywords: survival
>
> ### ** Examples
>
> afit <- survfit(Surv(time, status) ~1, data=aml)
> rwt <- rttright(Surv(time, status) ~1, data=aml)
>
> # Reproduce a Kaplan-Meier
> index <- order(aml$time)
> cdf <- cumsum(rwt[index]) # weighted CDF
> cdf <- cdf[!duplicated(aml$time[index], fromLast=TRUE)] # remove duplicate times
> cbind(time=afit$time, KM= afit$surv, RTTR= 1-cdf)
time KM RTTR
[1,] 5 0.91304348 0.91304348
[2,] 8 0.82608696 0.82608696
[3,] 9 0.78260870 0.78260870
[4,] 12 0.73913043 0.73913043
[5,] 13 0.69565217 0.69565217
[6,] 16 0.69565217 0.69565217
[7,] 18 0.64596273 0.64596273
[8,] 23 0.54658385 0.54658385
[9,] 27 0.49689441 0.49689441
[10,] 28 0.49689441 0.49689441
[11,] 30 0.44168392 0.44168392
[12,] 31 0.38647343 0.38647343
[13,] 33 0.33126294 0.33126294
[14,] 34 0.27605245 0.27605245
[15,] 43 0.22084196 0.22084196
[16,] 45 0.16563147 0.16563147
[17,] 48 0.08281573 0.08281573
[18,] 161 0.08281573 0.08281573
>
> # Hormonal patients have a diffent censoring pattern
> wt2 <- rttright(Surv(dtime, death) ~ hormon, rotterdam, times= 365*c(3, 5))
> dim(wt2)
[1] 2982 2
>
>
>
> cleanEx()
> nameEx("solder")
> ### * solder
>
> flush(stderr()); flush(stdout())
>
> ### Name: solder
> ### Title: Data from a soldering experiment
> ### Aliases: solder
> ### Keywords: datasets
>
> ### ** Examples
>
> # The balanced subset used by Chambers and Hastie
> # contains the first 180 of each mask and deletes mask A6.
> index <- 1 + (1:nrow(solder)) - match(solder$Mask, solder$Mask)
> solder.balance <- droplevels(subset(solder, Mask != "A6" & index <= 180))
>
>
>
> cleanEx()
> nameEx("statefig")
> ### * statefig
>
> flush(stderr()); flush(stdout())
>
> ### Name: statefig
> ### Title: Draw a state space figure.
> ### Aliases: statefig
> ### Keywords: survival hplot
>
> ### ** Examples
>
> # Draw a simple competing risks figure
> states <- c("Entry", "Complete response", "Relapse", "Death")
> connect <- matrix(0, 4, 4, dimnames=list(states, states))
> connect[1, -1] <- c(1.1, 1, 0.9)
> statefig(c(1, 3), connect)
>
>
>
> cleanEx()
> nameEx("strata")
> ### * strata
>
> flush(stderr()); flush(stdout())
>
> ### Name: strata
> ### Title: Identify Stratification Variables
> ### Aliases: strata
> ### Keywords: survival
>
> ### ** Examples
>
> a <- factor(rep(1:3,4), labels=c("low", "medium", "high"))
> b <- factor(rep(1:4,3))
> levels(strata(b))
[1] "1" "2" "3" "4"
> levels(strata(a,b,shortlabel=TRUE))
[1] "low, 1" "low, 2" "low, 3" "low, 4" "medium, 1" "medium, 2"
[7] "medium, 3" "medium, 4" "high, 1" "high, 2" "high, 3" "high, 4"
>
> coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ age + strata(rx), data = ovarian)
coef exp(coef) se(coef) z p
age 0.13735 1.14723 0.04741 2.897 0.00376
Likelihood ratio test=12.69 on 1 df, p=0.0003678
n= 26, number of events= 12
>
>
>
> cleanEx()
> nameEx("summary.aareg")
> ### * summary.aareg
>
> flush(stderr()); flush(stdout())
>
> ### Name: summary.aareg
> ### Title: Summarize an aareg fit
> ### Aliases: summary.aareg
> ### Keywords: survival
>
> ### ** Examples
>
> afit <- aareg(Surv(time, status) ~ age + sex + ph.ecog, data=lung,
+ dfbeta=TRUE)
> summary(afit)
slope coef se(coef) robust se z p
Intercept 5.05e-03 5.87e-03 4.74e-03 0.00477 1.23 0.219000
age 4.01e-05 7.15e-05 7.23e-05 0.00007 1.02 0.307000
sex -3.16e-03 -4.03e-03 1.22e-03 0.00123 -3.28 0.001030
ph.ecog 3.01e-03 3.67e-03 1.02e-03 0.00102 3.62 0.000299
Chisq=22.84 on 3 df, p=4.36e-05; test weights=aalen
> ## Not run:
> ##D slope test se(test) robust se z p
> ##D Intercept 5.05e-03 1.9 1.54 1.55 1.23 0.219000
> ##D age 4.01e-05 108.0 109.00 106.00 1.02 0.307000
> ##D sex -3.16e-03 -19.5 5.90 5.95 -3.28 0.001030
> ##D ph.ecog 3.01e-03 33.2 9.18 9.17 3.62 0.000299
> ##D
> ##D Chisq=22.84 on 3 df, p=4.4e-05; test weights=aalen
> ## End(Not run)
>
> summary(afit, maxtime=600)
slope coef se(coef) robust se z p
Intercept 4.16e-03 6.67e-03 4.62e-03 0.004580 1.450 0.146000
age 2.82e-05 5.74e-05 7.07e-05 0.000067 0.857 0.392000
sex -2.54e-03 -4.30e-03 1.17e-03 0.001180 -3.660 0.000256
ph.ecog 2.47e-03 3.54e-03 9.99e-04 0.000972 3.640 0.000271
Chisq=27.08 on 3 df, p=5.66e-06; test weights=aalen
> ## Not run:
> ##D slope test se(test) robust se z p
> ##D Intercept 4.16e-03 2.13 1.48 1.47 1.450 0.146000
> ##D age 2.82e-05 85.80 106.00 100.00 0.857 0.392000
> ##D sex -2.54e-03 -20.60 5.61 5.63 -3.660 0.000256
> ##D ph.ecog 2.47e-03 31.60 8.91 8.67 3.640 0.000271
> ##D
> ##D Chisq=27.08 on 3 df, p=5.7e-06; test weights=aalen
> ## End(Not run)
>
>
> cleanEx()
> nameEx("summary.coxph")
> ### * summary.coxph
>
> flush(stderr()); flush(stdout())
>
> ### Name: summary.coxph
> ### Title: Summary method for Cox models
> ### Aliases: summary.coxph
> ### Keywords: survival
>
> ### ** Examples
>
> fit <- coxph(Surv(time, status) ~ age + sex, lung)
> summary(fit)
Call:
coxph(formula = Surv(time, status) ~ age + sex, data = lung)
n= 228, number of events= 165
coef exp(coef) se(coef) z Pr(>|z|)
age 0.017045 1.017191 0.009223 1.848 0.06459 .
sex -0.513219 0.598566 0.167458 -3.065 0.00218 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0172 0.9831 0.9990 1.0357
sex 0.5986 1.6707 0.4311 0.8311
Concordance= 0.603 (se = 0.025 )
Likelihood ratio test= 14.12 on 2 df, p=9e-04
Wald test = 13.47 on 2 df, p=0.001
Score (logrank) test = 13.72 on 2 df, p=0.001
>
>
>
> cleanEx()
> nameEx("summary.survfit")
> ### * summary.survfit
>
> flush(stderr()); flush(stdout())
>
> ### Name: summary.survfit
> ### Title: Summary of a Survival Curve
> ### Aliases: summary.survfit
> ### Keywords: survival
>
> ### ** Examples
>
> summary( survfit( Surv(futime, fustat)~1, data=ovarian))
Call: survfit(formula = Surv(futime, fustat) ~ 1, data = ovarian)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
59 26 1 0.962 0.0377 0.890 1.000
115 25 1 0.923 0.0523 0.826 1.000
156 24 1 0.885 0.0627 0.770 1.000
268 23 1 0.846 0.0708 0.718 0.997
329 22 1 0.808 0.0773 0.670 0.974
353 21 1 0.769 0.0826 0.623 0.949
365 20 1 0.731 0.0870 0.579 0.923
431 17 1 0.688 0.0919 0.529 0.894
464 15 1 0.642 0.0965 0.478 0.862
475 14 1 0.596 0.0999 0.429 0.828
563 12 1 0.546 0.1032 0.377 0.791
638 11 1 0.497 0.1051 0.328 0.752
> summary( survfit( Surv(futime, fustat)~rx, data=ovarian))
Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
rx=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
59 13 1 0.923 0.0739 0.789 1.000
115 12 1 0.846 0.1001 0.671 1.000
156 11 1 0.769 0.1169 0.571 1.000
268 10 1 0.692 0.1280 0.482 0.995
329 9 1 0.615 0.1349 0.400 0.946
431 8 1 0.538 0.1383 0.326 0.891
638 5 1 0.431 0.1467 0.221 0.840
rx=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
353 13 1 0.923 0.0739 0.789 1.000
365 12 1 0.846 0.1001 0.671 1.000
464 9 1 0.752 0.1256 0.542 1.000
475 8 1 0.658 0.1407 0.433 1.000
563 7 1 0.564 0.1488 0.336 0.946
>
>
>
> cleanEx()
> nameEx("survSplit")
> ### * survSplit
>
> flush(stderr()); flush(stdout())
>
> ### Name: survSplit
> ### Title: Split a survival data set at specified times
> ### Aliases: survSplit
> ### Keywords: survival utilities
>
> ### ** Examples
>
> fit1 <- coxph(Surv(time, status) ~ karno + age + trt, veteran)
> plot(cox.zph(fit1)[1])
> # a cox.zph plot of the data suggests that the effect of Karnofsky score
> # begins to diminish by 60 days and has faded away by 120 days.
> # Fit a model with separate coefficients for the three intervals.
> #
> vet2 <- survSplit(Surv(time, status) ~., veteran,
+ cut=c(60, 120), episode ="timegroup")
> fit2 <- coxph(Surv(tstart, time, status) ~ karno* strata(timegroup) +
+ age + trt, data= vet2)
> c(overall= coef(fit1)[1],
+ t0_60 = coef(fit2)[1],
+ t60_120= sum(coef(fit2)[c(1,4)]),
+ t120 = sum(coef(fit2)[c(1,5)]))
overall.karno t0_60.karno t60_120 t120
-0.034443897 -0.049176157 -0.011031558 -0.007629841
>
> # Sometimes we want to split on one scale and analyse on another
> # Add a "current age" variable to the mgus2 data set.
> temp1 <- mgus2
> temp1$endage <- mgus2$age + mgus2$futime/12 # futime is in months
> temp1$startage <- temp1$age
> temp2 <- survSplit(Surv(age, endage, death) ~ ., temp1, cut=25:100,
+ start= "age1", end= "age2")
>
> # restore the time since enrollment scale
> temp2$time1 <- (temp2$age1 - temp2$startage)*12
> temp2$time2 <- (temp2$age2 - temp2$startage)*12
>
> # In this data set, initial age and current age have similar utility
> mfit1 <- coxph(Surv(futime, death) ~ age + sex, data=mgus2)
> mfit2 <- coxph(Surv(time1, time2, death) ~ age1 + sex, data=temp2)
>
>
>
> cleanEx()
> nameEx("survcondense")
> ### * survcondense
>
> flush(stderr()); flush(stdout())
>
> ### Name: survcondense
> ### Title: Shorten a (time1, time2) survival dataset
> ### Aliases: survcondense
> ### Keywords: survival
>
> ### ** Examples
>
> dim(aml)
[1] 23 3
> test1 <- survSplit(Surv(time, status) ~ ., data=aml,
+ cut=c(10, 20, 30), id="newid")
> dim(test1)
[1] 62 5
>
> # remove the added rows
> test2 <- survcondense(Surv(tstart, time, status) ~ x, test1, id=newid)
> dim(test2)
[1] 23 5
>
>
>
> cleanEx()
> nameEx("survdiff")
> ### * survdiff
>
> flush(stderr()); flush(stdout())
>
> ### Name: survdiff
> ### Title: Test Survival Curve Differences
> ### Aliases: survdiff print.survdiff
> ### Keywords: survival
>
> ### ** Examples
>
> ## Two-sample test
> survdiff(Surv(futime, fustat) ~ rx,data=ovarian)
Call:
survdiff(formula = Surv(futime, fustat) ~ rx, data = ovarian)
N Observed Expected (O-E)^2/E (O-E)^2/V
rx=1 13 7 5.23 0.596 1.06
rx=2 13 5 6.77 0.461 1.06
Chisq= 1.1 on 1 degrees of freedom, p= 0.3
>
> ## Stratified 7-sample test
>
> survdiff(Surv(time, status) ~ pat.karno + strata(inst), data=lung)
Call:
survdiff(formula = Surv(time, status) ~ pat.karno + strata(inst),
data = lung)
n=224, 4 observations deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
pat.karno=30 2 1 0.692 0.13720 0.15752
pat.karno=40 2 1 1.099 0.00889 0.00973
pat.karno=50 4 4 1.166 6.88314 7.45359
pat.karno=60 30 27 16.298 7.02790 9.57333
pat.karno=70 41 31 26.358 0.81742 1.14774
pat.karno=80 50 38 41.938 0.36978 0.60032
pat.karno=90 60 38 47.242 1.80800 3.23078
pat.karno=100 35 21 26.207 1.03451 1.44067
Chisq= 21.4 on 7 degrees of freedom, p= 0.003
>
> ## Expected survival for heart transplant patients based on
> ## US mortality tables
> expect <- survexp(futime ~ 1, data=jasa, cohort=FALSE,
+ rmap= list(age=(accept.dt - birth.dt), sex=1, year=accept.dt),
+ ratetable=survexp.us)
> ## actual survival is much worse (no surprise)
> survdiff(Surv(jasa$futime, jasa$fustat) ~ offset(expect))
Call:
survdiff(formula = Surv(jasa$futime, jasa$fustat) ~ offset(expect))
Observed Expected Z p
75.000 0.644 -92.681 0.000
>
> # The free light chain data set is close to the population.
> e2 <- survexp(futime ~ 1, data=flchain, cohort=FALSE,
+ rmap= list(age= age*365.25, sex=sex,
+ year=as.Date(paste0(sample.yr, "-07-01"))),
+ ratetable= survexp.mn)
> survdiff(Surv(futime, death) ~ offset(e2), flchain)
Call:
survdiff(formula = Surv(futime, death) ~ offset(e2), data = flchain)
Observed Expected Z p
2169.0000 2073.9133 -2.0880 0.0368
>
>
>
> cleanEx()
> nameEx("survexp")
> ### * survexp
>
> flush(stderr()); flush(stdout())
>
> ### Name: survexp
> ### Title: Compute Expected Survival
> ### Aliases: survexp print.survexp
> ### Keywords: survival
>
> ### ** Examples
>
> #
> # Stanford heart transplant data
> # We don't have sex in the data set, but know it to be nearly all males.
> # Estimate of conditional survival
> fit1 <- survexp(futime ~ 1, rmap=list(sex="male", year=accept.dt,
+ age=(accept.dt-birth.dt)), method='conditional', data=jasa)
> summary(fit1, times=1:10*182.5, scale=365) #expected survival by 1/2 years
Call: survexp(formula = futime ~ 1, data = jasa, rmap = list(sex = "male",
year = accept.dt, age = (accept.dt - birth.dt)), method = "conditional")
time n.risk survival
0.5 41 0.996
1.0 28 0.993
1.5 21 0.989
2.0 16 0.986
2.5 13 0.983
3.0 8 0.980
3.5 7 0.977
4.0 3 0.972
4.5 1 0.969
>
> # Estimate of expected survival stratified by prior surgery
> survexp(~ surgery, rmap= list(sex="male", year=accept.dt,
+ age=(accept.dt-birth.dt)), method='ederer', data=jasa,
+ times=1:10 * 182.5)
Call:
survexp(formula = ~surgery, data = jasa, rmap = list(sex = "male",
year = accept.dt, age = (accept.dt - birth.dt)), times = 1:10 *
182.5, method = "ederer")
age ranges from 8.8 to 64.4 years
male: 103 female: 0
date of entry from 1967-09-13 to 1974-03-22
time nrisk1 nrisk2 surgery=0 surgery=1
182 87 16 0.996 0.996
365 87 16 0.991 0.993
548 87 16 0.987 0.989
730 87 16 0.982 0.985
912 87 16 0.978 0.981
1095 87 16 0.973 0.977
1278 87 16 0.968 0.973
1460 87 16 0.963 0.969
1642 87 16 0.958 0.964
1825 87 16 0.952 0.960
>
> ## Compare the survival curves for the Mayo PBC data to Cox model fit
> ##
> pfit <-coxph(Surv(time,status>0) ~ trt + log(bili) + log(protime) + age +
+ platelet, data=pbc)
> plot(survfit(Surv(time, status>0) ~ trt, data=pbc), mark.time=FALSE)
> lines(survexp( ~ trt, ratetable=pfit, data=pbc), col='purple')
>
>
>
> cleanEx()
> nameEx("survexp.us")
> ### * survexp.us
>
> flush(stderr()); flush(stdout())
>
> ### Name: ratetables
> ### Title: Census Data Sets for the Expected Survival and Person Years
> ### Functions
> ### Aliases: survexp.us survexp.usr survexp.mn
> ### Keywords: survival datasets
>
> ### ** Examples
>
> survexp.uswhite <- survexp.usr[,,"white",]
>
>
>
> cleanEx()
> nameEx("survfit.formula")
> ### * survfit.formula
>
> flush(stderr()); flush(stdout())
>
> ### Name: survfit.formula
> ### Title: Compute a Survival Curve for Censored Data
> ### Aliases: survfit.formula [.survfit
> ### Keywords: survival
>
> ### ** Examples
>
> #fit a Kaplan-Meier and plot it
> fit <- survfit(Surv(time, status) ~ x, data = aml)
> plot(fit, lty = 2:3)
> legend(100, .8, c("Maintained", "Nonmaintained"), lty = 2:3)
>
> #fit a Cox proportional hazards model and plot the
> #predicted survival for a 60 year old
> fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian)
> plot(survfit(fit, newdata=data.frame(age=60)),
+ xscale=365.25, xlab = "Years", ylab="Survival")
>
> # Here is the data set from Turnbull
> # There are no interval censored subjects, only left-censored (status=3),
> # right-censored (status 0) and observed events (status 1)
> #
> # Time
> # 1 2 3 4
> # Type of observation
> # death 12 6 2 3
> # losses 3 2 0 3
> # late entry 2 4 2 5
> #
> tdata <- data.frame(time =c(1,1,1,2,2,2,3,3,3,4,4,4),
+ status=rep(c(1,0,2),4),
+ n =c(12,3,2,6,2,4,2,0,2,3,3,5))
> fit <- survfit(Surv(time, time, status, type='interval') ~1,
+ data=tdata, weight=n)
>
> #
> # Three curves for patients with monoclonal gammopathy.
> # 1. KM of time to PCM, ignoring death (statistically incorrect)
> # 2. Competing risk curves (also known as "cumulative incidence")
> # 3. Multi-state, showing Pr(in each state, at time t)
> #
> fitKM <- survfit(Surv(stop, event=='pcm') ~1, data=mgus1,
+ subset=(start==0))
> fitCR <- survfit(Surv(stop, event) ~1,
+ data=mgus1, subset=(start==0))
> fitMS <- survfit(Surv(start, stop, event) ~ 1, id=id, data=mgus1)
> ## Not run:
> ##D # CR curves show the competing risks
> ##D plot(fitCR, xscale=365.25, xmax=7300, mark.time=FALSE,
> ##D col=2:3, xlab="Years post diagnosis of MGUS",
> ##D ylab="P(state)")
> ##D lines(fitKM, fun='event', xmax=7300, mark.time=FALSE,
> ##D conf.int=FALSE)
> ##D text(3652, .4, "Competing risk: death", col=3)
> ##D text(5840, .15,"Competing risk: progression", col=2)
> ##D text(5480, .30,"KM:prog")
> ## End(Not run)
>
>
>
> cleanEx()
> nameEx("survfit.matrix")
> ### * survfit.matrix
>
> flush(stderr()); flush(stdout())
>
> ### Name: survfit.matrix
> ### Title: Create Aalen-Johansen estimates of multi-state survival from a
> ### matrix of hazards.
> ### Aliases: survfit.matrix
> ### Keywords: survival
>
> ### ** Examples
>
> etime <- with(mgus2, ifelse(pstat==0, futime, ptime))
> event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
> event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))
>
> cfit1 <- coxph(Surv(etime, event=="pcm") ~ age + sex, mgus2)
> cfit2 <- coxph(Surv(etime, event=="death") ~ age + sex, mgus2)
>
> # predicted competing risk curves for a 72 year old with mspike of 1.2
> # (median values), male and female.
> # The survfit call is a bit faster without standard errors.
> newdata <- expand.grid(sex=c("F", "M"), age=72, mspike=1.2)
>
> AJmat <- matrix(list(), 3,3)
> AJmat[1,2] <- list(survfit(cfit1, newdata, std.err=FALSE))
> AJmat[1,3] <- list(survfit(cfit2, newdata, std.err=FALSE))
> csurv <- survfit(AJmat, p0 =c(entry=1, PCM=0, death=0))
>
>
>
> cleanEx()
> nameEx("survobrien")
> ### * survobrien
>
> flush(stderr()); flush(stdout())
>
> ### Name: survobrien
> ### Title: O'Brien's Test for Association of a Single Variable with
> ### Survival
> ### Aliases: survobrien
> ### Keywords: survival
>
> ### ** Examples
>
> xx <- survobrien(Surv(futime, fustat) ~ age + factor(rx) + I(ecog.ps),
+ data=ovarian)
> coxph(Surv(time, status) ~ age + strata(.strata.), data=xx)
Call:
coxph(formula = Surv(time, status) ~ age + strata(.strata.),
data = xx)
coef exp(coef) se(coef) z p
age 0.5681 1.7649 0.1816 3.128 0.00176
Likelihood ratio test=10.55 on 1 df, p=0.001165
n= 230, number of events= 12
>
>
>
> cleanEx()
> nameEx("survreg")
> ### * survreg
>
> flush(stderr()); flush(stdout())
>
> ### Name: survreg
> ### Title: Regression for a Parametric Survival Model
> ### Aliases: survreg model.frame.survreg labels.survreg print.survreg.penal
> ### print.summary.survreg survReg anova.survreg anova.survreglist
> ### Keywords: survival
>
> ### ** Examples
>
> # Fit an exponential model: the two fits are the same
> survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='weibull',
+ scale=1)
Call:
survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian,
dist = "weibull", scale = 1)
Coefficients:
(Intercept) ecog.ps rx
6.9618376 -0.4331347 0.5815027
Scale fixed at 1
Loglik(model)= -97.2 Loglik(intercept only)= -98
Chisq= 1.67 on 2 degrees of freedom, p= 0.434
n= 26
> survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian,
+ dist="exponential")
Call:
survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian,
dist = "exponential")
Coefficients:
(Intercept) ecog.ps rx
6.9618376 -0.4331347 0.5815027
Scale fixed at 1
Loglik(model)= -97.2 Loglik(intercept only)= -98
Chisq= 1.67 on 2 degrees of freedom, p= 0.434
n= 26
>
> #
> # A model with different baseline survival shapes for two groups, i.e.,
> # two different scale parameters
> survreg(Surv(time, status) ~ ph.ecog + age + strata(sex), lung)
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + age + strata(sex),
data = lung)
Coefficients:
(Intercept) ph.ecog age
6.73234505 -0.32443043 -0.00580889
Scale:
sex=1 sex=2
0.7834211 0.6547830
Loglik(model)= -1137.3 Loglik(intercept only)= -1146.2
Chisq= 17.8 on 2 degrees of freedom, p= 0.000137
n=227 (1 observation deleted due to missingness)
>
> # There are multiple ways to parameterize a Weibull distribution. The survreg
> # function embeds it in a general location-scale family, which is a
> # different parameterization than the rweibull function, and often leads
> # to confusion.
> # survreg's scale = 1/(rweibull shape)
> # survreg's intercept = log(rweibull scale)
> # For the log-likelihood all parameterizations lead to the same value.
> y <- rweibull(1000, shape=2, scale=5)
> survreg(Surv(y)~1, dist="weibull")
Call:
survreg(formula = Surv(y) ~ 1, dist = "weibull")
Coefficients:
(Intercept)
1.604435
Scale= 0.4965001
Loglik(model)= -2199.4 Loglik(intercept only)= -2199.4
n= 1000
>
> # Economists fit a model called `tobit regression', which is a standard
> # linear regression with Gaussian errors, and left censored data.
> tobinfit <- survreg(Surv(durable, durable>0, type='left') ~ age + quant,
+ data=tobin, dist='gaussian')
>
>
>
> cleanEx()
> nameEx("survreg.distributions")
> ### * survreg.distributions
>
> flush(stderr()); flush(stdout())
>
> ### Name: survreg.distributions
> ### Title: Parametric Survival Distributions
> ### Aliases: survreg.distributions
> ### Keywords: survival
>
> ### ** Examples
>
> # time transformation
> survreg(Surv(time, status) ~ ph.ecog + sex, dist='weibull', data=lung)
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + sex, data = lung,
dist = "weibull")
Coefficients:
(Intercept) ph.ecog sex
5.8195907 -0.3557319 0.4013684
Scale= 0.7310495
Loglik(model)= -1133.1 Loglik(intercept only)= -1147.4
Chisq= 28.73 on 2 degrees of freedom, p= 5.76e-07
n=227 (1 observation deleted due to missingness)
> # change the transformation to work in years
> # intercept changes by log(365), everything else stays the same
> my.weibull <- survreg.distributions$weibull
> my.weibull$trans <- function(y) log(y/365)
> my.weibull$itrans <- function(y) 365*exp(y)
> survreg(Surv(time, status) ~ ph.ecog + sex, lung, dist=my.weibull)
Call:
survreg(formula = Surv(time, status) ~ ph.ecog + sex, data = lung,
dist = my.weibull)
Coefficients:
(Intercept) ph.ecog sex
-0.08030664 -0.35573188 0.40136844
Scale= 0.7310495
Loglik(model)= -1133.1 Loglik(intercept only)= -1147.4
Chisq= 28.73 on 2 degrees of freedom, p= 5.76e-07
n=227 (1 observation deleted due to missingness)
>
> # Weibull parametrisation
> y<-rweibull(1000, shape=2, scale=5)
> survreg(Surv(y)~1, dist="weibull")
Call:
survreg(formula = Surv(y) ~ 1, dist = "weibull")
Coefficients:
(Intercept)
1.604435
Scale= 0.4965001
Loglik(model)= -2199.4 Loglik(intercept only)= -2199.4
n= 1000
> # survreg scale parameter maps to 1/shape, linear predictor to log(scale)
>
> # Cauchy fit
> mycauchy <- list(name='Cauchy',
+ init= function(x, weights, ...)
+ c(median(x), mad(x)),
+ density= function(x, parms) {
+ temp <- 1/(1 + x^2)
+ cbind(.5 + atan(x)/pi, .5+ atan(-x)/pi,
+ temp/pi, -2 *x*temp, 2*temp*(4*x^2*temp -1))
+ },
+ quantile= function(p, parms) tan((p-.5)*pi),
+ deviance= function(...) stop('deviance residuals not defined')
+ )
> survreg(Surv(log(time), status) ~ ph.ecog + sex, lung, dist=mycauchy)
Call:
survreg(formula = Surv(log(time), status) ~ ph.ecog + sex, data = lung,
dist = mycauchy)
Coefficients:
(Intercept) ph.ecog sex
5.4517240 -0.3979387 0.4692383
Scale= 0.4788955
Loglik(model)= -274.6 Loglik(intercept only)= -294.9
Chisq= 40.75 on 2 degrees of freedom, p= 1.42e-09
n=227 (1 observation deleted due to missingness)
>
>
>
> cleanEx()
> nameEx("survregDtest")
> ### * survregDtest
>
> flush(stderr()); flush(stdout())
>
> ### Name: survregDtest
> ### Title: Verify a survreg distribution
> ### Aliases: survregDtest
> ### Keywords: survival
>
> ### ** Examples
>
> # An invalid distribution (it should have "init =" on line 2)
> # surveg would give an error message
> mycauchy <- list(name='Cauchy',
+ init<- function(x, weights, ...)
+ c(median(x), mad(x)),
+ density= function(x, parms) {
+ temp <- 1/(1 + x^2)
+ cbind(.5 + atan(temp)/pi, .5+ atan(-temp)/pi,
+ temp/pi, -2 *x*temp, 2*temp^2*(4*x^2*temp -1))
+ },
+ quantile= function(p, parms) tan((p-.5)*pi),
+ deviance= function(...) stop('deviance residuals not defined')
+ )
>
> survregDtest(mycauchy, TRUE)
[1] "Missing or invalid init function"
>
>
>
> cleanEx()
> nameEx("tcut")
> ### * tcut
>
> flush(stderr()); flush(stdout())
>
> ### Name: tcut
> ### Title: Factors for person-year calculations
> ### Aliases: tcut [.tcut levels.tcut
> ### Keywords: survival
>
> ### ** Examples
>
> # For pyears, all time variable need to be on the same scale; but
> # futime is in months and age is in years
> test <- mgus2
> test$years <- test$futime/30.5 # follow-up in years
>
> # first grouping based on years from starting age (= current age)
> # second based on years since enrollment (all start at 0)
> test$agegrp <- tcut(test$age, c(0,60, 70, 80, 100),
+ c("<=60", "60-70", "70-80", ">80"))
> test$fgrp <- tcut(rep(0, nrow(test)), c(0, 1, 5, 10, 100),
+ c("0-1yr", "1-5yr", "5-10yr", ">10yr"))
>
> # death rates per 1000, by age group
> pfit1 <- pyears(Surv(years, death) ~ agegrp, scale =1000, data=test)
> round(pfit1$event/ pfit1$pyears)
agegrp
<=60 60-70 70-80 >80
89 128 245 479
>
> #death rates per 100, by follow-up year and age
> # there are excess deaths in the first year, within each age stratum
> pfit2 <- pyears(Surv(years, death) ~ fgrp + agegrp, scale =1000, data=test)
> round(pfit2$event/ pfit2$pyears)
agegrp
fgrp <=60 60-70 70-80 >80
0-1yr 139 137 234 437
1-5yr 68 117 241 499
5-10yr 91 145 300 476
>10yr 0 424 0 2440
>
>
>
> cleanEx()
> nameEx("tmerge")
> ### * tmerge
>
> flush(stderr()); flush(stdout())
>
> ### Name: tmerge
> ### Title: Time based merge for survival data
> ### Aliases: tmerge
> ### Keywords: survival
>
> ### ** Examples
>
> # The pbc data set contains baseline data and follow-up status
> # for a set of subjects with primary biliary cirrhosis, while the
> # pbcseq data set contains repeated laboratory values for those
> # subjects.
> # The first data set contains data on 312 subjects in a clinical trial plus
> # 106 that agreed to be followed off protocol, the second data set has data
> # only on the trial subjects.
> temp <- subset(pbc, id <= 312, select=c(id:sex, stage)) # baseline data
> pbc2 <- tmerge(temp, temp, id=id, endpt = event(time, status))
> pbc2 <- tmerge(pbc2, pbcseq, id=id, ascites = tdc(day, ascites),
+ bili = tdc(day, bili), albumin = tdc(day, albumin),
+ protime = tdc(day, protime), alk.phos = tdc(day, alk.phos))
>
> fit <- coxph(Surv(tstart, tstop, endpt==2) ~ protime + log(bili), data=pbc2)
>
>
>
> cleanEx()
> nameEx("tobin")
> ### * tobin
>
> flush(stderr()); flush(stdout())
>
> ### Name: tobin
> ### Title: Tobin's Tobit data
> ### Aliases: tobin
> ### Keywords: datasets
>
> ### ** Examples
>
> tfit <- survreg(Surv(durable, durable>0, type='left') ~age + quant,
+ data=tobin, dist='gaussian')
>
> predict(tfit,type="response")
[1] -3.04968679 -4.31254182 -0.54163315 -0.25607164 -1.85017727 -2.40987803
[7] -3.50629220 -0.74041486 -4.05145594 -3.55880518 -0.32223237 -3.68044619
[13] -3.65997456 -2.63254564 0.22382063 0.02177674 -0.09571284 -3.17696755
[19] -0.61521215 -3.13913903
>
>
>
>
> cleanEx()
> nameEx("transplant")
> ### * transplant
>
> flush(stderr()); flush(stdout())
>
> ### Name: transplant
> ### Title: Liver transplant waiting list
> ### Aliases: transplant
> ### Keywords: datasets
>
> ### ** Examples
>
> #since event is a factor, survfit creates competing risk curves
> pfit <- survfit(Surv(futime, event) ~ abo, transplant)
> pfit[,2] #time to liver transplant, by blood type
Call: survfit(formula = Surv(futime, event) ~ abo, data = transplant)
n nevent rmean*
abo=A, death 325 21 164.9734
abo=B, death 103 10 202.4902
abo=AB, death 41 3 137.8293
abo=O, death 346 32 182.1075
*restricted mean time in state (max time = 2055 )
> plot(pfit[,2], mark.time=FALSE, col=1:4, lwd=2, xmax=735,
+ xscale=30.5, xlab="Months", ylab="Fraction transplanted",
+ xaxt = 'n')
> temp <- c(0, 6, 12, 18, 24)
> axis(1, temp*30.5, temp)
> legend(450, .35, levels(transplant$abo), lty=1, col=1:4, lwd=2)
>
> # competing risks for type O
> plot(pfit[4,], xscale=30.5, xmax=735, col=1:3, lwd=2)
> legend(450, .4, c("Death", "Transpant", "Withdrawal"), col=1:3, lwd=2)
>
>
>
> cleanEx()
> nameEx("udca")
> ### * udca
>
> flush(stderr()); flush(stdout())
>
> ### Name: udca
> ### Title: Data from a trial of usrodeoxycholic acid
> ### Aliases: udca udca1 udca2
> ### Keywords: datasets
>
> ### ** Examples
>
> # values found in table 8.3 of the book
> fit1 <- coxph(Surv(futime, status) ~ trt + log(bili) + stage,
+ cluster =id , data=udca1)
> fit2 <- coxph(Surv(futime, status) ~ trt + log(bili) + stage +
+ strata(endpoint), cluster=id, data=udca2)
>
>
>
>
> cleanEx()
> nameEx("untangle.specials")
> ### * untangle.specials
>
> flush(stderr()); flush(stdout())
>
> ### Name: untangle.specials
> ### Title: Help Process the 'specials' Argument of the 'terms' Function.
> ### Aliases: untangle.specials
> ### Keywords: survival
>
> ### ** Examples
>
> formula <- Surv(tt,ss) ~ x + z*strata(id)
> tms <- terms(formula, specials="strata")
> ## the specials attribute
> attr(tms, "specials")
$strata
[1] 4
> ## main effects
> untangle.specials(tms, "strata")
$vars
[1] "strata(id)"
$tvar
[1] 3
$terms
[1] 3
> ## and interactions
> untangle.specials(tms, "strata", order=1:2)
$vars
[1] "strata(id)"
$tvar
[1] 3
$terms
[1] 3 4
>
>
>
> cleanEx()
> nameEx("uspop2")
> ### * uspop2
>
> flush(stderr()); flush(stdout())
>
> ### Name: uspop2
> ### Title: Projected US Population
> ### Aliases: uspop2
> ### Keywords: datasets
>
> ### ** Examples
>
> us50 <- uspop2[51:101,, "2000"] #US 2000 population, 50 and over
> age <- as.integer(dimnames(us50)[[1]])
> smat <- model.matrix( ~ factor(floor(age/5)) -1)
> ustot <- t(smat) %*% us50 #totals by 5 year age groups
> temp <- c(50,55, 60, 65, 70, 75, 80, 85, 90, 95)
> dimnames(ustot) <- list(c(paste(temp, temp+4, sep="-"), "100+"),
+ c("male", "female"))
>
>
>
> cleanEx()
> nameEx("xtfrm.Surv")
> ### * xtfrm.Surv
>
> flush(stderr()); flush(stdout())
>
> ### Name: xtfrm.Surv
> ### Title: Sorting order for Surv objects
> ### Aliases: xtfrm.Surv sort.Surv order.Surv
> ### Keywords: survival
>
> ### ** Examples
>
> test <- c(Surv(c(10, 9,9, 8,8,8,7,5,5,4), rep(1:0, 5)), Surv(6.2, NA))
> test
[1] 10.0 9.0+ 9.0 8.0+ 8.0 8.0+ 7.0 5.0+ 5.0 4.0+ 6.2?
> sort(test)
[1] 4+ 5 5+ 7 8 8+ 8+ 9 9+ 10
>
>
>
> cleanEx()
> nameEx("yates")
> ### * yates
>
> flush(stderr()); flush(stdout())
>
> ### Name: yates
> ### Title: Population prediction
> ### Aliases: yates
> ### Keywords: models survival
>
> ### ** Examples
>
> fit1 <- lm(skips ~ Solder*Opening + Mask, data = solder)
> yates(fit1, ~Opening, population = "factorial")
Opening pmm std test chisq df ss Pr
L 3.2638 0.33460 global 573.4 2 15980 < 1e-08
M 3.5700 0.30480
S 12.3519 0.31251
>
> fit2 <- coxph(Surv(time, status) ~ factor(ph.ecog)*sex + age, lung)
> yates(fit2, ~ ph.ecog, predict="risk") # hazard ratio
factor(ph.ecog) pmm std test chisq df Pr
0 0.94238 0.46334 factor(ph.ecog) NA NA NA
1 1.42677 0.75697
2 1.74848 3.80017
3 NA 74.01221
>
>
>
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 4.665 0.351 5.018 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')