334 lines
10 KiB
Plaintext
334 lines
10 KiB
Plaintext
|
|
||
|
|
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
|