Original epiR package: http://cran.r-project.org/package=epiR
S3 method version of epiR::epi.2by2(): http://rpubs.com/kaz_yos/s3-epi-2by2
The new function returns the result as an object instead of printing the results. The S3 print method will take care of the printing. The summary method returns the data frame with data. The output object actually holds all the values created along the execution.
## Load epiR for comparison
library(epiR)
Loading required package: survival
Loading required package: splines
Package epiR 0.9-51 is loaded
Type help(epi.about) for summary information
### Define new function, class, and methods
## New function using S3 methods
source("./epi.tests.R")
## Print method
source("./epi.tests.print.R")
## Summary method to print the "verbose" part
source("./epi.tests.summary.R")
## Show available methods
methods(class = "epi.tests")
[1] print.epi.tests summary.epi.tests
### Implementation of S3 methods for epiR::epi.tests()
## Disease + Disease - Total
## Test + a b a + b
## Test - c d c + d
## Total a + c b + d a + b + c + d
## Create test data
dat <- as.table(matrix(c(670,202,74,640), nrow = 2, byrow = TRUE))
colnames(dat) <- c("Dis+","Dis-")
rownames(dat) <- c("Test+","Test-")
dat
Dis+ Dis-
Test+ 670 202
Test- 74 640
## Run the original
epi.tests(dat, conf.level = 0.95, verbose = FALSE)
Disease + Disease - Total
Test + 670 202 872
Test - 74 640 714
Total 744 842 1586
Point estimates and 95 % CIs:
---------------------------------------------------------
Apparent prevalence 0.55 (0.52, 0.57)
True prevalence 0.47 (0.44, 0.49)
Sensitivity 0.9 (0.88, 0.92)
Specificity 0.76 (0.73, 0.79)
Positive predictive value 0.77 (0.74, 0.8)
Negative predictive value 0.9 (0.87, 0.92)
---------------------------------------------------------
## Run the new one
out <- epi.tests.new(dat, conf.level = 0.95)
## Implicit print method. (The sprintf() formats the results better.)
out
Disease + Disease - Total
Test + 670 202 872
Test - 74 640 714
Total 744 842 1586
Point estimates and 95 % CIs:
---------------------------------------------------------
Apparent prevalence 0.55 (0.52, 0.57)
True prevalence 0.47 (0.44, 0.49)
Sensitivity 0.90 (0.88, 0.92)
Specificity 0.76 (0.73, 0.79)
Positive predictive value 0.77 (0.74, 0.80)
Negative predictive value 0.90 (0.87, 0.92)
---------------------------------------------------------
## Summary method returns a data frame
summary(out)
est lower upper
aprev 0.5498 0.5249 0.5745
tprev 0.4691 0.4443 0.4940
se 0.9005 0.8767 0.9211
sp 0.7601 0.7298 0.7886
diag.acc 0.8260 0.8064 0.8443
diag.or 28.6861 21.5182 38.2417
nnd 1.5137 1.4091 1.6487
youden 0.6606 0.6065 0.7097
ppv 0.7683 0.7389 0.7960
npv 0.8964 0.8716 0.9177
plr 3.7537 3.3207 4.2432
nlr 0.1309 0.1051 0.1630
## The main function
epi.tests.new
function (dat, conf.level = 0.95)
{
elements <- list()
elements <- within(elements, {
N. <- 1 - ((1 - conf.level)/2)
z <- qnorm(N., mean = 0, sd = 1)
.funincrisk <- function(cdat, conf.level) {
N. <- 1 - ((1 - conf.level)/2)
a <- cdat[, 1]
n <- cdat[, 2]
b <- n - a
p <- a/n
a. <- ifelse(a == 0, a + 1, a)
b. <- ifelse(b == 0, b + 1, b)
low <- a./(a. + (b. + 1) * (1/qf(1 - N., 2 * a.,
2 * b. + 2)))
up <- (a. + 1)/(a. + 1 + b./(1/qf(1 - N., 2 * b.,
2 * a. + 2)))
low <- ifelse(a == 0, 0, low)
up <- ifelse(a == n, 1, up)
rval <- data.frame(est = p, lower = low, upper = up)
rval
}
a <- dat[1]
b <- dat[3]
c <- dat[2]
d <- dat[4]
M1 <- a + c
M0 <- b + d
N1 <- a + b
N0 <- c + d
total <- a + b + c + d
tdat <- as.matrix(cbind(M1, total))
trval <- .funincrisk(tdat, conf.level)
tp <- trval$est
tp.low <- trval$lower
tp.up <- trval$upper
tprev <- data.frame(est = tp, lower = tp.low, upper = tp.up)
tdat <- as.matrix(cbind(N1, total))
trval <- .funincrisk(tdat, conf.level)
ap <- trval$est
ap.low <- trval$lower
ap.up <- trval$upper
aprev <- data.frame(est = ap, lower = ap.low, upper = ap.up)
tdat <- as.matrix(cbind(a, M1))
trval <- .funincrisk(tdat, conf.level)
se <- trval$est
se.low <- trval$lower
se.up <- trval$upper
sensitivity <- data.frame(est = se, lower = se.low, upper = se.up)
tdat <- as.matrix(cbind(d, M0))
trval <- .funincrisk(tdat, conf.level)
sp <- trval$est
sp.low <- trval$lower
sp.up <- trval$upper
specificity <- data.frame(est = sp, lower = sp.low, upper = sp.up)
tdat <- as.matrix(cbind(a, N1))
trval <- .funincrisk(tdat, conf.level)
ppv <- trval$est
ppv.low <- trval$lower
ppv.up <- trval$upper
pv.positive <- data.frame(est = ppv, lower = ppv.low,
upper = ppv.up)
tdat <- as.matrix(cbind(d, N0))
trval <- .funincrisk(tdat, conf.level)
npv <- trval$est
npv.low <- trval$lower
npv.up <- trval$upper
pv.negative <- data.frame(est = npv, lower = npv.low,
upper = npv.up)
lrpos <- (a/M1)/(1 - (d/M0))
lrpos.low <- exp(log(lrpos) - z * sqrt((1 - se)/(M1 *
se) + (sp)/(M0 * (1 - sp))))
lrpos.up <- exp(log(lrpos) + z * sqrt((1 - se)/(M1 *
se) + (sp)/(M0 * (1 - sp))))
lr.positive <- data.frame(est = lrpos, lower = lrpos.low,
upper = lrpos.up)
lrneg <- (1 - (a/M1))/(d/M0)
lrneg.low <- exp(log(lrneg) - z * sqrt((se)/(M1 * (1 -
se)) + (1 - sp)/(M0 * (sp))))
lrneg.up <- exp(log(lrneg) + z * sqrt((se)/(M1 * (1 -
se)) + (1 - sp)/(M0 * (sp))))
lr.negative <- data.frame(est = lrneg, lower = lrneg.low,
upper = lrneg.up)
tdat <- as.matrix(cbind((a + d), total))
trval <- .funincrisk(tdat, conf.level)
da <- trval$est
da.low <- trval$lower
da.up <- trval$upper
diag.acc <- data.frame(est = da, lower = da.low, upper = da.up)
dOR.p <- (a * d)/(b * c)
lndOR <- log(dOR.p)
lndOR.var <- 1/a + 1/b + 1/c + 1/d
lndOR.se <- sqrt(1/a + 1/b + 1/c + 1/d)
lndOR.l <- lndOR - (z * lndOR.se)
lndOR.u <- lndOR + (z * lndOR.se)
dOR.se <- exp(lndOR.se)
dOR.low <- exp(lndOR.l)
dOR.up <- exp(lndOR.u)
diag.or <- data.frame(est = dOR.p, lower = dOR.low, upper = dOR.up)
ndx <- 1/(se - (1 - sp))
ndx.1 <- 1/(se.low - (1 - sp.low))
ndx.2 <- 1/(se.up - (1 - sp.up))
ndx.low <- min(ndx.1, ndx.2)
ndx.up <- max(ndx.1, ndx.2)
nnd <- data.frame(est = ndx, lower = ndx.low, upper = ndx.up)
c.p <- se - (1 - sp)
c.1 <- se.low - (1 - sp.low)
c.2 <- se.up - (1 - sp.up)
c.low <- min(c.1, c.2)
c.up <- max(c.1, c.2)
youden <- data.frame(est = c.p, lower = c.low, upper = c.up)
})
rval <- list(aprev = elements$aprev, tprev = elements$tprev,
se = elements$sensitivity, sp = elements$specificity,
diag.acc = elements$diag.acc, diag.or = elements$diag.or,
nnd = elements$nnd, youden = elements$youden, ppv = elements$pv.positive,
npv = elements$pv.negative, plr = elements$lr.positive,
nlr = elements$lr.negative)
r1 <- with(elements, c(a, b, N1))
r2 <- with(elements, c(c, d, N0))
r3 <- with(elements, c(M1, M0, M0 + M1))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total")
rownames(tab) <- c("Test +", "Test -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(conf.level = conf.level, elements = elements,
rval = rval, tab = tab)
class(out) <- "epi.tests"
return(out)
}
## Print method
print.epi.tests
function (obj)
{
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%", "CIs:")
cat("\n---------------------------------------------------------")
with(obj$rval, {
cat(sprintf("\nApparent prevalence %.2f (%.2f, %.2f)",
aprev$est, aprev$lower, aprev$upper))
cat(sprintf("\nTrue prevalence %.2f (%.2f, %.2f)",
tprev$est, tprev$lower, tprev$upper))
cat(sprintf("\nSensitivity %.2f (%.2f, %.2f)",
se$est, se$lower, se$upper))
cat(sprintf("\nSpecificity %.2f (%.2f, %.2f)",
sp$est, sp$lower, sp$upper))
cat(sprintf("\nPositive predictive value %.2f (%.2f, %.2f)",
ppv$est, ppv$lower, ppv$upper))
cat(sprintf("\nNegative predictive value %.2f (%.2f, %.2f)",
npv$est, npv$lower, npv$upper))
})
cat("\n---------------------------------------------------------")
cat("\n")
}
## Summary method
summary.epi.tests
function (obj)
{
out <- do.call(rbind, obj$rval)
return(out)
}