Last active
September 2, 2015 14:58
-
-
Save carlbfrederick/c351c84ee9fca9888c57 to your computer and use it in GitHub Desktop.
I wrote these functions to calculate survival probabilities at various levels of the random effect estimate from a coxme.object. The code is heavily adapted from survfit.coxph(). It should work, but all errors and inelegant hacks I have introduced are certainly my own doing. Comments/improvements welcome, enjoy!
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#Internal Functions | |
MYagsurv <- function(y, x, wt, risk, survtype=3, vartype=3) { | |
nvar <- ncol(as.matrix(x)) | |
status <- y[, ncol(y)] | |
dtime <- y[, ncol(y) - 1] | |
death <- (status == 1) | |
time <- sort(unique(dtime)) | |
nevent <- as.vector(rowsum(wt * death, dtime)) | |
ncens <- as.vector(rowsum(wt * (!death), dtime)) | |
wrisk <- c(wt * risk) #Had to add c() to remove the dimnames so that the multiplication later would work | |
rcumsum <- function(x) rev(cumsum(rev(x))) | |
nrisk <- rcumsum(rowsum(wrisk, dtime)) | |
irisk <- rcumsum(rowsum(wt, dtime)) | |
if (ncol(y) == 2) { | |
temp2 <- rowsum(wrisk * x, dtime) | |
xsum <- apply(temp2, 2, rcumsum) | |
} else { | |
delta <- min(diff(time))/2 | |
etime <- c(sort(unique(y[, 1])), max(y[, 1]) + delta) | |
indx <- approx(etime, 1:length(etime), time, method = "constant", | |
rule = 2, f = 1)$y | |
esum <- rcumsum(rowsum(wrisk, y[, 1])) | |
nrisk <- nrisk - c(esum, 0)[indx] | |
irisk <- irisk - c(rcumsum(rowsum(wt, y[, 1])), 0)[indx] | |
xout <- apply(rowsum(wrisk * x, y[, 1]), 2, rcumsum) | |
xin <- apply(rowsum(wrisk * x, dtime), 2, rcumsum) | |
xsum <- xin - (rbind(xout, 0))[indx, , drop = F] | |
} | |
#Calculate number of deaths and number of times | |
ndeath <- rowsum(status, dtime) | |
ntime <- length(time) | |
#Final calcs, send all of this stuff to the proper C subroutine | |
if (survtype == 1) { | |
indx <- (which(status == 1))[order(dtime[status == 1])] | |
km <- .C(survival:::Cagsurv4, as.integer(ndeath), as.double(risk[indx]), | |
as.double(wt[indx]), as.integer(ntime), as.double(nrisk), | |
inc = double(ntime)) | |
} | |
if (survtype == 3 || vartype == 3) { | |
xsum2 <- rowsum((wrisk * death) * x, dtime) | |
erisk <- rowsum(wrisk * death, dtime) | |
tsum <- .C(survival:::Cagsurv5, as.integer(length(nevent)), as.integer(nvar), | |
as.integer(ndeath), as.double(nrisk), as.double(erisk), | |
as.double(xsum), as.double(xsum2), sum1 = double(length(nevent)), | |
sum2 = double(length(nevent)), xbar = matrix(0, | |
length(nevent), nvar)) | |
} | |
#Calc hazard rates, etc and wrap up results | |
haz <- switch(survtype, nevent/nrisk, nevent/nrisk, nevent * | |
tsum$sum1) | |
varhaz <- switch(vartype, nevent/(nrisk * ifelse(nevent >= | |
nrisk, nrisk, nrisk - nevent)), nevent/nrisk^2, nevent * | |
tsum$sum2) | |
xbar <- switch(vartype, (xsum/nrisk) * haz, (xsum/nrisk) * | |
haz, nevent * tsum$xbar) | |
result <- list(n = nrow(y), time = time, n.event = nevent, | |
n.risk = irisk, n.censor = ncens, hazard = haz, cumhaz = cumsum(haz), | |
varhaz = varhaz, ndeath = ndeath, xbar = apply(matrix(xbar, | |
ncol = nvar), 2, cumsum)) | |
if (survtype == 1) { | |
result$surv <- cumprod(km$inc) | |
} | |
else { | |
result$surv <- exp(-result$cumhaz) | |
} | |
return(result) | |
} | |
MYrisk <- function(object, newdata, rfx=0) { | |
dMat <- survival:::model.matrix.coxph(object, newdata) | |
xmeans <- object$means | |
names(xmeans) <- names(fixef(object)) | |
dMat <- scale(dMat, center=xmeans, scale=FALSE) | |
newrisk <- exp(dMat %*% fixef(object) + rfx) | |
return(newrisk) | |
} | |
predSurv <- function(survlist, newrisk) { | |
n.survlist <- length(survlist) | |
strata.names <- names(survlist) | |
n.newrisk <- nrow(newrisk) | |
for (i in 1:n.survlist) { | |
survlist[[i]]$surv <- as.data.frame(outer(survlist[[i]]$surv, newrisk, "^")) | |
colnames(survlist[[i]]$surv) <- paste("Surv_newdata_row", 1:n.newrisk, sep="") | |
survlist[[i]]$cumhaz <- as.data.frame(outer(survlist[[i]]$cumhaz, newrisk, "*")) | |
colnames(survlist[[i]]$cumhaz) <- paste("CumHaz_newdata_row", 1:n.newrisk, sep="") | |
survlist[[i]]$newrisk <- newrisk | |
dimnames(survlist[[i]]$newrisk)[[1]] <- list(paste("newdata_row", 1:n.newrisk, sep="")) | |
} | |
return(survlist) | |
} | |
survProbs <- function(survHat, newtime) { | |
out <- NULL | |
for (i in 1:length(survHat)) { | |
if (any(newtime == survHat[[i]]$time)) { | |
idx <- survHat[[i]]$time==newtime | |
out[[i]] <- data.frame( | |
Time=survHat[[i]]$time[idx], | |
`Survival Probability`=survHat[[i]]$surv[idx,1], | |
`Cumulative Incidence`=1-survHat[[i]]$surv[idx,1] | |
) | |
} else{ | |
btwnrows <- c(max(survHat[[i]]$time[survHat[[i]]$time<newtime]), | |
min(survHat[[i]]$time[survHat[[i]]$time>newtime])) | |
idx <- survHat[[i]]$time %in% btwnrows | |
interpDF <- data.frame( | |
time=survHat[[i]]$time[idx], | |
surv=survHat[[i]]$surv[idx,1] | |
) | |
out[[i]] <- data.frame( | |
Time=newtime, | |
`Survival Probability`=approx(x=interpDF[,"time"], y=interpDF[,"surv"], xout=newtime)$y | |
) | |
out[[i]]$`Cumulative Incidence`=1-out[[i]][,2] | |
} | |
} | |
return(out) | |
} | |
#Final Function | |
survProbs.coxme <- function(object, stratVar=NULL, rfx=0, groupFctr=NULL, level=NULL, newdata, newtime) { | |
#prep values for MYagsurv() | |
y <- object$y | |
x <- object$x | |
if (is.null(object$weights)) { | |
wt <- rep(1, nrow(y)) | |
} else wt <- object$weights | |
survtype <- vartype <- switch(object$ties, efron=3, breslow=2) | |
risk <- exp(x %*% fixef(object)) | |
#do MYagsurv() by strata | |
if (is.null(stratVar)) { | |
stratvar <- wt | |
} | |
if (is.factor(stratVar)) { | |
ustrata <- levels(stratVar) | |
} else ustrata <- sort(unique(stratVar)) | |
nstrata <- length(ustrata) | |
survlist <- vector("list", nstrata) | |
names(survlist) <- ustrata | |
for (i in 1:nstrata) { | |
indx <- which(stratVar == ustrata[i]) | |
survlist[[i]] <- MYagsurv(y[indx, , drop = F], x[indx, , drop = F], wt[indx], risk[indx], survtype, vartype) | |
} | |
#Calc newrisk for user-supplied newdata | |
if (!is.null(level)) { | |
if (!is.null(groupFctr)) { | |
rfx <- as.numeric(ranef(object)[[groupFctr]][level]) | |
} else rfx <- as.numeric(ranef(object)[[1]][level]) | |
} | |
newrisk <- MYrisk(object, newdata, rfx) | |
#Calculate survival function based on covariate profile in newdata | |
survHat <- predSurv(survlist, newrisk) | |
#Find value of survival functionat newtime | |
survProbs <- survProbs(survHat, newtime) | |
names(survProbs) <- names(survlist) | |
print(survProbs) | |
return(list(survList=survlist, survListHat=survHat, survProbs=survProbs)) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment