66 lines
2.5 KiB
R
66 lines
2.5 KiB
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)
|
|
aeq(c(fit1$scale, fit2$scale), fit3$scale)
|
|
aeq(fit1$loglik[2] + fit2$loglik[2], fit3$loglik[2])
|
|
|
|
#
|
|
# 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)
|
|
|
|
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)
|
|
all.equal(t(db1[-2,])%*% db1[-2,], fit1$var)
|
|
|
|
# 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)
|
|
|