Confidence Interval of the Difference between Two Correlated Proportions (McNemar Test)

CI.diff.two.correlated.proportions <- function(A, B, C, D, alpha = 0.05)
{
    ## Purpose: CI of the difference between two correlated proportions (McNemar Test)
    ## Arguments:
    ##   A, B, C, D: the 2x2 incidence count table for the matched pairs. 
    ##               A = both positive;
    ##               B = treatment positive, control negative;
    ##               C = treatment negative; control positive;
    ##               D = both negative;
    ##   alpha: type I error. Default to 0.05 (for two-sided CI).
    ## Return: 95% Confidence Interval of the Difference
    ## Author: Feiming Chen
    ## Method: Wald Z Method with Continuity Correction
    ## Reference: Newcombe (1998c), page 2638.
    ## ________________________________________________

    N <- A + B + C + D                  # total number of matched pairs
    Delta <- (B - C) / N
    sw <- sqrt(((A + D) * (B + C) + 4 * B * C) / N^3)
    z <- qnorm(1 - alpha / 2)           # default: 1.96
    corr.term <- 1 / N                  # continuity correction 
    L <- Delta - z * sw - corr.term
    U <- Delta + z * sw + corr.term
    x <- matrix(c(A, B, C, D), 2, 2, byrow = TRUE, 
                dimnames = list(c("Test.POS", "Test.NEG"), c("Control.POS", "Control.NEG")))
    cat("Raw Data (x):\n")
    print(x)
    cat("\nData Converted to Proportions:\n")
    print(round(x/N, 2))
    cat("\nDifference of Two Paired Proportions =", round(Delta, 3), "\n")
    p.discordant <- (B+C) / N
    cat("Proportion of Discordant Pairs =", round(p.discordant, 3), "\n")
    N.future <- sample.size.for.two.correlated.proportions.test(Delta, p.discordant)
    cat("\nCurrent Sample Size =", N, "\n")
    cat("Future  Sample Size =", N.future, "\n")

    print(mcnemar.test(x))
    cat("\nConfidence Interval = (", round(L, 4), ",", round(U, 4), ")\n\n")
}
if (F) {                                # Unit Test
    CI.diff.two.correlated.proportions(12, 9, 2, 21) 
    ## Confidence Interval = ( -0.0037 , 0.3219 )
    CI.diff.two.correlated.proportions(25, 20, 5, 25) 
    ## Raw Data (x):
    ##          Control.POS Control.NEG
    ## Test.POS          25          20
    ## Test.NEG           5          25

    ## Data Converted to Proportions:
    ##          Control.POS Control.NEG
    ## Test.POS        0.33        0.27
    ## Test.NEG        0.07        0.33

    ## Difference of Two Paired Proportions = 0.2 
    ## Proportion of Discordant Pairs = 0.333 

    ## Current Sample Size = 75 
    ## Future  Sample Size = 63 

    ##  McNemar's Chi-squared test with continuity correction

    ## data:  x
    ## McNemar's chi-squared = 7.84, df = 1, p-value = 0.00511


    ## Confidence Interval = ( 0.0641 , 0.3359 )
}

Sample Size for Two Correlated Proportions Based on McNemar Test

sample.size.for.two.correlated.proportions.test <- function(p.diff = 0.1, p.discordant = 0.2, alpha = 0.025, beta = 0.2)
{
    ## Purpose: Sample Size for Two Correlated Proportions Test based on McNemar Test
    ##   H0: p2 <= p1
    ##   H1: p2 >  p1
    ##   (Application: tests of sensitivity, specificity with three-way, matched-pairs study)
    ## Arguments:
    ##   p.diff: Effect size for the difference between two correlated proportions. Default to 10%. 
    ##   p.discordant: Proportion of pairs for which the responses differed. Default to 20%. 
    ##   alpha: type I error. Default to 0.025 (for one-sided test).
    ##   beta: type II error (1 - power). Default to 0.2 (so power is 80%). 
    ## Return:
    ##   Sample Size
    ## Author: Feiming Chen
    ## Reference: Machin, Campbell, Fayers, and Pinol (1997).
    ## ________________________________________________

    p10 <- (p.discordant + p.diff) / 2  # with default, p10 = 0.15
    p01 <- (p.discordant - p.diff) / 2  # default: p01 = 0.05
    OR <- p10 / p01                     # default: OR = 3
    OR.plus.one <- OR + 1               # default: 4
    OR.minus.one <- OR - 1              # default: 2
    
    N = (qnorm(1 - alpha) * OR.plus.one + qnorm(1 - beta) * sqrt(OR.plus.one^2 - OR.minus.one^2 * p.discordant))^2 / (OR.minus.one^2 * p.discordant)
    ceiling(N)
}
if (F) {                                # Unit Test
    sample.size.for.two.correlated.proportions.test() # 155 (PASS output: 155 under normal approx.)
}

Sample Size for One Proportion Test

sample.size.for.one.proportion.test <- function(p0  = 0.8, p1 = 0.9, alpha = 0.025, beta = 0.2)
{
    ## Purpose: Sample Size for One Proportion Test.
    ##   H0: p <= p0
    ##   H1: p >  p0
    ##   (Application: tests of sensitivity, specificity with performance goals)
    ## Arguments:
    ##   p0: performance goal (minimally acceptable proportion). Default to 80%. 
    ##   p1: expected performance (minimal proportion under the alternative hypothesis). Defaul to 90%. 
    ##   alpha: type I error. Default to 0.025 (for one-sided test).
    ##   beta: type II error (1 - power). Default to 0.2 (so power is 80%). 
    ## Return:
    ##   Sample Size
    ## Author: Feiming Chen
    ## Reference: Biswas, Bipasa. "Clinical performance evaluation of molecular diagnostic tests."
    ##            The Journal of Molecular Diagnostics 18.6 (2016): 803-812.
    ## ________________________________________________

    N = (qnorm(1 - alpha) * sqrt(p0 * (1 - p0)) + qnorm(1 - beta) * sqrt(p1 * (1 - p1)))^2 / (p1 - p0)^2
    ceiling(N)
}
if (F) {                                # Unit Test
    sample.size.for.one.proportion.test() # 108 (PASS output: 107 under exact test; 108 under normal approx.)
}

Reading From and Writing To Clipboard (usage: copy a data table from a spreadsheet or paste a table into a spreadsheet)

wc <- wc.linux <- function(x) {
    ## Write to Clipboard
    ## Write a table/data frame "x" to the Clipboard for Excel use.
    ff <- pipe("xclip -i -selection clipboard", "w")
    utils::write.table(x, file=ff, sep="\t", col.names=T, row.names=F, na="")
    close(ff)
}


rc <- function(header, p=TRUE, ...){                       
    ## Read from Clipboard
    ## Check is Header Line exists.
    ## Checking if the first element in the first line is a numeric type or not.
    ## if "p=TRUE", print out the vector definition for copying 
    if (missing(header)) {
        if (is.numeric(unlist(read.delim("clipboard", nrows=1, header=F))[1]))
            header=F
        else
            header=T
    }

    a <- utils::read.delim("clipboard", header=header, as.is = TRUE, ...)
    if (p) {
        for (i in seq(ncol(a))) pvec(a[[i]], var=letters[(22+i) %% 26 + 1])
        return(invisible(a))
    } else {
        if (ncol(a) == 1 || nrow(a) == 1) {    # convert to a vector if there is only one column
            a <- unlist(a)
            cat("\nClipboard is read into a vector of length:", length(a), "\n")
        } else cat("Clipboard is read into a data.frame of dimension:", dim(a), "\n")
        print(head(a, n=3))
        return(invisible(a))
    }
}

ROC Curve Analysis

ROC.curve <- function(R, D, n.thres = 100)
{
    ## Purpose: Perform ROC Curve Analysis
    ## Arguments:
    ##   R: Clinical Reference Standard (0 = Negative, 1 = Positive)
    ##   D: Device Diagnostic Output (a continuous variable)
    ##   n.thres: Number of Thresholds
    ## Return: ROC Curve and a Plot of Sensitivity and Specificity by Thresholds
    ## Author: Feiming Chen
    ## ________________________________________________

    N <- length(R)                      # sample size

    thres <- quantile(D, probs = seq(0, 1, 1 / n.thres)) # list of thresholds
    M <- length(thres)                            # number of thresholds
    sens <- spec <- accu <- rep(0, M)
    for (i in 1:M) {
        D1 <- ifelse(D > thres[i], 1, 0) # convert continuous output to binary output 0-1
        sens[i] <- sum(D1[R == 1]) / sum(R==1)
        spec[i] <- sum(D1[R == 0] == 0) / sum(R==0)
        accu[i] <- (sum(D1[R == 0] == 0) + sum(D1[R == 1])) / N # accuracy
    }
    J <- sens + spec - 1                # Youden's Index

    ## Calculate AUC (Area Under the Curve)
    f <- approxfun(1 - spec, sens, yleft = 0, yright = 1)
    AUC <- integrate(f, lower = 0, upper = 1)$value # c-statistic
    Gini <- 2 * AUC - 1

    plot(1 - spec, sens, type = "l", xlim = c(0, 1), ylim = c(0, 1), lwd = 2, col = "blue",
         xlab = "1 - Specificity (False Positive Rate)", ylab = "Sensitivity (True Positive Rate)")

    title(main = paste0("ROC Curve (AUC = ", round(AUC, 3), ", Gini = ", round(Gini, 3), ")"))

    ## Random Test
    abline(0, 1)                        # uninformative line
    abline(h=c(0, 1), col = "gray")
    abline(v=c(0, 1), col = "gray")

    ## Plot of Sensitivity and Specificity by Thresholds
    dev.new()
    plot(thres, sens, ylim = c(0, 1), main = "Sensitivity/Specificity by Thresholds", type = "l", lwd = 2, col = "blue", xlab = "Thresholds", ylab = "Performance")
    lines(thres, spec, lwd = 2, col = "red")
    lines(thres, J, lwd = 2, col = "orange")
    lines(thres, accu, lwd = 2, lty = 2, col = "black")
    abline(h=c(0,1), col = "gray")
    legend("right", legend= c("Sensitivity", "Specificity", "Youden's Index", "Accuracy"), bg="lightyellow", col= c("blue", "red", "orange", "black"), title="Performance Metrics", lwd=2, lty=c(rep(1, 3), 2))

    res <- data.frame(Threshold = round(thres, 2), Sensitivity = round(sens, 3), Specificity = round(spec, 3), J = round(J, 3), Accuracy = round(accu, 3))
    invisible(res)
}
if (F) {                                # Unit Test
    D <- runif(10000)
    R <- sapply(D, function(p) rbinom(1, size = 1, prob = p)) # perfect probability prediction
    ## R <- sapply(D, function(p) rbinom(1, size = 1, prob = 0.5)) # random probability prediction
    ROC.curve(R, D)
    ## (ROC.curve(R, D, n.thres = 4))
}

## Random Probability Prediction

## Perfect Probability Prediction

Decision Curve Analysis

Net.Benefit <- function(R, D, p.grid)
{
    ## Purpose: Calculate Net Benefit for Decision Curve 
    ## Arguments:
    ##   R: Clinical Reference Standard (0 = Negative, 1 = Positive)
    ##   D: Device Diagnostic Output (0 = Negative, 1 = Positive; OR D = Probability, 0 < D < 1)
    ##   p.grid: The probability levels at which net benefits are to be calculated. 
    ## Return: Net Benefits
    ## Author: Feiming Chen
    ## ________________________________________________

    N <- length(R)                      # sample size
    Net.Benefit <- rep(0, length(p.grid))
    for (i in seq_along(p.grid)) {
        p <- p.grid[i]
        N.TP <- sum( D > p & R == 1 )
        N.FP <- sum( D > p & R == 0 )
        Net.Benefit[i] <- (N.TP - N.FP * p / (1 - p)) / N
    }
    Net.Benefit
}
if (F) {                                # Unit Test
    D <- runif(100)
    R <- sapply(D, function(p) rbinom(1, size = 1, prob = p)) # perfect probability prediction
    p.grid <- seq(0, 0.99, 0.01)           # Grid of indifference probabilities
    Net.Benefit(R, D, p.grid)
}


decision.curve <- function(R, D)
{
    ## Purpose: Perform Decision Curve Analysis
    ## Arguments:
    ##   R: Clinical Reference Standard (0 = Negative, 1 = Positive)
    ##   D: Device Diagnostic Output (0 = Negative, 1 = Positive; OR D = Probability, 0 < D < 1)
    ## Return: Decision Curve
    ## Author: Feiming Chen
    ## ________________________________________________

    N <- length(R)                      # sample size
    p.grid <- seq(0, 0.99, 0.01)        # Grid of indifference probabilities
    NB <- Net.Benefit(R, D, p.grid)
    prevalence <- sum(R == 1) / N
    plot(p.grid, NB, type = "l", xlim = c(0, 1), ylim = c(0, prevalence), lwd = 2, col = "blue", main = "Decision Curve",
         xlab = "Preference (Indifference Probability)", ylab = "Net Benefit",
         sub = paste("Prevalence =", round(100*prevalence, 1), "%"))

    ## Intervention for all
    NB.all <- Net.Benefit(R, rep(1, N), p.grid)
    lines(p.grid, NB.all, type = "l", col = "red", lwd = 1.5)

    ## Perfect Binary Test
    NB.perfect.binary <- Net.Benefit(R, R, p.grid)
    lines(p.grid, NB.perfect.binary, type = "l", col = "orange", lwd = 1.5)
}
if (F) {                                # Unit Test
    D <- runif(100000)
    R <- sapply(D, function(p) rbinom(1, size = 1, prob = p)) # perfect probability prediction
    decision.curve(R, D)
}

 

    
  

if (F) {                                # Simulation Code
    ## Perfect Probability Prediction (D0)
    D <- runif(100000)
    R <- sapply(D, function(p) rbinom(1, size = 1, prob = p)) 
    decision.curve(R, D)

    ## Binary test (B1) with 50% sensitivity and 100% specificity.
    p.grid <- seq(0, 0.99, 0.01)           # Grid of indifference probabilities
    B1 <- sapply(R, function(x) ifelse(x == 1, rbinom(1, 1, 0.5), 0))
    NB.high.spec <- Net.Benefit(R, B1, p.grid)
    lines(p.grid, NB.high.spec, type = "l", col = "orange", lwd = 1.5)

    ## Binary test (B2) with 100% sensitivity and 50% specificity.
    B2 <- sapply(R, function(x) ifelse(x == 0, rbinom(1, 1, 0.5), 1))
    NB.high.sens <- Net.Benefit(R, B2, p.grid)
    lines(p.grid, NB.high.sens, type = "l", col = "orange", lwd = 1.5)

    ## Random prediction model (D1) for:  $D \sim \mbox{Uniform}(0, 1), R \sim \mbox{Bernoulli}(0.5)$
    NB.random <- Net.Benefit(rbinom(100000, 1, 0.5), D, p.grid)
    lines(p.grid, NB.random, type = "l", col = "red", lwd = 1.5)
}

Convert a data set to 0-1 format for binary diagnostic outcomes

binary.recode <- function(dat, pos.label = "Positive")
{
    ## Purpose: Convert a data set to 0-1 format for binary diagnostic outcomes
    ## Arguments:
    ##   dat: a data frame with binary diagnostic outcomes
    ##   pos.label: a label for the positive outcome that is used in the "dat"
    ## Return: a data frame with binary diagnostic outcomes that are coded as numeric 0 or 1. 
    ## Author: Feiming Chen
    ## ________________________________________________

    dat[] <- lapply(dat, function(x, p=pos.label) ifelse(x == p, 1, 0))
    dat
}
if (F) {                                # Unit Test
    binary.recode(data.frame(x = rep("Positive", 2), y = rep("Negative", 2)))
    ##   x y
    ## 1 1 0
    ## 2 1 0
}