Preprocess clustered binary data before inference on sample proportion

preprocess.clustered.binary.data <- function(x, n)
{
    ## Purpose: Preprocess clustered binary data before inference on sample proportion
    ##          Reference: Rao, J. N. K., & Scott, A. J. (1992). A simple method for
    ##                     the analysis of clustered binary data. Biometrics, 577-585.
    ## Keywords: variance adjustment, ratio estimator, correlated data
    ## Arguments:
    ##   x: a vector for the numerators across clusters.
    ##   n: a vector for the denominators across clusters.
    ## Return: a pre-processed numerator and denominator for overall sample proportion
    ## Author: Feiming Chen
    ## ________________________________________________

    m <- length(x)
    cat("Number of Clusters =", m, "\n")

    n0 <- sum(n)
    cat("Raw Sample Size =", n0, "\n")

    x0 <- sum(x)
    cat("Raw Incidence Count =", x0, "\n")

    cat("\nEstimates Regarding the Overall Sample Proportion:\n")

    p <- x0 / n0
    cat("  Point Estimate =", round(p, 4), "\n")

    v0 <- p * (1 - p) / n0
    cat("  Naive Binomial Variance =", round(v0, 4), "\n")

    r <- x - n * p
    v <- m * sum(r^2) / (m - 1) / n0^2
    cat("  Correct Variance  =", round(v, 4), "\n")

    d <- v / v0
    cat("\nDesign Effect (Variance Inflation Factor due to Clustering) =", round(d, 4), "\n")

    n1 <- n0 / d
    cat("Effective Sample Size (n) =", round(n1), "\n")

    x1 <- x0 / d
    cat("Effective Incidence Count (x) =", round(x1), "\n")

    ## Return the pre-processed numerator (x) and denominator (n) for
    ## the overall sample proportion (p = x / n). 
    list(x = x1, n = n1)

}
if (F) {                                # Unit Test
    x = c(1,1,2,0,5,0,1,4,0,1,0,0,0,0,0,4,3,0,1,1,1,2,0,1,0,2,1,0,1,5,0,0,0,0,0,0,2,0,2,0,0,2,1,2,2,0,0,1,0,1) 
    n = c(2,2,3,0,7,1,2,5,1,2,0,0,4,2,0,7,4,1,1,1,4,2,3,1,0,2,1,0,1,5,3,0,1,3,0,0,2,1,2,0,0,5,3,2,2,1,0,2,1,1) 

    a <- preprocess.clustered.binary.data(x, n)
    ## Number of Clusters = 50 
    ## Raw Sample Size = 93 
    ## Raw Incidence Count = 50 

    ## Estimates Regarding the Overall Sample Proportion:
    ##   Point Estimate = 0.5376 
    ##   Naive Binomial Variance = 0.0027 
    ##   Correct Variance  = 0.004 

    ## Design Effect (Variance Inflation Factor due to Clustering) = 1.4895 
    ## Effective Sample Size (n) = 62 
    ## Effective Incidence Count (x) = 34 

    prop.test(a$x, a$n)
    ## 95 percent confidence interval:
    ##  0.40774 0.66290
    ## sample estimates:
    ##       p 
    ## 0.53763 

    ## Compare to Bootstrap Confidence Intervals
    library(boot)
    r <- boot(data.frame(x=x, n=n), function(dat, idx) { d <- dat[idx,]; sum(d$x)/sum(d$n)}, R = 10000)
    boot.ci(r)
    ## Bootstrap Percentile CI: ( 0.4118,  0.6569 ), which is a bit tighter than the VIF method. 

    ## Compare to Unadjusted (Wrong) CI: 
    prop.test(sum(x), sum(n))
    ## 95 percent confidence interval:
    ##  0.43159 0.64053, which is too narrow and is wrong. 

}

Bootstrap Confidence Interval for Data with Repeated Measures

stat.bootstrap.cluster <- function(id, val, func = mean, boot.size = 10000)
{
    ## Purpose: Calculate bootstrap-based 95% confidence interval for data with repeated measures
    ## Arguments:
    ##   id: uniquely identifies a subject (patient)
    ##   val: a numeric vector to be summarized.  It contains repeated measures per subject.
    ##   func: a function for calculating the summary statistic from a numeric vector.  Default to "mean". 
    ##   boot.size: number of bootstrap samples.  Default to 10000. 
    ## Return: Point estimate, bootstrap percentile 95% CI, histogram for the bootstrap distribution of the target statistic
    ## Author: Feiming Chen
    ## ________________________________________________

    fname <- deparse(substitute(func))
    ans <- func(val)
    cat("Summary Statistic:", fname, "=", ans, "\n")

    unique.ID <- unique(id)

    set.seed(1)
    replicate(n = boot.size, {
        s <- sample(unique.ID, replace = TRUE)  # a bootstrap sample of patient ID's
        ## find all rows with the ID's in the bootstrap sample of ID's 
        v <- c()
        for (j in s) v <- c(v, val[id == j]) # a new bootstrap sample
        func(v)             #  and its statistic
    }) -> est.boot

    hist(est.boot, xlab = fname, main = paste("Bootstrap Distribution: ", fname))
    cat("95% Bootstrap Percentile Confidence Interval:\n")
    quantile(est.boot, c(0.025, 0.975))
}
if (F) {                                # Unit Test
    id <- c(1, 1, 1, 2, 2, 3)
    val <- c(3, 3, 3, 4, 4, 5)
    stat.bootstrap.cluster(id, val)          # compare with t-test based CI: (2.8098 4.5235)
    ## Summary Statistic: mean = 3.6667 
    ## 95% Bootstrap Percentile Confidence Interval:
    ##  2.5% 97.5% 
    ##     3     5 
    stat.bootstrap.cluster(id, val, func = sd)
    ## Summary Statistic: sd = 0.8165 
    ## 95% Bootstrap Percentile Confidence Interval:
    ##   2.5%  97.5% 
    ## 0.0000 1.095
}

Calculate the 95% CI of the ratio of two proportions, with application to diagnostic likelihood ratios

ratio.of.two.prop <- function(x, n)
{
    ## Purpose: Calculate the 95% CI of the ratio of two proportions, with application to diagnostic likelihood ratios. 
    ## Arguments:
    ##   x: a vector of two integers (numerator of the two proportions). 
    ##   n: a vector of two integers (denominator of the two proportions).
    ## Return: The point estimate and the 95% confidence interval of the ratio of two proportions, where the ratio is the first proportion divided by the second proportion. 
    ## Author: Feiming Chen
    ## ________________________________________________

    if (x[1] == 0 || x[1] == n[1] || x[2] == 0 || x[2] == n[2]) {
        x <- x + 0.5
        n <- n + 1
    }

    p <- x / n

    se <- sqrt((1 - p[1]) / (p[1] * n[1]) + (1 - p[2]) / (p[2] * n[2]))
    
    pp <- p[1] / p[2]                   # point estimate
    LB <- exp(log(pp) - 1.96 * se)
    UB <- exp(log(pp) + 1.96 * se)

    c(pp, LB, UB)
}
if (F) {                                # Unit Test
    ratio.of.two.prop(c(431, 30), c(460, 146))
    ## [1] 4.5599 3.3116 6.2786
    ratio.of.two.prop(c(0, 146), c(460, 146))
    ## [1] 0.00108830 0.00006817 0.01737423   
}

Calculate the correlation within subjects

repeated.measures.correlation <- function(subject, x, y)
{
    ## Purpose: Calculate the correlation within subjects based on the paper:
    ##    Bland, J. Martin, and Douglas G. Altman. "Statistics notes: Calculating
    ##    correlation coefficients with repeated observations: Part 1—correlation
    ##    within subjects." Bmj 310.6977 (1995): 446.
    ## Arguments:
    ##    subject: ID for subject.  Each subject has multiple correlated measurements (x, y)
    ##    x: the variable that correlates with y.
    ##    y: the variable that correlates with x. 
    ## Return: 
    ## Author: Feiming Chen
    ## ________________________________________________

    dat <- data.frame(subject = factor(subject), x = x, y = y)
    b <- summary(aov(y ~ subject + x, dat))
    print(b)
    s <- b[[1]][[2]][-1]    
    r0 <- sqrt(s[1]/sum(s))             # within-subject correlation
    cat("\nWithin-Subject Correlation =", round(r0, 4), "\n")
    r0
}
if (F) {                                # Unit Test
    ## simulated data
    r <- 0.7                            # True Within-Subject Correlation
    z <- atanh(r)
    N <- 30
    se <- 1/sqrt(N-3)
    zz <- rnorm(N, z, se)
    rr <- tanh(zz)

    library(MASS)
    set.seed(1)
    dat <- data.frame(subject = factor(rep(1:N, each = 100)))
    dat$x <- dat$y <- NA
    for (i in 1:N) {
        a <- mvrnorm(n = 100, mu = c(0, 0), Sigma = matrix(c(1, rr[i], rr[i], 1), 2))
        dat[dat$subject == i, c(2, 3)] <- a
    }

    ## ANOVA analysis
    b <- summary(aov(y ~ subject + x, dat))
    s <- b[[1]][[2]][-1]    
    r0 <- sqrt(s[1]/sum(s))             # 0.70471

    ## using average of correlations
    d1 <- as.numeric(by(dat, dat$subject, function(d) cor(d$x, d$y)))
    d2 <- atanh(d1)
    d3 <- sample.mean.CI(d2)
    tanh(d3)                            # mean = 0.72 (0.68,0.75) (with Fisher's Z)
    sample.mean.CI(d1)                  # mean = 0.71 (0.67,0.74) (without Fisher's Z)
    repeated.measures.correlation(dat$subject, dat$x, dat$y) # 0.7047
}

Fisher's Z Transform

Fisher.Z <- function(x, inverse = FALSE)
{
    ## Purpose: Fisher's Z Transform (application to correlation confidence interval)
    ## Arguments:
    ##   x: a positive numeric vector
    ##   inverse: if TRUE, do the inverse transform. 
    ## Return: a transformed (or back-transformed) vector
    ## Author: Feiming Chen
    ## ________________________________________________

    if (inverse) {
        y <- tanh(x)
        ## equivalent to the following:
        ## e <- exp(2 * x)
        ## y <- (e - 1) / (e + 1)
    } else {
        y <- atanh(x)
        ## equivalent to the following:
        ## y <- 1/2 * log((1 + x) / (1 - x))
    }
    y
}
if (F) {                                # Unit Test
    Fisher.Z(0.8)                       # 1.0986
    Fisher.Z(1.2, inverse = TRUE)       # 0.83365

    x <- rbeta(10, 1, 1)                # simulate correlations
    Fisher.Z(x)
    Fisher.Z(x, inverse = TRUE)
    all.equal(x, Fisher.Z(Fisher.Z(x), inverse = TRUE)) # TRUE
}

Calculate bootstrap confidence intervals for a statistic of a numeric vector

Using the library "boot"

my.boot <- function(x, func = mean, R = 10000) {
    ## Purpose: Calculate bootstrap confidence interval for a statistic of a numeric vector
    ## Arguments:
    ##   x: a vector of numerical values
    ##   func: a function to calculate the statistic (default to "mean")
    ##   R: number of replicate bootstrap samples (default to 10000)
    ## Return: the bootstrap confidence intervals (percentile and BCa)
    ## Author: Feiming Chen
    ## ________________________________________________

    library(boot)
    b <- boot(x, function(x, i) func(x[i]), R = R)
    print(boot.ci(b, type = c("perc", "bca")))
}
if (F) {
    x <- rnorm(100)
    my.boot(x)
    my.boot(x, func = median)
    ## Level     Percentile            BCa          
    ## 95%   (-0.1284,  0.2726 )   (-0.1371,  0.2533 )  
}
Not using the library "boot"

stat.bootstrap.independent <- function(val, func = mean, N.bootstrap = 10000)
{
    ## Purpose: Calculate bootstrap-based 95% confidence interval for independent univariate data
    ## Arguments:
    ##   val: a numeric vector to be summarized. 
    ##   func: a function for calculating the summary statistic from a numeric vector.  Default to "mean". 
    ##   N.bootstrap: number of bootstrap samples.  Default to 10000. 
    ## Return: Point estimate, bootstrap percentile 95% CI, histogram for the bootstrap distribution of the target statistic
    ## Author: Feiming Chen
    ## ________________________________________________

    fname <- deparse(substitute(func))
    ans <- func(val)
    cat("Summary Statistic:", fname, "=", ans, "\n")

    N <- length(val)
    r <- rep(NA, N.bootstrap)
    set.seed(1)
    for (i in 1:N.bootstrap) {
        s <- sample(val, size = N, replace = TRUE)
        r[i] <- func(s)
    }

    hist(r, xlab = fname, main = paste("Bootstrap Distribution: ", fname))

    cat("95% Bootstrap Percentile Confidence Interval:\n")
    quantile(r, c(0.025, 0.975))
}
if (F) {                                # Unit Test
    val <- c(3, 3, 3, 4, 4, 5)
    stat.bootstrap.independent(val)          # compare with t-test based CI: (2.8098 4.5235)
    ## Summary Statistic: mean = 3.6667 
    ## 95% Bootstrap Percentile Confidence Interval:
    ##   2.5%  97.5% 
    ## 3.1667 4.3333 
}

Simple Survival Analysis of time-to-event data with missing data

my.survival.analysis <- function(dat, study.end.time = NULL)
{
    ## Purpose: Survival Analysis of time-to-event data with missing data. Assume two groups. 
    ## Arguments:
    ##   dat: time-to-event data that look like (variable name is fixed):
    ## 
    ##         time status  arm
    ##           9      1    1
    ##          13      1    0
    ##          13      0    1
    ## 
    ##     The meaning of each variable is as follows: 
    ##       time: the follow-up time (numeric) for right censored data.
    ##       status: the status indicator. 0 = right censored, 1 = event. 
    ##       arm: the group indicator for comparison. 0 = control group, 1 = treatment group. 
    ## 
    ##   study.end.time: truncation time (the end of the time window in which to compare two groups).
    ##                   It needs to be smaller than the default value (minimum of the largest observed time in each of the two groups)
    ## 
    ## Return: Survival Data Analysis Summary and Plots
    ## Author: Feiming Chen
    ## ________________________________________________

    library(survival)

    dat <<- dat                         # so that "cox.zph(f2)" below works. 

    fmla <- Surv(time, status) ~ arm    # formula for survival time analysis

    f <- survfit(formula = fmla, data = dat)

    print(f)

    ## Simple plot of survival curves
    plot(f, lwd = 2, col = c("red", "blue"), mark.time = TRUE, xlab = "Time", ylab = "Survival Probability")
    legend("topright", legend = names(f$strata), bg = "lightyellow", col = c("red", "blue"), title = "Groups", lwd = 2)

    cat("\n------------------------------\nLog-Rank test for comparing two survival curves:\n")
    print(survdiff(fmla, dat))

    cat("\n------------------------------\nSurvival Probability at Study End:\n")
    if (missing(study.end.time)) {
        study.end.time <- max(f$time)
        tau <- NULL
    } else  tau <- study.end.time

    ans <- summary(f, times = study.end.time, extend = TRUE, rmean = study.end.time)
    print(ans)
    ## Post-test disease risk = 1 - survival.probability.at.end

    ## ARR (Absolute Risk Reduction) - Check the direction of the risk reduction manually! 
    ARR <- diff(ans$surv)
    cat("ARR (Absolute Risk Reduction) =", round(ARR, 2), "\n")
    NNT <- ceiling(1 / ARR)
    cat("NNT (Number Needed to Treat) =", NNT, "\n")
    OR <- ans$surv[1] / (1 - ans$surv[1]) / (ans$surv[2] / (1 - ans$surv[2]))
    cat("Odds Ratio:", round(OR, 2), "\n")

    cat("\n------------------------------\nUsing Cox Regression Model:\n")
    f2 <- coxph(formula = fmla, data = dat)
    print(f2)
    cat("Hazard Ratio =", round(exp(f2$coefficients), 2), "\n")

    cat("\n------------------------------\nAssessing Proportional Hazards Assumptions:\n")
    print(cox.zph(f2))

    cat("\n------------------------------\nRMST (Restricted Mean Survival Time) Analysis\n")
    library(survRM2)
    r <- rmst2(time = dat$time, status = dat$status, arm = dat$arm, tau = tau)
    print(r)
    ## plot(r, xlab = "Time", ylab = "Survival Probability")

    f
}

if (F) {                                # Unit Test
    d <- data.frame(time = aml$time, status = aml$status, arm = ifelse(aml$x == "Maintained", 1, 0))
    f <- my.survival.analysis(d, study.end.time = 35)

    d <- survRM2::rmst2.sample.data()
    f <- my.survival.analysis(d, study.end.time = 10)

    ## cannot do it inside a function.
    ## make survival curve with cumulative events/incidences
    f <- survfit(Surv(time, status) ~ x, data = aml)
    library(survminer)
    ggsurvplot(f, data = aml, fun = "event", conf.int = TRUE, pval = TRUE, risk.table = "abs_pct", ggtheme = theme_bw(), palette = c("red", "blue"), risk.table.y.text.col = TRUE, risk.table.y.text = FALSE, ncensor.plot = TRUE)

}

Call: survfit(formula = fmla, data = dat)

       n events median 0.95LCL 0.95UCL
arm=0 12     11     23       8      NA
arm=1 11      7     31      18      NA

------------------------------
Log-Rank test for comparing two survival curves:
Call:
survdiff(formula = fmla, data = dat)

       N Observed Expected (O-E)^2/E (O-E)^2/V
arm=0 12       11     7.31      1.86       3.4
arm=1 11        7    10.69      1.27       3.4

 Chisq= 3.4  on 1 degrees of freedom, p= 0.07 

------------------------------
Survival Probability at Study End:
Call: survfit(formula = fmla, data = dat)

                arm=0 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     35.0000       2.0000       9.0000       0.1944       0.1219       0.0569       0.6642 

                arm=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      35.000        3.000        6.000        0.368        0.163        0.155        0.875 

ARR (Absolute Risk Reduction) = 0.17 
NNT (Number Needed to Treat) = 6 
Odds Ratio: 0.41 

------------------------------
Using Cox Regression Model:
Call:
coxph(formula = fmla, data = dat)

     coef exp(coef) se(coef)    z    p
arm -0.92      0.40     0.51 -1.8 0.07

Likelihood ratio test=3.4  on 1 df, p=0.066
n= 23, number of events= 18 
Hazard Ratio = 0.4 

------------------------------
Assessing Proportional Hazards Assumptions:
         chisq df    p
arm    0.00788  1 0.93
GLOBAL 0.00788  1 0.93

------------------------------
RMST (Restricted Mean Survival Time) Analysis

The truncation time: tau = 35  was specified. 

Restricted Mean Survival Time (RMST) by arm 
               Est.    se lower .95 upper .95
RMST (arm=1) 27.057 2.906    21.361    32.753
RMST (arm=0) 20.958 3.456    14.185    27.732


Restricted Mean Time Lost (RMTL) by arm 
               Est.    se lower .95 upper .95
RMTL (arm=1)  7.943 2.906     2.247    13.639
RMTL (arm=0) 14.042 3.456     7.268    20.815


Between-group contrast 
                      Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) 6.098    -2.752    14.949 0.177
RMST (arm=1)/(arm=0) 1.291     0.878     1.899 0.194
RMTL (arm=1)/(arm=0) 0.566     0.238     1.342 0.196