Problemas con un estudio de simulación de la explicación de experimentos repetidos de un intervalo de confianza del 95%: ¿en qué me equivoco?

Question

Estoy tratando de escribir un script en R para simular la interpretación de experimentos repetidos de un intervalo de confianza del 95%. He descubierto que sobrestima 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 aun así me interesó.

Mi función toma una muestra samp_n de una distribución Bernoulli con probabilidad pop_p y luego calcula un intervalo de confianza del 95% con prop.test() utilizando la corrección de continuidad, o más exactamente con binom.test() . Devuelve 1 si la verdadera proporción de población pop_p está contenido en el IC del 95%. He escrito dos funciones, una que utiliza prop.test() y uno que utiliza binom.test() y han tenido resultados similares con ambos:

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)
    }
 }

He comprobado que cuando se repite el experimento unos cuantos miles de veces, la proporción de veces en que el pop_p está dentro del IC del 95% de la muestra está más cerca de 0,96 que 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é esto puede ser el caso son

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

¿Alguna sugerencia?

Answer 1

No te equivocas. Simplemente no es posible construir un intervalo de confianza para una proporción binomial que siempre tiene cobertura de exactamente 95% debido a la naturaleza discreta del resultado. Se garantiza que el intervalo Clopper-Pearson ("exacto") tiene una cobertura de al menos 95%. Otros intervalos tienen una cobertura más cercana al 95%. de media cuando se promedia la proporción real.

Yo mismo me inclino por el intervalo de Jeffreys, ya que tiene una cobertura cercana al 95% de media y (a diferencia del intervalo de puntuación de Wilson) una cobertura aproximadamente igual en ambas colas.

Con sólo un pequeño cambio en el código de la pregunta, podemos calcular la cobertura exacta sin necesidad de 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))

El resultado es el siguiente.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087