I would like to perform power and sample size calculations for comparison of unpaired receiver-operating characteristic (ROC) curves. I have tried using the power.roc.test function from pROC package on R, but realise that it is meant for paired ROC curves only. Are there any functions on R, or other software that I can use to conduct such analysis? Thank you.
1 Answer
One way is to get the variances of ROC curves by pROC::var
. Followed by calculating the critical value of null hypothesis roc_null
and the power of alternative hypothesis roc_alt
under normal distribution.
critical_value = qnorm(1-0.05/2, mean=auc(roc_null), sd=sqrt(var(roc_null)))
power = 1 - pnorm(critical_value, mean=auc(roc_alt), sd=sqrt(var(roc_alt)))
The simulation below (100 cases and 100 controls with AUC around 0.8) shows that the AUC is almost normally distributed and the estimated variances are pretty accurate. However, the estimated power is not very stable.
# Summary table
#############################
SD_of_auc1 Power
#############################
Empirical 0.03049 0.7959
Estimated 0.03084 0.7479
library(pROC)
set.seed(1)
N_case <- 100 #case
N_control <- 100 #control
y = rep(c(1,0), c(N_case, N_control))
N = 10000
AUC_null = AUC_alt = rep(NA, N)
sd_null = sd_alt = rep(NA, N)
power_pair = power = rep(NA, N)
for (i in 1:N) {
x_null = rnorm(N_case+N_control, mean=y*1.2) #null hypothesis
x_alt = rnorm(N_case+N_control, mean=y*1.65) #alternative hypothesis
roc_null = roc(y, x_null, direction='<', quiet=TRUE)
roc_alt = roc(y, x_alt, direction='<', quiet=TRUE)
AUC_null[i] = auc(roc_null)
AUC_alt[i] = auc(roc_alt)
sd_null[i] = sqrt(var(roc_null))
sd_alt[i] = sqrt(var(roc_alt))
crit = qnorm(1-0.05/2, mean=auc(roc_null), sd=sqrt(var(roc_null)))
power[i] = 1 - pnorm(crit, mean=auc(roc_alt), sd=sqrt(var(roc_alt)))
}
# Empirical results
power_true = sum(AUC_alt > quantile(AUC_null, 1-0.05/2)) / length(AUC_alt) #Power 0.7959
sd_null_true = sd(AUC_null) #0.03049554
sd_alt_true = sd(AUC_alt) #0.02404829
# Estimated results by pROC
quantile(power, c(0.25,0.5,0.75))
# 25% 50% 75%
#0.3631244 0.7479268 0.9654586
quantile(sd_null, c(0.25,0.5,0.75))
# 25% 50% 75%
#0.02921476 0.03084504 0.03232556
#normal distribution
hist(AUC_null, breaks=100, probability=TRUE)
x_ = seq(0.5, 1, by=0.001)
lines(x_, dnorm(x_, mean(AUC_null), sd(AUC_null)), col='red')