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

9

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?

Andrew
quelle
(+1) Übrigens habe ich Ihren Code times=100000ein paar Mal wiederholt und das gleiche Ergebnis gesehen. Ich bin gespannt, ob jemand eine Erklärung dafür hat. Der Code ist so einfach, dass ich mir ziemlich sicher bin, dass kein Codierungsfehler vorliegt. Auch ein Lauf mit times=1000000gab .954931als Ergebnis.
Makro
3
np
2
Zur Unterstützung von Kardinalskommentaren sind exakte Binomialwahrscheinlichkeiten genau, da sie auf einer exakten Wahrscheinlichkeitsberechnung basieren, aber nicht unbedingt das genaue Konfidenzniveau angeben. Das liegt daran, dass das Binom eine diskrete Verteilung ist. Daher wählt Clopper-Pearson den Endpunkt für das Intervall so aus, dass Sie die Wahrscheinlichkeit haben, die dem Konfidenzniveau bei oder darüber am nächsten kommt. Dies erzeugt auch ein Sägezahnverhalten für die Potenzfunktion eines exakten Binomialtests. Dieses seltsame, aber grundlegende Ergebnis wird in meiner Arbeit mit Christine Liu in der American Statistician (2002) diskutiert.
Michael R. Chernick
1
Details zu meinem Artikel unter folgendem
Michael R. Chernick
3
1α1α 1α1α
whuber

Antworten:

9

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
ein Stop
quelle
1
" Es ist einfach nicht möglich, ein Konfidenzintervall für einen Binomialanteil zu konstruieren, der aufgrund der diskreten Natur des Ergebnisses immer eine Abdeckung von genau 95% aufweist " - abgesehen vielleicht von der (etwas merkwürdigen) Möglichkeit randomisierter Intervalle . (Zumindest auf diese Weise kann es getan werden, obwohl es gut sein kann, dass es normalerweise nicht sollte .)
Glen_b -State Monica
2
@Glen_b Ich war lange neugierig auf die Einwände gegen zufällige Entscheidungen. Ich glaube, Jack Kiefer bemerkte, dass Sie kein Problem damit haben sollten, Ihre Stichproben im Entscheidungsprozess zu verwenden, wenn Sie in Ordnung sind, Randomisierung zum Sammeln Ihrer Proben zu verwenden. Wenn Sie ein Entscheidungsverfahren benötigen, das reproduzierbar, dokumentiert und schwer zu betrügen ist, generieren Sie einfach alle für das zufällige Intervall erforderlichen Zufallswerte, bevor Sie die Daten erfassen - machen Sie es Teil des Entwurfs.
whuber