61 lines
2.3 KiB
R
61 lines
2.3 KiB
R
|
options(na.action=na.exclude) # preserve missings
|
||
|
|
options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
|
||
|
|
library(survival)
|
||
|
|
|
||
|
|
#
|
||
|
|
# Test out the strata capabilities
|
||
|
|
#
|
||
|
|
tol <- survreg.control()$rel.tolerance
|
||
|
|
aeq <- function(x,y,...) all.equal(as.vector(x), as.vector(y), ...)
|
||
|
|
|
||
|
|
# intercept only models
|
||
|
|
fit1 <- survreg(Surv(time, status) ~ strata(sex), lung)
|
||
|
|
fit2 <- survreg(Surv(time, status) ~ strata(sex) + sex, lung)
|
||
|
|
fit3a<- survreg(Surv(time,status) ~1, lung, subset=(sex==1))
|
||
|
|
fit3b<- survreg(Surv(time,status) ~1, lung, subset=(sex==2))
|
||
|
|
|
||
|
|
fit1
|
||
|
|
fit2
|
||
|
|
aeq(fit2$scale, c(fit3a$scale, fit3b$scale), tolerance=tol)
|
||
|
|
aeq(fit2$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol)
|
||
|
|
aeq(fit2$coef[1] + 1:2*fit2$coef[2], c(fit3a$coef, fit3b$coef), tolerance=tol)
|
||
|
|
|
||
|
|
#penalized models
|
||
|
|
fit1 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+
|
||
|
|
strata(sex), lung)
|
||
|
|
fit2 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+
|
||
|
|
strata(sex) + sex, lung)
|
||
|
|
fit1
|
||
|
|
fit2
|
||
|
|
|
||
|
|
age1 <- ifelse(lung$sex==1, lung$age, mean(lung$age))
|
||
|
|
age2 <- ifelse(lung$sex==2, lung$age, mean(lung$age))
|
||
|
|
fit3 <- survreg(Surv(time,status) ~ pspline(age1, theta=.92) +
|
||
|
|
pspline(age2, theta=.95) + sex + strata(sex), lung)
|
||
|
|
fit3a<- survreg(Surv(time,status) ~pspline(age, theta=.92), lung,
|
||
|
|
subset=(sex==1))
|
||
|
|
fit3b<- survreg(Surv(time,status) ~pspline(age, theta=.95), lung,
|
||
|
|
subset=(sex==2))
|
||
|
|
fit3b<- survreg(Surv(time,status) ~pspline(age, theta=.95),
|
||
|
|
lung[lung$sex==2,], x=T)
|
||
|
|
#
|
||
|
|
# The above line is tricky, and it took me a long time to realize
|
||
|
|
# it's necessity. The range of age1 = range(age) = 39-82. That for
|
||
|
|
# age2 = range of females = 41-77. The basis functions for pspline are
|
||
|
|
# based on age. If I used data=lung, subset=(sex==2) in fit3b (earlier
|
||
|
|
# form of the test, the pspline function is called before the subset
|
||
|
|
# occurs, and fit3b has a different basis for the second spline than
|
||
|
|
# fit3 does; leading to failure of the all.equal tests below. A theta
|
||
|
|
# of .95 on one basis is not exactly the same as a theta of .95 on the
|
||
|
|
# other. Coefficients were within 1%, but not the same.
|
||
|
|
|
||
|
|
aeq(fit3$scale, c(fit3a$scale, fit3b$scale))
|
||
|
|
aeq(fit3$loglik[2], (fit3a$loglik + fit3b$loglik)[2])
|
||
|
|
pred <- predict(fit3)
|
||
|
|
aeq(pred[lung$sex==1] , predict(fit3a))
|
||
|
|
aeq(pred[lung$sex==2], predict(fit3b))
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|