File: C:/Program Files/R/R-4.4.0/library/survival/tests/book6.Rout.save
R Under development (unstable) (2020-02-28 r77869) -- "Unsuffered Consequences"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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 -- no yet found in any of the
> # printings though, it awaits the next edition
> #
> # Efron approximation
> #
> 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))
> xx <- testw1$wt
>
> # Efron estimate
> byhand <- function(beta, newx=0) {
+ r <- exp(beta)
+ a <- 7*r +3; b<- 4*r+2
+ loglik <- 11*beta - (log(r^2 + 11*r +7) + 10*log(11*r +5)/3 +
+ 10*log(a*2/3 +b)/3 + 10*log(a/3 +b)/3 +2*log(2*r+1))
+
+ hazard <- c(1/(r^2 + 11*r +7),
+ 10/(3*c(11*r +5, a*2/3 +b, a/3+b)), 2/(2*r+1))
+ temp <- c(hazard[1], hazard[1]+hazard[2] + hazard[3]*2/3 + hazard[4]/3,
+ cumsum(hazard)[4:5])
+ risk <- c(r^2, 1,r,r,1,r,1,r,1)
+ expected <- risk* temp[c(1,1,2,2,2,3,3,4,4)]
+
+ # 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,2/3,2/3,2/3,1,1,1,1,
+ 0,0,1/3,1/3,1/3,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)
+ xbar <- colSums(x*wtmat)/ colSums(wtmat)
+ imat <- (4*r^2 + 11*r)*hazard[1] - xbar[1]^2 +
+ 10* mean(xbar[2:4] - xbar[2:4]^2) + 2*(xbar[5] - xbar[5]^2)
+
+ 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 resids
+ hazmat <- riskmat %*% diag(c(1,10/3,10/3, 10/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[3:5, i] <- dM[3:5,i] + 1/3
+ dM[8,5] <- dM[8,5] +1
+ mart <- rowSums(dM)
+ resid <-dM * outer(x, xbar ,'-')
+
+ # Increments to the variance of the hazard
+ var.g <- cumsum(hazard^2* c(1,3/10, 3/10, 3/10, 1/2))
+ var.d <- cumsum((xbar-newx)*hazard)
+
+ sxbar <- c(xbar[1], mean(xbar[2:4]), xbar[5]) #xbar for Schoen
+ list(loglik=loglik, imat=imat, hazard=hazard, xbar=xbar,
+ mart=status-expected, expected=expected,
+ score=rowSums(resid), schoen=c(2,1,1,0,1) - sxbar[c(1,2,2,2,3)],
+ varhaz=((var.g + var.d^2/imat)* exp(2*beta*newx))[c(1,4,5)])
+ }
>
> # Verify
> temp <- byhand(0,0)
> aeq(temp$xbar, c(13/19, 11/16, 26/38, 19/28, 2/3))
[1] TRUE
> aeq(temp$hazard, c(1/19, 5/24, 5/19, 5/14, 2/3))
[1] TRUE
>
> fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt, iter=0)
> fit <- coxph(Surv(time, status) ~x, testw1, weights=wt)
>
> 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$hazard)[c(1,4,5)])
[1] TRUE
>
> truth <- byhand(fit$coef, .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$hazard)* exp(fit$coef*.3))[c(1,4,5)])
[1] TRUE
>
>
> fit0
Call:
coxph(formula = Surv(time, status) ~ x, data = testw1, weights = wt,
iter = 0)
coef exp(coef) se(coef) z p
x 0.0000 1.0000 0.5843 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)
n= 9, number of events= 5
coef exp(coef) se(coef) z Pr(>|z|)
x 0.8726 2.3931 0.7126 1.225 0.221
exp(coef) exp(-coef) lower .95 upper .95
x 2.393 0.4179 0.5921 9.672
Concordance= 0.637 (se = 0.161 )
Likelihood ratio test= 1.75 on 1 df, p=0.2
Wald test = 1.5 on 1 df, p=0.2
Score (logrank) test = 1.58 on 1 df, p=0.2
> resid(fit0, type='score')
1 2 3 4 5 6
1.24653740 0.03601108 0.14118105 0.14118105 -0.30336782 -0.27962308
7 8 9
0.60164259 -0.16851197 1.04608703
> resid(fit0, type='scho')
1 2 2 2 4
1.3157895 0.3165727 0.3165727 -0.6834273 0.3333333
>
> resid(fit, type='score')
1 2 3 4 5 6
0.88116056 0.02477248 0.06057806 0.06057806 -0.59724033 -0.16737066
7 8 9
0.38040295 -0.13750290 0.66631324
> resid(fit, type='scho')
1 2 2 2 4
1.0325955 0.1621759 0.1621759 -0.8378241 0.1728229
>
> 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
0.891 0.068 0.952