Power analysis for Receiver Operating Characteristic (ROC) curves in R

Question

I need to perform a priori power analysis to determine the minimum sample size for the study we have designed.

In this study, we will utilize Receiver Operating Characteristic (ROC) curves to assess the effectiveness of three distinct biomarkers (A, B, C) in differentiating the stages of a given disease: normal (without the disease), mild stage, and severe stage.

In addition, these curves will be analyzed separately for participants who carry a specific protein gene (carriers) and those without this risk factor (non-carriers). It is important to note that the biomarkers (A, B, C) will be measured as continuous variables.

I would greatly appreciate it if you could analyze the function implemented in R for calculating power analysis.

My main concern is that, in the implemented function, using 50 or 5 participants for each group yields high power results (> 80). I would need to select very low effect sizes (let's say, 0.0, 0.1, and 0.2) to achieve low power in my analysis. In my field of research, the values of 0.00 (baseline), 0.2 (mild), and 0.40 (severe) are deemed a conservative approach,

Thank you for your assistance!

library(pROC)

set.seed(123)

# Define parameters
n_sim <- 1000            # Number of simulations
n <- 50                  # Small sample size per group for testing
mu <- c(0.00, 0.2, 0.40) # Effect sizes for Normal, Mild, Severe / We set an effect size of zero for Normal because it works as the reference group (baseline)
sigma <- 0.1             # Standard deviation (assuming equal variances)
alpha <- 0.05            # Significance level

# Storage for results
power_results <- data.frame(protein_status = rep(c("carrier", "non-carrier"), each = 3),
                            comparison = rep(c("Normal vs Mild", "Normal vs Severe", "Mild vs Severe"), 2),
                            biomarker = rep(c("A", "B", "C"), each = 6),
                            power = NA)

# Loop over protein status, biomarkers, and pairwise comparisons
for (protein in c("carrier", "non-carrier")) {
  for (biomarker in c("A", "B", "C")) {
    for (comparison in c("Normal vs Mild", "Normal vs Severe", "Mild vs Severe")) {
      
      sig_count <- 0
      
      # Run simulations
      for (i in 1:n_sim) {
        
        # Simulate data for each group
        data_normal <- rnorm(n, mean = mu[1], sd = sigma)
        data_mild <- rnorm(n, mean = mu[2], sd = sigma)
        data_severe <- rnorm(n, mean = mu[3], sd = sigma)
        
        # Select data based on comparison
        if (comparison == "Normal vs Mild") {
          data <- data.frame(
            value = c(data_normal, data_mild),
            group = factor(rep(c("normal", "mild"), each = n))
          )
        } else if (comparison == "Normal vs Severe") {
          data <- data.frame(
            value = c(data_normal, data_severe),
            group = factor(rep(c("normal", "severe"), each = n))
          )
        } else if (comparison == "Mild vs Severe") {
          data <- data.frame(
            value = c(data_mild, data_severe),
            group = factor(rep(c("mild", "severe"), each = n))
          )
        }
        
        # Calculate ROC curve and suppress messages
        roc_res <- suppressMessages(roc(data$group, data$value, levels = rev(levels(data$group))))
        
        # Calculate the confidence interval of the AUC
        ci <- ci.auc(roc_res)
        
        # Debugging: Print AUC and confidence interval for the first few simulations
        if (i <= 10) {
          cat("Sim:", i, "AUC:", auc(roc_res), "CI:", ci, "\n")
        }
        
        # Check if the lower bound of the confidence interval is greater than 0.5
        if (ci[1] > 0.5) {
          sig_count <- sig_count + 1
        }
      }
      
      # Calculate power
      power <- sig_count / n_sim
      
      # Store results
      power_results[power_results$protein_status == protein & power_results$biomarker == biomarker & power_results$comparison == comparison, "power"] <- power
    }
  }
}

# Print results
print(power_results)

Answer 1

Concordance probabilities, AKA $c$-index or AUROC are decent measures of pure predictive discrimination but are not recommended for comparison of two models, because of low power. Taking the difference between two rank measures is not a good idea in general, e.g., we never take the difference in two Wilcoxon statistics, which is what you are doing by comparing two AUROCs. Instead use full-information measures such as likelihood ratio $\chi^2$ statistics (see here) or methods discussed here.