92 lines
3.5 KiB
Plaintext
92 lines
3.5 KiB
Plaintext
|
|
R Under development (unstable) (2019-06-28 r76752) -- "Unsuffered Consequences"
|
|
Copyright (C) 2019 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.
|
|
|
|
> # A short test on coxph.detail, to ensure that the computed hazard is
|
|
> # equal to the theoretical value
|
|
> library(survival)
|
|
> aeq <- function(a,b) all.equal(as.vector(a), as.vector(b))
|
|
>
|
|
> # taken from book4.R
|
|
> test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
|
|
+ stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
|
|
+ event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
|
|
+ x =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )
|
|
>
|
|
> byhand <- function(beta, newx=0) {
|
|
+ r <- exp(beta)
|
|
+ loglik <- 4*beta - (log(r+1) + log(r+2) + 2*log(3*r+2) + 2*log(3*r+1) +
|
|
+ log(2*r +2))
|
|
+ u <- 1/(r+1) + 1/(3*r+1) + 2*(1/(3*r+2) + 1/(2*r+2)) -
|
|
+ ( r/(r+2) +3*r/(3*r+2) + 3*r/(3*r+1))
|
|
+ imat <- r*(1/(r+1)^2 + 2/(r+2)^2 + 6/(3*r+2)^2 +
|
|
+ 6/(3*r+1)^2 + 6/(3*r+2)^2 + 4/(2*r +2)^2)
|
|
+
|
|
+ hazard <-c( 1/(r+1), 1/(r+2), 1/(3*r+2), 1/(3*r+1), 1/(3*r+1),
|
|
+ 1/(3*r+2), 1/(2*r +2) )
|
|
+
|
|
+
|
|
+ # The matrix of weights, one row per obs, one col per time
|
|
+ # deaths at 2,3,6,7,8,9
|
|
+ wtmat <- matrix(c(1,0,0,0,1, 0, 0,0,0,0,
|
|
+ 0,1,0,1,1, 0, 0,0,0,0,
|
|
+ 0,0,1,1,1, 0, 1,1,0,0,
|
|
+ 0,0,0,1,1, 0, 1,1,0,0,
|
|
+ 0,0,0,0,1, 1, 1,1,0,0,
|
|
+ 0,0,0,0,0, 1, 1,1,1,1,
|
|
+ 0,0,0,0,0,.5,.5,1,1,1), ncol=7)
|
|
+ wtmat <- diag(c(r,1,1,r,1,r,r,r,1,1)) %*% wtmat
|
|
+
|
|
+ x <- c(1,0,0,1,0,1,1,1,0,0)
|
|
+ status <- c(1,1,1,1,1,1,1,0,0,0)
|
|
+ xbar <- colSums(wtmat*x)/ colSums(wtmat)
|
|
+ n <- length(x)
|
|
+
|
|
+ # Table of sums for score and Schoenfeld resids
|
|
+ hazmat <- wtmat %*% diag(hazard) #each subject's hazard over time
|
|
+ dM <- -hazmat #Expected part
|
|
+ for (i in 1:5) dM[i,i] <- dM[i,i] +1 #observed
|
|
+ dM[6:7,6:7] <- dM[6:7,6:7] +.5 # observed
|
|
+ mart <- rowSums(dM)
|
|
+
|
|
+ # Table of sums for score and Schoenfeld resids
|
|
+ # Looks like the last table of appendix E.2.1 of the book
|
|
+ resid <- dM * outer(x, xbar, '-')
|
|
+ score <- rowSums(resid)
|
|
+ scho <- colSums(resid)
|
|
+
|
|
+ # We need to add the ties back up (they are symmetric)
|
|
+ scho[6:7] <- rep(mean(scho[6:7]), 2)
|
|
+
|
|
+ list(loglik=loglik, u=u, imat=imat, xbar=xbar, haz=hazard* exp(beta*newx),
|
|
+ mart=mart, score=score, rmat=resid,
|
|
+ scho=scho)
|
|
+ }
|
|
>
|
|
> # The actual coefficient of the fit is close to zero. Using a larger
|
|
> # number pushes the test harder, but it should still work without
|
|
> # the init and iter arguments, i.e., for any coefficient.
|
|
> fit1 <- coxph(Surv(start, stop, event) ~x, test2,init=-1, iter=0)
|
|
> temp <- coxph.detail(fit1)
|
|
> temp2 <- byhand(fit1$coef, fit1$means)
|
|
> aeq(temp$haz, c(temp2$haz[1:5], sum(temp2$haz[6:7])))
|
|
[1] TRUE
|
|
>
|
|
>
|
|
> proc.time()
|
|
user system elapsed
|
|
1.633 0.128 1.747
|