S3 method version of epiR::epi.tests()

References

Improvements

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 the original package and the S3 version

## 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

Test run along with the original

### 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

Code

## 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)
}