avg.surv <- function(cfit, var.name, var.values, data, weights) { if(missing(data)) { if(!is.null(cfit$model)) mframe <- cfit$model else mframe <- model.frame(cfit, sys.parent()) } else mframe <- model.frame(cfit, data) var.num <- match(var.name, names(mframe)) data.patterns <- apply(data.matrix(mframe[, - c(1, var.num)]), 1, paste, collapse = ",") data.patterns <- factor(data.patterns,levels=unique(data.patterns)) if(missing(weights)) weights <- table(data.patterns) else weights <- tapply(weights, data.patterns, sum) kp <- !duplicated(data.patterns) mframe <- mframe[kp,] obs.var <- mframe[,var.num] lps <- (cfit$linear.predictor)[kp] tframe <- mframe[rep(1,length(var.values)),] tframe[,var.num] <- var.values xmat <- model.matrix(cfit,data=tframe)[,-1] tlps <- as.vector(xmat%*%cfit$coef) names(tlps) <- var.values obs.off <- tlps[as.character(obs.var)] explp.off <- exp(outer(lps - obs.off ,tlps,"+")) bfit <- survfit.coxph(cfit, se.fit = F) fits <- outer(bfit$surv,explp.off,function(x,y) x^y) avg.fits <- apply(sweep(fits,2,weights,"*"),c(1,3),sum)/sum(weights) dimnames(avg.fits) <- list(NULL,var.values) list(time=bfit$time,fits=avg.fits) }