Problemy z badaniem symulacyjnym powtórzonych eksperymentów z wyjaśnieniem 95% przedziału ufności - gdzie się mylę?

9

Próbuję napisać skrypt R, aby zasymulować interpretację powtarzanych eksperymentów z 95% przedziałem ufności. Przekonałem się, że przecenia on odsetek przypadków, w których prawdziwa wartość populacyjna części jest zawarta w 95% CI próby. Niewielka różnica - około 96% vs 95%, ale to mnie jednak zainteresowało.

Moja funkcja pobiera próbkę samp_nz rozkładu Bernoulliego z prawdopodobieństwem pop_p, a następnie oblicza 95% przedział ufności za prop.test()pomocą korekcji ciągłości, a dokładniej za pomocą binom.test(). Zwraca 1, jeśli prawdziwy odsetek populacji pop_pzawiera się w 95% CI. Napisałem dwie funkcje, jedną, która używa prop.test()i jedną, która używa binom.test()i przyniosła podobne wyniki dla obu:

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

Przekonałem się, że kiedy powtórzysz eksperyment kilka tysięcy razy, odsetek przypadków, gdy pop_pmieści się w 95% CI próbki, jest bliższy 0,96 niż 0,95.

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

Moje dotychczasowe przemyślenia na temat tego, dlaczego tak się dzieje

  • mój kod jest nieprawidłowy (ale często go sprawdzałem)
  • Początkowo myślałem, że to z powodu normalnego problemu z przybliżeniem, ale potem znalazłem binom.test()

Jakieś sugestie?

Andrzej
źródło
(+1) Przy okazji, ponownie uruchomiłem twój kod times=100000kilka razy i zobaczyłem ten sam wynik. Jestem ciekawy, czy ktoś ma na to wytłumaczenie. Kod jest wystarczająco prosty, że jestem prawie pewien, że nie ma błędu kodowania. Ponadto jeden przebieg z times=1000000podał .954931jako wynik.
Makro
3
(+1) Ale dlaczego spodziewasz się uzyskać dokładnie 95%? Na przykład Clopper Pearson ma gwarancję zachowania ostrożności. Dla Państwa i , mam, że CI powinny obejmować rzeczywistej wartości 95.3648% czasu. np
kardynał
2
Aby poprzeć komentarz kardynałów, dokładne dwumianowe prawdopodobieństwa są dokładne, ponieważ są oparte na dokładnym obliczeniu prawdopodobieństwa, ale niekoniecznie dają dokładny poziom ufności. Jest tak, ponieważ dwumian jest rozkładem dyskretnym. Tak więc Clopper-Pearson wybiera punkt końcowy dla tego przedziału, aby uzyskać najbliższe prawdopodobieństwo do poziomu ufności na poziomie powyżej lub powyżej. Tworzy to również piłokształtne zachowanie funkcji mocy dokładnego testu dwumianowego. Ten dziwny, ale podstawowy wynik omówiono w mojej pracy z Christine Liu w American Statistician (2002).
Michael R. Chernick,
1
Szczegóły mojej pracy pod tym linkiem: citeulike.org/user/austin987/article/7571878
Michael R.
3
Dokładne dwumianowe elementy CI są „dokładne”, ponieważ ich rzeczywista wydajność jest równa ich nominalnej wydajności, a nie dlatego, że obliczenie prawdopodobieństwa jest „dokładne”! Należy rozumieć, że CI CI musi mieć co najmniej szansę na pokrycie prawdziwego parametru, bez względu na to, jaki jest rozkład podstawowy (w obrębie założonej rodziny). „Dokładne” oznacza, że ​​minimum wszystkich tych zasięgów, przejętych przez całą rodzinę dystrybucji, wynosi . Aby to osiągnąć, faktyczny zasięg wielu możliwych dystrybucji często musi być większy niż . 1-α1-α 1-α1-α
whuber

Odpowiedzi:

9

Nie pomylisz się. Po prostu niemożliwe jest zbudowanie przedziału ufności dla proporcji dwumianowej, która zawsze ma dokładnie 95% pokrycie z powodu dyskretnej natury wyniku. Gwarantowany interwał Cloppera-Pearsona („dokładny”) wynosi co najmniej 95%. Inne przedziały mają zasięg zbliżony średnio do 95% , gdy uśrednia się je względem rzeczywistej proporcji.

Ja sam preferuję interwał Jeffreysa, ponieważ ma on zasięg średnio blisko 95% i (w przeciwieństwie do przedziału punktacji Wilsona) w przybliżeniu jednakowy zasięg w obu ogonach.

Przy niewielkiej zmianie kodu w pytaniu możemy obliczyć dokładny zasięg bez symulacji.

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

Daje to następujące dane wyjściowe.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087
jeden przystanek
źródło
1
Po prostu nie jest możliwe zbudowanie przedziału ufności dla proporcji dwumianowej, która zawsze ma pokrycie dokładnie 95% ze względu na dyskretny charakter wyniku ” - być może, pomijając (nieco dziwną) możliwość losowych przedziałów . (Przynajmniej w ten sposób można to zrobić, choć może się zdarzyć, że zwykle nie powinno .)
Glen_b
2
@Glen_b Od dawna ciekawi mnie sprzeciw wobec losowych decyzji. Uważam, że Jack Kiefer zauważył, że jeśli jesteś w porządku przy użyciu randomizacji do pobierania próbek, nie powinieneś mieć problemu z wykorzystaniem jej w procesie decyzyjnym. Jeśli potrzebujesz procedury decyzyjnej, która jest odtwarzalna, udokumentowana i trudna do oszukiwania, po prostu wygeneruj dowolne wartości potrzebne dla losowego interwału przed zebraniem danych - włącz je do projektu.
whuber