File: C:/Program Files/R/R-4.4.0/library/survival/tests/testreg.Rout.save
R Under development (unstable) (2019-06-28 r76752) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> options(na.action=na.exclude) #preserve length of missings
> library(survival)
>
> #
> # Run a test that can be verified using other packages
> #
> test1 <- data.frame(time= c(9, 3,1,1,6,6,8),
+ status=c(1,NA,1,0,1,1,0),
+ x= c(0, 2,1,1,1,0,0))
> fit1w <- survreg(Surv(time, status) ~x, test1, dist='weibull')
> fit1w
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "weibull")
Coefficients:
(Intercept) x
2.2373335 -0.7442249
Scale= 0.4563163
Loglik(model)= -10.3 Loglik(intercept only)= -11.4
Chisq= 2.22 on 1 degrees of freedom, p= 0.136
n=6 (1 observation deleted due to missingness)
> summary(fit1w)
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "weibull")
Value Std. Error z p
(Intercept) 2.237 0.330 6.78 1.2e-11
x -0.744 0.486 -1.53 0.13
Log(scale) -0.785 0.433 -1.81 0.07
Scale= 0.456
Weibull distribution
Loglik(model)= -10.3 Loglik(intercept only)= -11.4
Chisq= 2.22 on 1 degrees of freedom, p= 0.14
Number of Newton-Raphson Iterations: 8
n=6 (1 observation deleted due to missingness)
>
> fit1e <- survreg(Surv(time, status) ~x, test1, dist='exponential')
> fit1e
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "exponential")
Coefficients:
(Intercept) x
2.442347 -1.056053
Scale fixed at 1
Loglik(model)= -11.7 Loglik(intercept only)= -12.2
Chisq= 1.07 on 1 degrees of freedom, p= 0.302
n=6 (1 observation deleted due to missingness)
> summary(fit1e)
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "exponential")
Value Std. Error z p
(Intercept) 2.442 0.707 3.45 0.00055
x -1.056 1.000 -1.06 0.29094
Scale fixed at 1
Exponential distribution
Loglik(model)= -11.7 Loglik(intercept only)= -12.2
Chisq= 1.07 on 1 degrees of freedom, p= 0.3
Number of Newton-Raphson Iterations: 4
n=6 (1 observation deleted due to missingness)
>
> fit1l <- survreg(Surv(time, status) ~x, test1, dist='loglogistic')
> fit1l
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "loglogistic")
Coefficients:
(Intercept) x
2.177208 -1.195672
Scale= 0.3847582
Loglik(model)= -10.7 Loglik(intercept only)= -12
Chisq= 2.7 on 1 degrees of freedom, p= 0.1
n=6 (1 observation deleted due to missingness)
> summary(fit1l)
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "loglogistic")
Value Std. Error z p
(Intercept) 2.177 0.365 5.96 2.5e-09
x -1.196 0.711 -1.68 0.093
Log(scale) -0.955 0.396 -2.41 0.016
Scale= 0.385
Log logistic distribution
Loglik(model)= -10.7 Loglik(intercept only)= -12
Chisq= 2.7 on 1 degrees of freedom, p= 0.1
Number of Newton-Raphson Iterations: 4
n=6 (1 observation deleted due to missingness)
>
> fit1g <- survreg(Surv(time, status) ~x, test1, dist='lognormal')
> summary(fit1g)
Call:
survreg(formula = Surv(time, status) ~ x, data = test1, dist = "lognormal")
Value Std. Error z p
(Intercept) 2.210 0.404 5.48 4.4e-08
x -1.268 0.585 -2.17 0.03
Log(scale) -0.446 0.342 -1.30 0.19
Scale= 0.64
Log Normal distribution
Loglik(model)= -10.5 Loglik(intercept only)= -12.1
Chisq= 3.26 on 1 degrees of freedom, p= 0.071
Number of Newton-Raphson Iterations: 5
n=6 (1 observation deleted due to missingness)
> #
> # Do a test with the ovarian data
> #
> fitfw <- survreg(Surv(futime, fustat) ~ age + ecog.ps, ovarian,
+ dist='weibull')
> fitfw
Call:
survreg(formula = Surv(futime, fustat) ~ age + ecog.ps, data = ovarian,
dist = "weibull")
Coefficients:
(Intercept) age ecog.ps
12.28496723 -0.09702669 0.09977342
Scale= 0.6032744
Loglik(model)= -90 Loglik(intercept only)= -98
Chisq= 15.98 on 2 degrees of freedom, p= 0.000339
n= 26
>
> fitfl <- survreg(Surv(futime, fustat) ~ age + ecog.ps, ovarian,
+ dist='loglogistic')
> fitfl
Call:
survreg(formula = Surv(futime, fustat) ~ age + ecog.ps, data = ovarian,
dist = "loglogistic")
Coefficients:
(Intercept) age ecog.ps
11.50853384 -0.08876814 0.09033348
Scale= 0.4464064
Loglik(model)= -89.5 Loglik(intercept only)= -97.4
Chisq= 15.67 on 2 degrees of freedom, p= 0.000395
n= 26
>
> #test out interval censoring, using some dummy time values
>
> idat <- read.table('data.interval', skip=3, header=T, sep=',')
> flsurv<- Surv(idat$ltime, idat$rtime, type='interval2')
>
> fitfw2 <- survreg(flsurv ~ age + ecog.ps, idat, dist='weibull')
> summary(fitfw2)
Call:
survreg(formula = flsurv ~ age + ecog.ps, data = idat, dist = "weibull")
Value Std. Error z p
(Intercept) 12.3886 1.6027 7.73 1.1e-14
age -0.0986 0.0254 -3.89 0.0001
ecog.ps 0.0971 0.3776 0.26 0.7971
Log(scale) -0.4773 0.2583 -1.85 0.0647
Scale= 0.62
Weibull distribution
Loglik(model)= -56.2 Loglik(intercept only)= -64
Chisq= 15.57 on 2 degrees of freedom, p= 0.00042
Number of Newton-Raphson Iterations: 6
n= 26
>
> fitfl2 <- survreg(flsurv ~ age + ecog.ps, idat, dist='loglogistic')
> summary(fitfl2)
Call:
survreg(formula = flsurv ~ age + ecog.ps, data = idat, dist = "loglogistic")
Value Std. Error z p
(Intercept) 11.5268 1.5283 7.54 4.6e-14
age -0.0888 0.0240 -3.70 0.00021
ecog.ps 0.0818 0.3642 0.22 0.82219
Log(scale) -0.8023 0.2706 -2.96 0.00303
Scale= 0.448
Log logistic distribution
Loglik(model)= -55.9 Loglik(intercept only)= -63.5
Chisq= 15.35 on 2 degrees of freedom, p= 0.00046
Number of Newton-Raphson Iterations: 5
n= 26
>
> fitfg2 <- survreg(flsurv ~ age + ecog.ps, idat, dist='lognormal')
> summary(fitfg2)
Call:
survreg(formula = flsurv ~ age + ecog.ps, data = idat, dist = "lognormal")
Value Std. Error z p
(Intercept) 11.1548 1.4347 7.77 7.6e-15
age -0.0855 0.0238 -3.60 0.00032
ecog.ps 0.2066 0.3828 0.54 0.58945
Log(scale) -0.2297 0.2508 -0.92 0.35972
Scale= 0.795
Log Normal distribution
Loglik(model)= -56 Loglik(intercept only)= -63.5
Chisq= 14.94 on 2 degrees of freedom, p= 0.00057
Number of Newton-Raphson Iterations: 5
n= 26
>
> logt <- c(survreg.distributions$t,
+ survreg.distributions$weibull[c('trans', 'itrans', 'dtrans')])
> logt$name <- 'log(t)'
>
> fitft2 <- survreg(Surv(ltime, rtime, type='interval2') ~ age + ecog.ps,
+ idat, dist=logt, parm=100)
> summary(fitft2) #should be quite close to fitfg2
Call:
survreg(formula = Surv(ltime, rtime, type = "interval2") ~ age +
ecog.ps, data = idat, dist = logt, parms = 100)
Value Std. Error z p
(Intercept) 11.1856 1.4419 7.76 8.7e-15
age -0.0858 0.0238 -3.61 0.00031
ecog.ps 0.1978 0.3814 0.52 0.60401
Log(scale) -0.2394 0.2522 -0.95 0.34254
Scale= 0.787
log(t) distribution: parmameters= 100
Loglik(model)= -56 Loglik(intercept only)= -63.5
Chisq= 14.97 on 2 degrees of freedom, p= 0.00056
Number of Newton-Raphson Iterations: 5
n= 26
>
> #
> # Check out the survreg density and probability functions
> #
>
> # Gaussian
> x <- -10:10
> p <- seq(.1, .95, length=25)
> all.equal(dsurvreg(x, 1, 5, 'gaussian'), dnorm(x, 1, 5))
[1] TRUE
> all.equal(psurvreg(x, 1, 5, 'gaussian'), pnorm(x, 1, 5))
[1] TRUE
> all.equal(qsurvreg(p, 1, 5, 'gaussian'), qnorm(p, 1, 5))
[1] TRUE
>
> # Lognormal
> x <- 1:10
> all.equal(dsurvreg(x, 1, 5, 'lognormal'), dlnorm(x, 1, 5))
[1] TRUE
> all.equal(psurvreg(x, 1, 5, 'lognormal'), plnorm(x, 1, 5))
[1] TRUE
> all.equal(qsurvreg(p, 1, 5, 'lognormal'), qlnorm(p, 1, 5))
[1] TRUE
>
> # Weibull
> lambda <- exp(-2)
> rho <- 1/3
> temp <- (lambda*x)^rho
> all.equal(psurvreg(x, 2, 3), 1- exp(-temp))
[1] TRUE
> all.equal(dsurvreg(x, 2, 3), lambda*rho*(lambda*x)^(rho-1)*exp(-temp))
[1] TRUE
>
>
> # verify labeling in the vcov function, with 0, 1, or 2 scale factors
> fit0 <- survreg(Surv(time, status) ~ age + ph.ecog, lung, dist='exponential')
> vcov(fit0)
(Intercept) age ph.ecog
(Intercept) 0.330860114 -0.0051541551 0.001441498
age -0.005154155 0.0000854214 -0.000240252
ph.ecog 0.001441498 -0.0002402520 0.013051391
> fit1 <- survreg(Surv(time, status) ~ age + ph.ecog, lung, dist='weibull')
> vcov(fit1)
(Intercept) age ph.ecog Log(scale)
(Intercept) 0.1837842638 -2.858745e-03 0.0003634856 2.800938e-03
age -0.0028587447 4.731485e-05 -0.0001279049 -3.517281e-05
ph.ecog 0.0003634856 -1.279049e-04 0.0073616277 -7.091072e-04
Log(scale) 0.0028009379 -3.517281e-05 -0.0007091072 3.837229e-03
>
> fit2 <- survreg(Surv(time, status) ~ age + ph.ecog + strata(sex), lung)
> vcov(fit2)
(Intercept) age ph.ecog Log(scale[sex=1])
(Intercept) 0.179740626 -2.845266e-03 0.0016839233 -1.147693e-03
age -0.002845266 4.803051e-05 -0.0001515280 4.208343e-05
ph.ecog 0.001683923 -1.515280e-04 0.0074813166 -6.291573e-04
Log(scale[sex=1]) -0.001147693 4.208343e-05 -0.0006291573 6.272257e-03
Log(scale[sex=2]) 0.009087054 -1.600040e-04 -0.0008306451 -4.558932e-04
Log(scale[sex=2])
(Intercept) 0.0090870545
age -0.0001600040
ph.ecog -0.0008306451
Log(scale[sex=1]) -0.0004558932
Log(scale[sex=2]) 0.0113833096
>
> proc.time()
user system elapsed
0.804 0.044 0.853