94 lines
3.4 KiB
Plaintext
94 lines
3.4 KiB
Plaintext
|
|
||
|
|
R Under development (unstable) (2020-03-27 r78086) -- "Unsuffered Consequences"
|
||
|
|
Copyright (C) 2020 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.
|
||
|
|
|
||
|
|
> library(survival)
|
||
|
|
> options(na.action=na.exclude, contrasts=c('contr.treatment', 'contr.poly'))
|
||
|
|
>
|
||
|
|
> # Verify stratified fits in a simple way, but combining two data
|
||
|
|
> # sets and doing a single fit
|
||
|
|
> #
|
||
|
|
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
|
||
|
|
>
|
||
|
|
> tdata <- data.frame(time=c(lung$time, ovarian$futime),
|
||
|
|
+ status=c(lung$status-1, ovarian$fustat),
|
||
|
|
+ group =rep(0:1, c(nrow(lung), nrow(ovarian))))
|
||
|
|
> fit1 <- survreg(Surv(time, status) ~ 1, lung)
|
||
|
|
> fit2 <- survreg(Surv(futime, fustat) ~ 1, ovarian)
|
||
|
|
> fit3 <- survreg(Surv(time, status) ~ group + strata(group), tdata)
|
||
|
|
>
|
||
|
|
> aeq(c(fit1$coef, fit2$coef-fit1$coef), fit3$coef)
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(c(fit1$scale, fit2$scale), fit3$scale)
|
||
|
|
[1] TRUE
|
||
|
|
> aeq(fit1$loglik[2] + fit2$loglik[2], fit3$loglik[2])
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> #
|
||
|
|
> # Test out the cluster term in survreg, which means first a test
|
||
|
|
> # of the dfbeta residuals
|
||
|
|
> # I also am checking that missing values propogate
|
||
|
|
> 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),
|
||
|
|
+ id= 1:7)
|
||
|
|
> fit1 <- survreg(Surv(time, status) ~ x, cluster = id, test1)
|
||
|
|
> fit2 <- survreg(Surv(time, status) ~ x + cluster(id), test1) #old form
|
||
|
|
> all.equal(fit1, fit2)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> db1 <- resid(fit1, 'dfbeta')
|
||
|
|
> ijack <-db1
|
||
|
|
> eps <- 1e-7
|
||
|
|
> for (i in 1:7) {
|
||
|
|
+ temp <- rep(1.0,7)
|
||
|
|
+ temp[i] <- 1-eps
|
||
|
|
+ tfit <- survreg(Surv(time, status) ~ x, test1, weight=temp)
|
||
|
|
+ ijack[i,] <- c(tfit$coef, log(tfit$scale))
|
||
|
|
+ }
|
||
|
|
> ijack[2,] <- NA # stick the NA back in
|
||
|
|
> ijack <- (rep(c(fit1$coef, log(fit1$scale)), each=nrow(db1)) - ijack)/eps
|
||
|
|
> all.equal(db1, ijack, tolerance=eps)
|
||
|
|
[1] TRUE
|
||
|
|
> all.equal(t(db1[-2,])%*% db1[-2,], fit1$var)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
> # This is a harder test since there are multiple strata and multiple
|
||
|
|
> # obs/subject. Use of enum + strata(enum) in essenence fits a different
|
||
|
|
> # baseline Weibull to each strata, with common coefficients for rx, size, and
|
||
|
|
> # number.
|
||
|
|
> fit1 <- survreg(Surv(stop-start, event) ~ rx + size + number +
|
||
|
|
+ factor(enum) + strata(enum), data=bladder2)
|
||
|
|
>
|
||
|
|
> db1 <- resid(fit1, type='dfbeta', collapse=bladder2$id)
|
||
|
|
> ijack <- db1 # a matrix of the same size
|
||
|
|
> for (i in 1:nrow(db1)) {
|
||
|
|
+ twt <- rep(1., nrow(bladder2))
|
||
|
|
+ twt[bladder2$id==i] <- 1-eps
|
||
|
|
+ tfit <- survreg(Surv(stop-start, event) ~ rx + size + number +
|
||
|
|
+ factor(enum) + strata(enum), data=bladder2,
|
||
|
|
+ weight=twt)
|
||
|
|
+ ijack[i,] <- c(coef(tfit), log(tfit$scale))
|
||
|
|
+ }
|
||
|
|
> ijack <- (rep(c(fit1$coef, log(fit1$scale)), each=nrow(db1)) - ijack)/eps
|
||
|
|
> all.equal(db1, ijack, tolerance=eps*2)
|
||
|
|
[1] TRUE
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> proc.time()
|
||
|
|
user system elapsed
|
||
|
|
1.207 0.065 1.262
|