swh:1:snp:ff0951ca787d0b7f47dc2335f47fed43820a6324
Tip revision: 8aca325bb4949c0a75f3e27fabfa944a7a8a503d authored by Venkatraman E. Seshan on 18 January 2012, 00:00:00 UTC
version 0.9.9
version 0.9.9
Tip revision: 8aca325
roc.area.test.R
roc.area.test <- function(markers, status) {
if (any(!is.finite(markers))) stop("Marker values should be finite")
if (any(!is.finite(status))) stop("All status should be finite")
markers <- as.matrix(markers)[order(status), , drop=FALSE]
nvar <- ncol(markers)
nn <- sum(status == 0)
nd <- sum(status == 1)
if (nn*nd == 0) stop("Status vector should have least one each of 0 & 1")
n <- nn + nd
zzz <- .Fortran("rocarea",
as.integer(n),
as.integer(nvar),
as.integer(nn),
as.integer(nd),
as.double(markers),
area=as.double(rep(0,nvar)),
jkarea=as.double(matrix(0,n,nvar)),
PACKAGE="clinfun")
out <- NULL
out$area <- zzz$area
areajk <- matrix(zzz$jkarea, n, nvar)
out$var <- ((nn - 1)^2 * var(areajk[1:nn, ]))/nn + ((nd - 1)^2 *
var(areajk[nn + (1:nd), ]))/nd
if (nvar == 2) {
if (out$area[2] != out$area[1]) {
out$stat <- (out$area[2] - out$area[1])/sqrt(out$var[1,1] +
out$var[2,2] - 2*out$var[2,1])
out$p.value <- 2*pnorm(-abs(out$stat))
} else {
out$stat <- 0
out$p.value <- 1
}
}
if (nvar > 2) {
A <- diag(1, nvar) - matrix(1, nvar, nvar)/nvar
x <- (A%*%as.matrix(out$area))[-1]
v <- (A%*%out$var%*%A)[-1,-1,drop=FALSE]
if (qr(v)$rank == nvar - 1) {
out$stat <- sum(solve(v,x)*x)
out$p.value <- 1 - pchisq(out$stat, df=nvar -1)
out$df <- nvar - 1
} else {
warning("Some markers are perfectly correlated")
out$stat <- out$p.value <- out$df <- NA
}
}
class(out) <- "roc.area.test"
out
}
print.roc.area.test <- function(x, ...) {
if (!inherits(x, 'roc.area.test')) stop("Object not of class roc.area.test")
k <- length(x$area)
if (k == 1) {
cat(" AUC =", x$area, "with s.d", sqrt(x$var), "\n")
} else {
if (k == 2) {
msg <- "from standard normal reference"
} else {
msg <- paste("from chi-square (df = ", x$df, ") reference", sep="")
}
cat(" ", k, "markers with AUC", x$area, "\n")
cat(" test statistic =", x$stat, "\n")
cat(" p-value =", x$p.value, msg, "\n")
}
}