R Under development (unstable) (2023-08-12 r84939 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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.90 0.12 1.01