Problemas con un estudio de simulación de la explicación de los experimentos repetidos de un intervalo de confianza del 95%: ¿dónde me estoy equivocando?

Estoy tratando de escribir un script R para simular la interpretación de experimentos repetidos de un intervalo de confianza del 95%. Descubrí que sobreestima la proporción de veces en que el verdadero valor poblacional de una proporción está contenido dentro del IC del 95% de la muestra. No es una gran diferencia, alrededor del 96% frente al 95%, pero de todos modos esto me interesó.

Mi función toma una muestra samp_nde una distribución de Bernoulli con probabilidad pop_p, y luego calcula un intervalo de confianza del 95% con el prop.test()uso de corrección de continuidad, o más exactamente con binom.test(). Devuelve 1 si la verdadera proporción de la población pop_pestá contenida dentro del IC del 95%. He escrito dos funciones, una que usa prop.test()y otra que usa binom.test()y ha tenido resultados similares con ambas:

in_conf_int_normal <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses normal approximation to calculate confidence interval
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- prop.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2]
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
}
in_conf_int_binom <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses Clopper and Pearson method
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- binom.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2] 
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
 }

Descubrí que cuando repites el experimento unas pocas miles de veces, la proporción de veces que pop_pestá dentro del IC del 95% de la muestra está más cerca de 0,96 en lugar de 0,95.

set.seed(1234)
times = 10000
results <- replicate(times, in_conf_int_binom())
sum(results) / times
[1] 0.9562

Mis pensamientos hasta ahora sobre por qué este puede ser el caso son

• mi código está mal (pero lo he comprobado mucho)
• Inicialmente pensé que esto se debía al problema normal de aproximación, pero luego descubrí binom.test()

¿Alguna sugerencia?

No te vas a equivocar. Simplemente no es posible construir un intervalo de confianza para una proporción binomial que siempre tenga una cobertura de exactamente el 95% debido a la naturaleza discreta del resultado. Se garantiza que el intervalo Clopper-Pearson ('exacto') tiene una cobertura de al menos el 95%. Otros intervalos tienen una cobertura cercana al 95% en promedio , cuando se promedian sobre la proporción real.

Tiendo a favorecer el intervalo Jeffreys, ya que tiene una cobertura cercana al 95% en promedio y (a diferencia del intervalo de puntaje de Wilson) una cobertura aproximadamente igual en ambas colas.

Con solo un pequeño cambio en el código en la pregunta, podemos calcular la cobertura exacta sin simulación.

p <- 0.3
n <- 1000

# Normal test
CI <- sapply(0:n, function(m) prop.test(m,n)$conf.int[1:2])
caught.you <- which(CI[1,] <= p & p <= CI[2,])
coverage.pr <- sum(dbinom(caught.you - 1, n, p))

# Clopper-Pearson
CI <- sapply(0:n, function(m) binom.test(m,n)$conf.int[1:2])
caught.you.again <- which(CI[1,] <= p & p <= CI[2,])
coverage.cp <- sum(dbinom(caught.you.again - 1, n, p))

Esto produce el siguiente resultado.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087