R Under development (unstable) (2023-08-12 r84939 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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) > options(na.action=na.exclude) # preserve missings > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type > > # Tests of the weighted Cox model > # This is section 1.3 of my appendix -- not yet found in the book > # though, it awaits the next edition > # > # Similar data set to test1, but add weights, > # a double-death/censor tied time > # a censored last subject > # The latter two are cases covered only feebly elsewhere. > # > # The data set testw2 has the same data, but done via replication > # > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y)) > > testw1 <- data.frame(time= c(1,1,2,2,2,2,3,4,5), + status= c(1,0,1,1,1,0,0,1,0), + x= c(2,0,1,1,0,1,0,1,0), + wt = c(1,2,3,4,3,2,1,2,1), + id = 1:9) > # Expanded data set > testw2 <- testw1[rep(1:9, testw1$wt), -4] > row.names(testw2) <- NULL > indx <- match(1:9, testw2$id) > > # Breslow estimate > byhand <- function(beta, newx=0) { + r <- exp(beta) + loglik <- 11*beta - (log(r^2 + 11*r +7) + 10*log(11*r +5) +2*log(2*r+1)) + hazard <- c(1/(r^2 + 11*r +7), 10/(11*r +5), 2/(2*r+1)) + xbar <- c((2*r^2 + 11*r)*hazard[1], 11*r/(11*r +5), r*hazard[3]) + U <- 11- (xbar[1] + 10*xbar[2] + 2*xbar[3]) + imat <- (4*r^2 + 11*r)*hazard[1] - xbar[1]^2 + + 10*(xbar[2] - xbar[2]^2) + 2*(xbar[3] - xbar[3]^2) + + temp <- cumsum(hazard) + risk <- c(r^2, 1,r,r,1,r,1,r,1) + expected <- risk* temp[c(1,1,2,2,2,2,2,3,3)] + + # The matrix of weights, one row per obs, one col per death + # deaths at 1,2,2,2, and 4 + riskmat <- matrix(c(1,1,1,1,1,1,1,1,1, + 0,0,1,1,1,1,1,1,1, + 0,0,1,1,1,1,1,1,1, + 0,0,1,1,1,1,1,1,1, + 0,0,0,0,0,0,0,1,1), ncol=5) + wtmat <- diag(c(r^2, 2, 3*r, 4*r, 3, 2*r, 1, 2*r, 1)) %*% riskmat + + x <- c(2,0,1,1,0,1,0,1,0) + status <- c(1,0,1,1,1,0,0,1,0) + wt <- c(1,2,3,4,3,2,1,2,1) + # Table of sums for score and Schoenfeld resids + hazmat <- riskmat %*% diag(c(1,3,4,3,2)/colSums(wtmat)) + dM <- -risk*hazmat #Expected part + dM[1,1] <- dM[1,1] +1 # deaths at time 1 + for (i in 2:4) dM[i+1, i] <- dM[i+1,i] +1 + dM[8,5] <- dM[8,5] +1 + mart <- rowSums(dM) + resid <-dM * outer(x, xbar[c(1,2,2,2,3)] ,'-') + + # Increments to the variance of the hazard + var.g <- cumsum(hazard^2/ c(1,10,2)) + var.d <- cumsum((xbar-newx)*hazard) + + list(loglik=loglik, U=U, imat=imat, hazard=hazard, xbar=xbar, + mart=c(1,0,1,1,1,0,0,1,0)-expected, expected=expected, + score=rowSums(resid), schoen=c(2,1,1,0,1) - xbar[c(1,2,2,2,3)], + varhaz=(var.g + var.d^2/imat)* exp(2*beta*newx)) + } > > aeq(byhand(0)$expected, c(1/19, 1/19, rep(103/152, 5), rep(613/456,2))) #verify [1] TRUE > > fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt, + method='breslow', iter=0) > fit0b <- coxph(Surv(time, status) ~x, testw2, method='breslow', iter=0) > fit <- coxph(Surv(time, status) ~x, testw1, weights=wt, method='breslow') > fitb <- coxph(Surv(time, status) ~x, testw2, method='breslow') > > aeq(resid(fit0, type='mart'), (resid(fit0b, type='mart'))[indx]) [1] TRUE > aeq(resid(fit0, type='scor'), (resid(fit0b, type='scor'))[indx]) [1] TRUE > aeq(unique(resid(fit0, type='scho')), unique(resid(fit0b, type='scho'))) [1] TRUE > > truth0 <- byhand(0,pi) > aeq(fit0$loglik[1], truth0$loglik) [1] TRUE > aeq(1/truth0$imat, fit0$var) [1] TRUE > aeq(truth0$mart, fit0$resid) [1] TRUE > aeq(truth0$scho, resid(fit0, 'schoen')) [1] TRUE > aeq(truth0$score, resid(fit0, 'score')) [1] TRUE > sfit <- survfit(fit0, list(x=pi), censor=FALSE) > aeq(sfit$std.err^2, truth0$var) [1] TRUE > aeq(-log(sfit$surv), cumsum(truth0$haz)) [1] TRUE > > truth <- byhand(0.85955744, .3) > aeq(truth$loglik, fit$loglik[2]) [1] TRUE > aeq(1/truth$imat, fit$var) [1] TRUE > aeq(truth$mart, fit$resid) [1] TRUE > aeq(truth$scho, resid(fit, 'schoen')) [1] TRUE > aeq(truth$score, resid(fit, 'score')) [1] TRUE > > sfit <- survfit(fit, list(x=.3), censor=FALSE) > aeq(sfit$std.err^2, truth$var) [1] TRUE > aeq(-log(sfit$surv), (cumsum(truth$haz)* exp(fit$coef*.3))) [1] TRUE > > > fit0 Call: coxph(formula = Surv(time, status) ~ x, data = testw1, weights = wt, method = "breslow", iter = 0) coef exp(coef) se(coef) z p x 0.0000 1.0000 0.5858 0 1 Likelihood ratio test=0 on 1 df, p=1 n= 9, number of events= 5 > summary(fit) Call: coxph(formula = Surv(time, status) ~ x, data = testw1, weights = wt, method = "breslow") n= 9, number of events= 5 coef exp(coef) se(coef) z Pr(>|z|) x 0.8596 2.3621 0.7131 1.205 0.228 exp(coef) exp(-coef) lower .95 upper .95 x 2.362 0.4233 0.5839 9.556 Concordance= 0.637 (se = 0.161 ) Likelihood ratio test= 1.69 on 1 df, p=0.2 Wald test = 1.45 on 1 df, p=0.2 Score (logrank) test = 1.52 on 1 df, p=0.2 > resid(fit0, type='score') 1 2 3 4 5 6 1.24653740 0.03601108 0.10056700 0.10056700 -0.22180142 -0.21193300 7 8 9 0.46569858 -0.10082189 0.91014302 > resid(fit0, type='scho') 1 2 2 2 4 1.3157895 0.3125000 0.3125000 -0.6875000 0.3333333 > > resid(fit, type='score') 1 2 3 4 5 6 0.88681615 0.02497653 0.03608964 0.03608964 -0.54297652 -0.12528780 7 8 9 0.29564605 -0.09476911 0.58400064 > resid(fit, type='scho') 1 2 2 2 4 1.0368337 0.1613774 0.1613774 -0.8386226 0.1746960 > aeq(resid(fit, type='mart'), (resid(fitb, type='mart'))[indx]) [1] TRUE > aeq(resid(fit, type='scor'), (resid(fitb, type='scor'))[indx]) [1] TRUE > aeq(unique(resid(fit, type='scho')), unique(resid(fitb, type='scho'))) [1] TRUE > rr1 <- resid(fit, type='mart') > rr2 <- resid(fit, type='mart', weighted=T) > aeq(rr2/rr1, testw1$wt) [1] TRUE > > rr1 <- resid(fit, type='score') > rr2 <- resid(fit, type='score', weighted=T) > aeq(rr2/rr1, testw1$wt) [1] TRUE > > > proc.time() user system elapsed 1.42 0.15 1.56