Problèmes avec une étude de simulation de l'explication des expériences répétées d'un intervalle de confiance à 95% - où vais-je mal?

J'essaie d'écrire un script R pour simuler l'interprétation des expériences répétées d'un intervalle de confiance à 95%. J'ai trouvé qu'il surestime la proportion de fois où la valeur réelle de la population d'une proportion est contenue dans l'IC à 95% de l'échantillon. Pas une grande différence - environ 96% vs 95% mais cela m'intéresse néanmoins.

Ma fonction prend un échantillon samp_nd'une distribution de Bernoulli avec probabilité pop_p, puis calcule un intervalle de confiance à 95% en prop.test()utilisant la correction de continuité, ou plus exactement avec binom.test(). Il renvoie 1 si la véritable proportion de la population pop_pest contenue dans l'IC à 95%. J'ai écrit deux fonctions, une qui utilise prop.test()et une qui utilise binom.test()et a eu des résultats similaires avec les deux:

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

J'ai constaté que lorsque vous répétez l'expérience plusieurs milliers de fois, la proportion de fois où le pop_pest dans l'IC à 95% de l'échantillon est plus proche de 0,96 plutôt que de 0,95.

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

Jusqu'à présent, mes réflexions sur les raisons pour lesquelles cela peut être le cas

• mon code est faux (mais je l'ai beaucoup vérifié)
• J'ai d'abord pensé que cela était dû au problème d'approximation normal, mais j'ai trouvé binom.test()

Aucune suggestion?

Vous ne vous trompez pas. Il n'est tout simplement pas possible de construire un intervalle de confiance pour une proportion binomiale qui a toujours une couverture d' exactement 95% en raison de la nature discrète du résultat. L'intervalle Clopper-Pearson («exact») est garanti d'avoir une couverture d' au moins 95%. D'autres intervalles ont une couverture plus proche de 95% en moyenne , lorsqu'ils sont moyennés sur la vraie proportion.

J'ai tendance à favoriser moi-même l'intervalle de Jeffreys, car il a une couverture proche de 95% en moyenne et (contrairement à l'intervalle de score de Wilson) une couverture à peu près égale dans les deux queues.

Avec seulement un petit changement dans le code de la question, nous pouvons calculer la couverture exacte sans simulation.

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

Cela donne la sortie suivante.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087