### stat functions
# extractAIC.crr, AIC.crr, BIC.crr, extractIC, logLik.crr, deviance.crr,
# coef.crr, coefficients.crr, terms.crr2, terms.crr, crrFits, crrWald.test
#
# unexported: wald.test
###
#' \code{crr} model fit statistics
#'
#' @description
#' Functions to assess the quality of fitted \code{\link[cmprsk]{crr}}
#' and \code{\link{crr2}} objects.
#'
#' Select multiple types of Akaike or Bayesian (Schwarz) information
#' criterion from a \code{\link[cmprsk]{crr}} object to assess the relative
#' quality of models for a given data set.
#'
#' @details
#' \code{AIC} and \code{BIC} are calculated in the usual way. \code{AICc}
#' is the AIC with a correction for finite sample sizes. This assumes that
#' the model is univariate, linear, and has normally-distributed residuals
#' (conditional upon regressors).
#'
#' \code{BICc} is a newly proposed criteria that is a modified BIC for
#' competing risks data subject to right censoring (Kuk, 2013).
#'
#' @param object an object of class \code{\link[cmprsk]{crr}}
#' @param ic information criterion, one of \code{"AIC"}, \code{"BIC"},
#' \code{"AICc"}, or \code{"BICc"}; see details
#' @param p an optional penalty term to be multiplied by, \code{k}, the
#' number of free parameters estimated in each model including the intercept
#' term
#' @param fit,scale,k see \code{\link{extractAIC}}
#' @param ... additional arguments passed to or from other methods
#'
#' @seealso
#' \code{\link{crrFits}}; \code{\link{crrwald.test}}
#'
#' @references
#' Kuk D, Varadhan R. Model selection in competing risks regression.
#' \emph{Stat Med}. 2013 Aug \strong{15};32.
#'
#' @examples
#' crrs <- crr2(Surv(futime, event(censored) == death) ~ age, transplant)
#' crr1 <- crrs[[1L]]
#'
#' crr1$coef
#' coef(crr1)
#' coefficients(crr1)
#'
#' extractAIC(crr1)
#' sapply(crrs, AIC)
#'
#' ## these are equivalent
#' extractIC(crrs[[1L]], p = 0)
#' -2 * logLik(crrs[[1]])[1]
#' -2 * crrs[[1]]$loglik
#'
#' deviance(crr1)
#'
#' crrFits(crr1)
#' crrwald.test(crr1)
#'
#' @name crrfit
NULL
#' @rdname crrfit
#' @export
extractIC <- function(object, ic = c('AIC', 'BIC', 'AICc', 'BICc'), p) {
cl <- object$call
ic <- match.arg(ic)
n <- object[['n']]
k <- length(object[['coef']])
## get number of events of interest from function call
ne <- tryCatch(
sum({object$nuftime} %||% {eval(cl$fstatus) %in% eval(cl$failcode)}),
error = function(e) {
if (ic == 'BIC')
warning(
'BIC is approximate: use p = # events of interest', call. = FALSE
)
length(object[['uftime']])
})
p <- if (!missing(p)) {
ic <- paste('p', p, sep = '=')
p
} else
switch(
ic,
AIC = 2,
BIC = log(ne),
AICc = 2,
BICc = log(length(object[['res']][, 1L]))
)
res <- p * k - 2 * object[['loglik']]
if (ic == 'AICc')
res <- res + 2 * (k + 1) * (k + 2) / (n - k - 2)
setNames(res, ic)
}
#' @rdname crrfit
#' @export
extractAIC.crr <- function(fit, scale, k = 2, ...) {
extractIC(fit, p = k)
}
#' @rdname crrfit
#' @export
AIC.crr <- function(object, ...) {
extractIC(object, 'AIC')
}
#' @rdname crrfit
#' @export
BIC.crr <- function(object, ...) {
extractIC(object, 'BIC')
}
#' @rdname crrfit
#' @export
logLik.crr <- function(object, ...) {
structure(
object[['loglik']], df = length(object[['coef']]),
class = 'logLik'
)
}
#' @rdname crrfit
#' @export
deviance.crr <- function(object, ...) {
cat('Call:\n')
dput(object$call)
cat('\n')
x <- summary(object)
dv <- unname(x$logtest[['test']])
df <- x$logtest[['df']]
pv <- pchisq(dv, df, lower.tail = FALSE)
cat('Deviance =', signif(dv, 3L), 'on', df, 'df,',
format.pval(pv, show.p = TRUE), '\n\n')
invisible(c(chi.sq = dv, df = df, p.value = pv))
}
#' \code{crr} model coefficients
#'
#' @param object an object of class \code{\link[cmprsk]{crr}}
#' @param ... ignored
#'
#' @export
coef.crr <- function(object, ...) {
object[['coef']]
}
#' @rdname coef.crr
#' @export
coefficients.crr <- coef.crr
#' \code{crr2} model terms
#'
#' @param x an object of class \code{\link{crr2}}
#' @param ... ignored
#'
#' @export
terms.crr2 <- function(x, ...) {
if (!inherits(x, 'crr2')) {
message('\'terms\' is not available for a \'crr\' object - try \'?crr2\'')
return(invisible(NULL))
}
terms(attr(x, 'model.frame'))
}
#' @rdname terms.crr2
#' @export
terms.crr <- terms.crr2
#' \code{crr} model selection table
#'
#' Return several types of fit statistics for a \code{\link[cmprsk]{crr}}
#' saturated model compared to the null model. Note that comparisons among
#' models only make sense for ones fit to the same data set.
#'
#' @param ... one or more objects of class \code{\link[cmprsk]{crr}}
#' @param p an optional penalty term to be multiplied by, \code{k}, the
#' number of free parameters estimated in each model including the intercept
#' term
#'
#' @seealso
#' \code{\link{crrfit}}; \code{\link{crrwald.test}}
#'
#' @references
#' Scrucca L, Santucci A, Aversa F (2009). Regression Modeling of Competing
#' Risk Using R: An In Depth Guide for Clinicians. \emph{Bone Marrow
#' Transplantation} (2010) \strong{45}, 1388-1395.
#'
#' @examples
#' \dontrun{
#' ## example, figures, tables from
#' ## http://www.nature.com/bmt/journal/v45/n9/full/bmt2009359a.html
#'
#' bmt <- read.csv('http://www.stat.unipg.it/luca/misc/bmtcrr.csv')
#' bmt <- within(bmt, {
#' Sex <- relevel(Sex, 'M')
#' Phase <- relevel(Phase, 'Relapse')
#' })
#' cov1 <- with(bmt, model.matrix(~ Age + Sex + D + Phase + Source)[, -1])
#'
#' m1 <- with(bmt, crr(ftime, Status, cov1))
#' summary(m1)
#' crrwald.test(m1, c(Phase = 'Phase'))
#'
#' ## model selection
#' m2 <- with(bmt, crr(ftime, Status, cov1[, c(4:6)]))
#' m3 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7)]))
#' m4 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 1)]))
#' m5 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 2)]))
#' m6 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 3)]))
#'
#' crrFits(m1, m2, m3, m4, m5, m6, p = 3)
#'
#'
#' par(mfrow = c(2, 2))
#' with(m2, {
#' for (ii in seq.int(ncol(res)))
#' scatter.smooth(uftime, res[, ii], main = names(coef)[ii],
#' xlab = 'Failure time', ylab = 'Schoenfeld residuals')
#' })
#'
#' pred <- with(bmt, model.matrix(~ levels(Phase)))[, -1L]
#' pred <- predict(m2, pred)
#'
#' plot(
#' pred, col = 1:4, lty = 1, ylim = c(0, 1),
#' xlab = 'Failure time', ylab = 'CIF'
#' )
#' legend(
#' 'top', lty = 1L, col = 1:4, horiz = TRUE, bty = 'n',
#' title = 'Phase', legend = levels(bmt$Phase),
#' x.intersp = 0.1, y.intersp = 0.5
#' )
#' }
#'
#' @export
crrFits <- function(..., p) {
l <- if (islist(..1)) {
if (inherits(..1[[1L]], 'crr2') & is.null(coef(..1[[1L]])))
unlist(c(...), recursive = FALSE) else c(...)
} else {
if (inherits(..1, 'crr2') & is.null(coef(..1)))
c(...) else list(...)
}
## drop cox models from crr2 objects
l <- Filter(function(x) !inherits(x, 'coxph'), l)
stopifnot(
all(sapply(l, inherits, 'crr')),
## assert same data was used to fit
# all(diff(sapply(l, function(x) x[['loglik.null']])) == 0),
all(diff(sapply(l, function(x) x[['n']])) == 0L)
)
null <- l[[1L]]
null$loglik <- null$loglik.null
null$coef <- null$call$cov1 <- null$call$cov2 <- NULL
l <- c(list(null), l)
n <- length(l)
# model <- c('Null model', sapply(l[-1L], function(x) x[['call']]))
model <- sapply(seq_along(l), function(ii)
capture.output({
if (ii == 1L)
cat('0: Null Model', '\n')
else {
cat(sprintf('%s: Model %1$s call:\n', ii - 1L))
dput(l[[ii]][['call']])
}
cat('\n')
}))
res <- list()
res$n <- vapply(l, function(x) x$n, integer(1L))
res$loglik <- vapply(l, logLik, numeric(1L))
res$df <- vapply(l, function(x) length(x$coef), integer(1L))
res$k <- res$df + 1L
res$`-2logLik` <- vapply(l, extractIC, p = 0, numeric(1L))
res$`-2logLik diff` <- res$`-2logLik` - min(res$`-2logLik`)
res$AIC <- vapply(l, AIC, numeric(1L))
res$`AIC diff` <- res$AIC - min(res$AIC)
res$BIC <- vapply(l, BIC, numeric(1L))
res$`BIC diff` <- res$BIC - min(res$BIC)
if (!missing(p)) {
res[[sprintf('p=%s', p)]] <- pv <- vapply(l, extractIC, p = p, numeric(1L))
res[[sprintf('p=%s diff', p)]] <- pv - min(pv)
}
res <- as.data.frame(res, check.names = FALSE)
rownames(res) <- seq.int(n) - 1L
title <- 'Model selection table\n'
# model <- sprintf('Model %s: %s', rownames(res), model)
structure(
res, heading = c(title, unlist(model), ''),
class = c('anova', 'data.frame')
)
}
#' Wald test for \code{crr} model coefficients
#'
#' @param object an object of class \code{\link[cmprsk]{crr}}
#' @param terms a list (optionally named) giving the indices of the
#' coefficients to test; a list of length \code{n} will perform \code{n}
#' tests; default is an overall test and each coefficient individually
#'
#' Alternatively, a (named) vector (or list) of character strings can be
#' used which will \code{\link{grep}} for each pattern in the coefficient
#' names; patterns that match multiple terms will be tested as a group, and
#' only one character string per test will be used
#'
#' @seealso
#' \code{\link[aod]{wald.test}}; \code{\link{crrfit}}; \code{\link{crrFits}}
#'
#' @examples
#' crrs <- crr2(Surv(futime, event(censored) == death) ~
#' age + sex + abo, transplant)
#'
#' ## default is to test all estimates overall and individually
#' crrwald.test(crrs[[1L]])
#' summary(crrs[[1]])$coef[, 'p-value']
#'
#'
#' ## identical ways to call crrwald.test
#' crrwald.test(crrs[[1L]], list(overall = 1:5, Age = 1,
#' Sex = 2, ABO = 3:5))
#'
#' crrwald.test(crrs[[1L]], c(overall = '.', Age = 'age', ## or list()
#' Sex = 'sex', ABO = 'abo'))
#'
#' @export
crrwald.test <- function(object, terms) {
assert_class(object, 'crr')
co <- coef(object)
nn <- names(co)
nt <- seq_along(co)
if (missing(terms)) {
terms <- c(list(nt), nt)
names(terms) <- c('Overall', nn)
}
if (is.character(unlist(terms)))
terms <- lapply(terms, function(x) grep(x, nn))
wald <- lapply(terms, function(x)
wald.test(object[['var']], co, x)$result$chi2)
res <- do.call('rbind', wald)
rownames(res) <- names(terms)
res
}
wald.test <- function (Sigma, b, Terms = NULL, L = NULL, H0 = NULL,
df = NULL, verbose = FALSE) {
## aod::wald.test
if (is.null(Terms) & is.null(L))
stop('One of the arguments Terms or L must be used.')
if (!is.null(Terms) & !is.null(L))
stop('Only one of the arguments Terms or L must be used.')
w <- if (is.null(Terms)) {
Terms <- seq(length(b))[colSums(L) > 0]
nrow(L)
} else length(Terms)
if (is.null(H0))
H0 <- rep_len(0, w)
if (w != length(H0))
stop('Vectors of tested coefficients and of null hypothesis have different lengths\n')
if (is.null(L)) {
L <- matrix(rep_len(0, length(b) * w), ncol = length(b))
for (i in seq.int(w))
L[i, Terms[i]] <- 1
}
dimnames(L) <- list(paste0('L', as.character(seq.int(NROW(L)))), names(b))
f <- L %*% b
V <- Sigma
mat <- qr.solve(L %*% V %*% t(L))
stat <- t(f - H0) %*% mat %*% (f - H0)
p <- pchisq(stat, df = w, lower.tail = FALSE)
res <- if (is.null(df))
list(chi2 = c(chi2 = stat, df = w, P = p))
else {
fstat <- stat / nrow(L)
df1 <- nrow(L)
df2 <- df
list(chi2 = c(chi2 = stat, df = w, P = p),
Ftest = c(Fstat = fstat, df1 = df1, df2 = df2,
P = pf(fstat, df1, df2, lower.tail = FALSE)))
}
structure(
list(Sigma = Sigma, b = b, Terms = Terms, H0 = H0,
L = L, result = res, verbose = verbose, df = df),
class = 'wald.test'
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.