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.