packrat/lib-R/survival/tests/detail.R

In UBC-MDS/Karl:

# 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])))