R version 4.4.0 beta (2024-04-12 r86413 ucrt) -- "Puppy Cup" Copyright (C) 2024 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. > library(survival) > # > # Tests with the pspline function, to verify the prediction aspects > # > options(na.action=na.exclude) > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...) > > spfit <- coxph(Surv(time, status) ~ pspline(age) + ph.ecog, lung) > > spfit2 <- coxph(Surv(time, status) ~ pspline(age) + ph.ecog, lung, x=TRUE) > x2 <- model.matrix(spfit) > all.equal(spfit2$x, x2) [1] TRUE > > keep <- (lung$age < 60) > x3 <- model.matrix(spfit, data=lung[keep,]) > attr(x3, 'assign') <- NULL #subscripting loses the assign attr below > all.equal(napredict(spfit$na.action,x2)[keep,], x3) [1] TRUE > > p2 <- predict(spfit, newdata=lung[keep,]) > aeq(p2, predict(spfit)[keep]) [1] TRUE > > > p3 <- survfit(spfit) > p4 <- survfit(spfit, newdata=lung[1:2,]) > temp <- scale(x2[1:2,], center=spfit$means, scale=FALSE)%*% coef(spfit) > aeq(p3$time, p4$time) [1] TRUE > aeq(outer(-log(p3$surv), exp(temp), '*'), -log(p4$surv)) [1] TRUE > > # Check out model.frame > spfit3 <- coxph(Surv(time, status) ~ pspline(age) + sex, lung, + model=TRUE) #avoid the missing value > m2 <- model.frame(spfit3, data=lung[keep,]) > all.equal(m2, spfit3$model[keep,], check.attributes=FALSE) [1] TRUE > > # > # Test of residuals, in response to a reported bug. > # These are three progam paths that should all lead to the same C routine > fit <- coxph(Surv(tstart, tstop, status) ~ sex + treat + pspline(age), cgd) > fit2 <- coxph(Surv(tstart, tstop, status) ~ fit$linear, cgd, iter=0, init=1) > fit3 <- coxph(Surv(tstart, tstop, status) ~ offset(fit$linear), cgd) > all.equal(fit$resid, fit2$resid) [1] TRUE > all.equal(fit$resid, fit3$resid) [1] TRUE > > # > # Check using coxph.detail. The matrix multiply below only is > # valid for the breslow approximation. > fit4 <- coxph(Surv(tstart, tstop, status) ~ sex + treat + pspline(age), + cgd, ties='breslow') > dt <- coxph.detail(fit4, riskmat=TRUE) > > # the results of coxph.detail used to be in time order, now are in data set > # order > rscore <- exp(fit4$linear) > exp4 <- (rscore *dt$riskmat) %*% dt$hazard > r4 <- cgd$status - exp4 > aeq(r4, fit4$resid) [1] TRUE > > proc.time() user system elapsed 0.81 0.14 0.93