[R] help-- Cox ph models

[John Maindonald]

The following code should do the job.  You can if you want rename 
sum.coxph to
summary.coxph.  print.summary.coxph is needed so that printing of the 
objects
that are thus created has the same effect as calling the existing 
summary.coxph
It needs a modified help file to go with it.  I have been meaning to 
try to get this,
or something like it, incorporated into the survival package.

"sum.coxph" <-
function (object, table = TRUE, coef = TRUE, conf.int = 0.95,
     scale = 1,  ...)
{
     cox <- object
coxsum <- list(call=cox$call, fail=cox$fail, na.action=cox$na.action,
   n=cox$n, icc=cox$icc, naive.var=cox$naive.var)
     if (!is.null(cox$fail)) {
         return()
     }
     omit <- cox$na.action
     if (is.null(cox$coef)) {
         return()
     }
     beta <- cox$coef
     nabeta <- !(is.na(beta))
     beta2 <- beta[nabeta]
     if (is.null(beta) | is.null(cox$var)){
         coxsum$invalid <- TRUE
         return()
         }
         else coxsum$invalid<-FALSE
     se <- sqrt(diag(cox$var))
     if (!is.null(cox$naive.var))
         nse <- sqrt(diag(cox$naive.var))
     if (coef) {
         if (is.null(cox$naive.var)) {
             tmp <- cbind(beta, exp(beta), se, beta/se, 1 -
                 pchisq((beta/se)^2, 1))
             dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)",
                 "se(coef)", "z", "p"))
         }
         else {
             tmp <- cbind(beta, exp(beta), nse, se, beta/se, 1 -
                 pchisq((beta/se)^2, 1))
             dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)",
                 "se(coef)", "robust se", "z", "p"))
         }
     } else tmp <- NULL
     if (conf.int) {
         z <- qnorm((1 + conf.int)/2, 0, 1)
         beta <- beta * scale
         se <- se * scale
         tmpc <- cbind(exp(beta), exp(-beta), exp(beta - z * se),
             exp(beta + z * se))
         dimnames(tmpc) <- list(names(beta), c("exp(coef)", "exp(-coef)",
             paste("lower .", round(100 * conf.int, 2), sep = ""),
             paste("upper .", round(100 * conf.int, 2), sep = "")))
     } else tmpc <- NULL
     logtest <- -2 * (cox$loglik[1] - cox$loglik[2])
     sctest <- cox$score
     df <- length(beta2)
     Rsquare <- 1 - exp(-logtest/cox$n)
     maxRsquare <- 1 - exp(2 * cox$loglik[1]/cox$n)
     plogtest <- 1 - pchisq(logtest, df)
     pwald.test <- 1 - pchisq(cox$wald.test, df)
     pscore <- 1 - pchisq(sctest, df)
     if(!is.null(cox$rscore)){
       rscore <- cox$rscore
       prscore <- 1 - pchisq(cox$rscore, df)
       robust <- list(val=rscore, pval=prscore)
       }
     else robust <- NULL
     wald.test <- cox$wald.test
     coxsum <- c(list(call=cox$call, coefficients=tmp, intervals=tmpc, 
df=df,
       rsquare=c(Rsquare, maxRsquare),
       tests=c(logtest=logtest, wald.test=wald.test, score=sctest),
       p.tests = c(logtest=plogtest, wald.test=pwald.test,
         score=pscore), robust=robust), coxsum)
     class(coxsum) <- "summary.coxph"
     coxsum
}
"print.summary.coxph" <-
function (object, table = TRUE, coef = TRUE, conf.int = 0.95,
     scale = 1, digits = max(options()$digits - 4, 3), ...)
{
     coxsum <- object
     if (!is.null(cl <- coxsum$call)) {
         cat("Call:\n")
         dput(cl)
         cat("\n")
     }
     if (!is.null(coxsum$fail)) {
         cat(" Coxreg failed.", coxsum$fail, "\n")
         return()
     }
     savedig <- options(digits = digits)
     on.exit(options(savedig))
     omit <- coxsum$na.action
     if (length(omit))
         cat("  n=", coxsum$n, " (", naprint(omit), ")\n", sep = "")
     else cat("  n=", coxsum$n, "\n")
     if (length(coxsum$icc))
         cat("  robust variance based on", coxsum$icc[1], "groups, 
intra-class correlation =",
             format(coxsum$icc[2:3]), "\n")
     if (is.null(coxsum$coef)) {
         cat("   Null model\n")
         return()
     }
     if (coxsum$invalid)stop("input is invalid")
       if (!is.null(coxsum$coef)) {
         cat("\n")
         coxsum$coef[,"p"] <- signif(coxsum$coef[,"p"],digits-1)
         prmatrix(coxsum$coef)
     }
     if (!is.null(coxsum$intervals)) {
         cat("\n")
         prmatrix(coxsum$intervals)
     }
     logtest <- coxsum$tests["logtest"]
     sctest <- coxsum$tests["score"]
     df <- coxsum$df
     cat("\n")
     cat("Rsquare=", format(round(coxsum$rsquare[1], 3)),
         "  (max possible=", format(round(coxsum$rsquare[2],3)), ")\n")
     cat("Likelihood ratio test= ", 
format(round(coxsum$tests["logtest"], 2)),
         "  on ", df, " df,", "   p=", format(coxsum$p.tests["logtest"]),
         "\n", sep = "")
     cat("Wald test            = ", format(coxsum$tests["wald.test.x"]),
         "  on ", df, " df,", "   p=", format(
         coxsum$p.tests["wald.test.x"]), "\n", sep = "")
     cat("Score (logrank) test = ", format(round(coxsum$tests["score"], 
2)),
         "  on ", df, " df,", "   p=", format(coxsum$p.tests["score"]),
         sep = "")
     if (is.null(coxsum$robust))
         cat("\n\n")
     else cat(",   Robust = ", format(round(coxsum$robust$val, 2)), "  
p=",
         format(coxsum$robust$pval), "\n\n", sep = "")
     if (!is.null(coxsum$tests))
         cat("  (Note: the likelihood ratio and score tests",
             "assume independence of\n     observations within a 
cluster,",
             "the Wald and robust score tests do not).\n")
    invisible()
}

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Bioinformation Science, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.