Here is a tool used to argue that fault detection algorithm health should be measured using lower confidence interval given sample size instead of prevalence. It is a lower estimate, but it takes into account sample size, and protects from customer-facing falsehoods - which can be expensive.
Note the legend.
Plot:
Code:
# Purpose: Show relationship between prevalence and CI for binomial
# -------------
# housekeeping
# rm(list=ls())
#libraries
require(pacman)
p_load(PropCIs) #for Clopper-Pearson (exact)
#reproducibility
set.seed(09302016)
N_loops <- 100000 #decent sample size
#CI threshold
ci_threshold = 0.95
#stage store for loop
p_true <- numeric(length = N_loops)
n_samp <- numeric(length = N_loops)
prev <- numeric(length = N_loops)
L_CI <- numeric(length = N_loops)
store <- data.frame(p_true, n_samp, prev, L_CI)
### main loop ###
for (i in 1:N_loops){
#how many samples show up (uniform random)
n_samp1 <- sample(x=c(1:300), size = 1,
replace = FALSE)
#true prevalence (uniform random)
p_true1 <- runif(n = 1, min = 0, max = 1)
#draw samples (binomial)
mysamples <- rbinom(n = n_samp1, size=1,
prob = p_true1)
#compute prevalence
prev1 <- sum(mysamples)/length(mysamples)
#compute the lower conf interval
# Agresti-Couli exact as "finish line"
btest <- exactci(x = sum(mysamples),
n = length(mysamples),
conf.level=ci_threshold)
L_CI1 <- btest$conf.int[1]
#pack into store
store$p_true[i] <- p_true1
store$n_samp[i] <- n_samp1
store$prev[i] <- prev1
store$L_CI[i] <- L_CI1
}
#make scatterplots
#take off extreme tails
idx_kill <- which(store$p_true < 0.01 | store$p_true > 0.99)
store <- store[-idx_kill,]
#compute variance per row
var_samp <- store$n_samp*store$p_true*(1-store$p_true)
#convert to visually intuitive form
radius <- sqrt(var_samp/pi)/3
## make plot
par(mfrow=c(2,1))
#first subplot
symbols(x = store$p_true,
y = store$prev,
circles = radius,inches=0.2,
fg="LightBlue",
bg=rgb(red=0,green=0,blue=1,alpha=0.5),
xlim=c(0,1),ylim=c(0,1),
main="Prevalence based estimate",
xlab="True Capability",
ylab="Estimated Capability")
grid()
abline(a=0, b=1, col="Green", lwd=3)
#some annotations
text(0.25,0.85,"Over-estimate zone",pos = 4,bg="white")
text(0.4,0.8,"(ending up here is expensive)",bg="White")
text(0.7,0.15,"Under-estimate zone",pos = 4,bg="white")
text(0.75,0.1,"(ending up here is cheap)",bg="White")
legend(
"topleft",
title = "Uncertainty scaled by 3-sigma",
legend=c("0.3", "1.5", "3.0","6.0","9.0"),
pch = 21,
bty = "n",
col = "black",
pt.bg = "Blue",
pt.cex = c(0.1,0.5,1,2,3)
)
#second subplot
symbols(store$p_true,store$L_CI,radius,inches=0.2,fg="yellow",
bg=rgb(red=1,green=0,blue=0,alpha=0.6),xlim=c(0,1),ylim=c(0,1),
main="Lower CI based estimate",
xlab="True Capability",
ylab="Estimated Capability")
grid()
abline(a=0,b=1,col="Green",lwd=3)
#some annotations
text(0.25,0.85,"Over-estimate zone",pos = 4,bg="white")
text(0.4,0.8,"(ending up here is expensive)",bg="White")
text(0.7,0.15,"Under-estimate zone",pos = 4,bg="white")
text(0.75,0.1,"(ending up here is cheap)",bg="White")
legend(
"topleft",
title = "Uncertainty scaled by 3-sigma",
legend=c("0.3", "1.5", "3.0","6.0","9.0"),
pch = 21,
bty = "n",
col = "black",
pt.bg = "Red",
pt.cex = c(0.1,0.5,1,2,3)
)
par(mfrow=c(1,1))