Prevalence Scaling

prevalence.scaling <- function(model.P, model.prevalence = 0.5, target.prevalence = 0.1)
{
    ## Purpose: Convert a probability based on training data prevalence to
    ##          a probability based on target population prevalence. 
    ## Arguments:
    ##   model.P: a vector of model-based probability (risk score)
    ##   model.prevalence: prevalence of the model's training data
    ##   target.prevalence: prevalence of the target population
    ## Return:
    ##   a vector of scaled probabilities reflecting the prevalence of the target population
    ## Author: Feiming Chen
    ## ________________________________________________

    ## convert a probability to the odds
    O <- function(p)  p / (1-p)

    1 / ( 1 + O(model.prevalence) / O(target.prevalence) / O(model.P) )
}
if (F) {                                # Unit Test
    set.seed(1)
    prevalence.scaling(runif(5))
    ## [1] 0.038614 0.061784 0.129688 0.523663 0.027304
}

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