Probleme mit einer Simulationsstudie der wiederholten Experimente Erklärung eines 95% -Konfidenzintervalls - wo mache ich Fehler?

Ich versuche, ein R-Skript zu schreiben, um die wiederholte Interpretation eines 95% -Konfidenzintervalls zu simulieren. Ich habe festgestellt, dass es den Anteil der Zeiten überschätzt, in denen der wahre Populationswert eines Anteils im 95% -KI der Stichprobe enthalten ist. Kein großer Unterschied - ungefähr 96% gegenüber 95%, aber das hat mich trotzdem interessiert.

Meine Funktion entnimmt samp_nmit Wahrscheinlichkeit eine Stichprobe aus einer Bernoulli-Verteilung pop_pund berechnet dann ein 95% -Konfidenzintervall unter prop.test()Verwendung der Kontinuitätskorrektur oder genauer mit binom.test(). Es wird 1 zurückgegeben, wenn der wahre Bevölkerungsanteil pop_pim 95% -KI enthalten ist. Ich habe zwei Funktionen geschrieben, eine, die verwendet, prop.test()und eine, die binom.test()ähnliche Ergebnisse verwendet und mit beiden erzielt hat:

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

Ich habe festgestellt, dass, wenn Sie das Experiment einige tausend Mal wiederholen, der Anteil der Fälle, in denen der Wert pop_pinnerhalb des 95% -KI der Probe liegt, eher bei 0,96 als bei 0,95 liegt.

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

Meine bisherigen Gedanken darüber, warum dies der Fall sein könnte, sind

• Mein Code ist falsch (aber ich habe ihn oft überprüft)
• Ich dachte zunächst, dass dies auf das normale Approximationsproblem zurückzuführen ist, fand es dann aber binom.test()

Irgendwelche Vorschläge?

Du machst nichts falsch. Es ist einfach nicht möglich, ein Konfidenzintervall für einen Binomialanteil zu erstellen, der aufgrund der diskreten Natur des Ergebnisses immer eine Abdeckung von genau 95% aufweist. Das Clopper-Pearson-Intervall ("genau") hat garantiert eine Abdeckung von mindestens 95%. In anderen Intervallen liegt die Abdeckung im Durchschnitt näher bei 95% , wenn sie über den tatsächlichen Anteil gemittelt wird.

Ich tendiere dazu, das Jeffreys-Intervall selbst zu bevorzugen, da es im Durchschnitt eine Abdeckung von fast 95% aufweist und (im Gegensatz zum Wilson-Score-Intervall) in beiden Schwänzen ungefähr die gleiche Abdeckung aufweist.

Mit nur einer kleinen Änderung des Codes in der Frage können wir die genaue Abdeckung ohne Simulation berechnen.

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

Dies ergibt die folgende Ausgabe.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087