# tests/r_resid.R In survival: Survival Analysis

```options(na.action=na.exclude) # preserve missings
options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
library(survival)

aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)

fit1 <- survreg(Surv(futime, fustat) ~ age + ecog.ps, ovarian)
fit4 <- survreg(Surv(log(futime), fustat) ~age + ecog.ps, ovarian,
dist='extreme')

print(fit1)
summary(fit4)

# Hypothesis (and I'm fairly sure): censorReg shares the fault of many
#  iterative codes -- it returns the loglik and variance for iteration k
#  but the coef vector of iteration k+1.  Hence the "all.equal" tests
#  below don't come out perfect.
#
if (exists('censorReg')) {  #true for Splus, not R
fit2 <- censorReg(censor(futime, fustat) ~ age + ecog.ps, ovarian)
fit3 <- survreg(Surv(futime, fustat) ~ age + ecog.ps, ovarian,
iter=0, init=c(fit2\$coef,   log(fit2\$scale)))

aeq(resid(fit2, type='working')[,1], resid(fit3, type='working'))
aeq(resid(fit2, type='response')[,1], resid(fit3, type='response'))

temp <- sign(resid(fit3, type='working'))
aeq(resid(fit2, type='deviance')[,1],
temp*abs(resid(fit3, type='deviance')))
aeq(resid(fit2, type='deviance')[,1], resid(fit3, type='deviance'))
}
#
# Now check fit1 and fit4, which should follow identical iteration paths
#   These tests should all be true
#
aeq(fit1\$coef, fit4\$coef)

resid(fit1, type='working')
resid(fit1, type='response')
resid(fit1, type='deviance')
resid(fit1, type='dfbeta')
resid(fit1, type='dfbetas')
resid(fit1, type='ldcase')
resid(fit1, type='ldresp')
resid(fit1, type='ldshape')
resid(fit1, type='matrix')

aeq(resid(fit1, type='working'),resid(fit4, type='working'))
#aeq(resid(fit1, type='response'), resid(fit4, type='response'))#should differ
aeq(resid(fit1, type='deviance'), resid(fit4, type='deviance'))
aeq(resid(fit1, type='dfbeta'),   resid(fit4, type='dfbeta'))
aeq(resid(fit1, type='dfbetas'),  resid(fit4, type='dfbetas'))
aeq(resid(fit1, type='ldcase'),   resid(fit4, type='ldcase'))
aeq(resid(fit1, type='ldresp'),   resid(fit4, type='ldresp'))
aeq(resid(fit1, type='ldshape'),  resid(fit4, type='ldshape'))
aeq(resid(fit1, type='matrix'),   resid(fit4, type='matrix'))
#
# Some tests of the quantile residuals
#
# These should agree exactly with Ripley and Venables' book
fit1 <- survreg(Surv(time, status) ~ temp, data=imotor)
summary(fit1)

#
# The first prediction has the SE that I think is correct
#  The third is the se found in an early draft of Ripley; fit1 ignoring
#  the variation in scale estimate, except via it's impact on the
#  upper left corner of the inverse information matrix.
# Numbers 1 and 3 differ little for this dataset
#
predict(fit1, data.frame(temp=130), type='uquantile', p=c(.5, .1), se=T)

fit2 <- survreg(Surv(time, status) ~ temp, data=imotor, scale=fit1\$scale)
predict(fit2, data.frame(temp=130), type='uquantile', p=c(.5, .1), se=T)

fit3 <- fit2
fit3\$var <- fit1\$var[1:2,1:2]
predict(fit3, data.frame(temp=130), type='uquantile', p=c(.5, .1), se=T)

pp <- seq(.05, .7, length=40)
xx <- predict(fit1, data.frame(temp=130), type='uquantile', se=T,
p=pp)
#matplot(pp, cbind(xx\$fit, xx\$fit+2*xx\$se, xx\$fit - 2*xx\$se), type='l')

#
# Now try out the various combinations of strata, #predicted, and
#  number of quantiles desired
#
fit1 <- survreg(Surv(time, status) ~ inst + strata(inst) + age + sex, lung)
qq1 <- predict(fit1, type='quantile', p=.3, se=T)
qq2 <- predict(fit1, type='quantile', p=c(.2, .3, .4), se=T)
aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
aeq(qq1\$fit, qq2\$fit[,2])
aeq(qq1\$se.fit, qq2\$se.fit[,2])

qq3 <- predict(fit1, type='quantile', p=c(.2, .3, .4), se=T,
newdata= lung[1:5,])
aeq(qq3\$fit, qq2\$fit[1:5,])

qq4 <- predict(fit1, type='quantile', p=c(.2, .3, .4), se=T, newdata=lung[7,])
aeq(qq4\$fit, qq2\$fit[7,])

qq5 <- predict(fit1, type='quantile', p=c(.2, .3, .4), se=T, newdata=lung)
aeq(qq2\$fit, qq5\$fit)
aeq(qq2\$se.fit, qq5\$se.fit)
```

## Try the survival package in your browser

Any scripts or data that you put into this service are public.

survival documentation built on Feb. 16, 2023, 7:34 p.m.