Szacowanie dwumianowego przedziału ufności - dlaczego nie jest symetryczne?

30

Użyłem następującego kodu r do oszacowania przedziałów ufności proporcji dwumianowej, ponieważ rozumiem, że zastępuje to „obliczanie mocy” podczas projektowania krzywych charakterystycznych dla odbiornika działającego na wykrywanie chorób w populacji.

n wynosi 150, a choroba, naszym zdaniem, występuje w populacji w 25%. Obliczyłem wartości 75% czułości i 90% swoistości (ponieważ wydaje się, że ludzie tak robią).

    binom.test(c(29,9), p=0.75, alternative=c("t"), conf.level=0.95)

    binom.test(c(100, 12), p=0.90, alternative=c("t"), conf.level=0.95)

Odwiedziłem również tę stronę:

http://statpages.org/confint.html

Która jest stroną Java, która oblicza dwumianowe przedziały ufności i daje tę samą odpowiedź.

W każdym razie po tej długiej konfiguracji chcę zapytać, dlaczego przedziały ufności nie są symetryczne, np. Czułość jest

   95 percent confidence interval:
   0.5975876 0.8855583 

   sample estimate probability: 0.7631579

Przepraszam, jeśli to głupie pytanie, ale wszędzie, na co patrzę, wydaje się sugerować, że będą symetryczne, a mój kolega wydaje się, że też będą.

confidence-interval binomial

— Chris Beeley
źródło

20

Uważa się je za symetryczne, ponieważ dość często stosuje się normalne przybliżenie. Ten działa wystarczająco dobrze, jeśli p wynosi około 0,5. binom.testz drugiej strony podaje „dokładne” przedziały Cloppera-Pearsona, które są oparte na rozkładzie F (zobacz tutaj dokładne formuły obu podejść). Gdybyśmy zaimplementowali interwał Cloppera-Pearsona w R, byłoby to coś w rodzaju (patrz uwaga ):

Clopper.Pearson <- function(x, n, conf.level){
    alpha <- (1 - conf.level) / 2
    QF.l <- qf(1 - alpha, 2*n - 2*x + 2, 2*x)
    QF.u <- qf(1 - alpha, 2*x + 2, 2*n - 2*x)

    ll <- if (x == 0){
          0
    } else { x / ( x + (n-x+1)*QF.l ) }

    uu <- if (x == 0){
          0
    } else { (x+1)*QF.u / ( n - x + (x+1)*QF.u ) }

    return(c(ll, uu))
}

Zarówno w linku, jak i implementacji widać, że formuła dla górnej i dolnej granicy jest zupełnie inna. Jedynym przypadkiem symetrycznego przedziału ufności jest sytuacja, gdy p = 0,5. Używając wzorów z linku i biorąc pod uwagę, że w tym przypadku łatwo jest ustalić, jak to jest. $n = 2\times x$

Osobiście zrozumiałem to lepiej, patrząc na przedziały ufności oparte na podejściu logistycznym. Dane dwumianowe są na ogół modelowane za pomocą funkcji łącza logit, zdefiniowanej jako:

l o g i t (x) = \log (\frac{x}{1 - x})

${\rm logit}(x) = \log\! \bigg( \frac{x}{1-x} \bigg)$

Ta funkcja łącza „odwzorowuje” błąd w regresji logistycznej na rozkład normalny. W konsekwencji przedziały ufności w ramach logistycznych są symetryczne wokół wartości logit, podobnie jak w klasycznych ramach regresji liniowej. Transformacja logit jest używana dokładnie, aby umożliwić wykorzystanie całej teorii opartej na normalności wokół regresji liniowej.

Po wykonaniu odwrotnej transformacji:

{l o g i t}^{- 1} (x) = \frac{e^{x}}{1 + e^{x}}

${\rm logit}^{-1}(x) = \frac{e^x}{1+e^{x}}$

Znowu dostajesz asymetryczny interwał. Teraz te przedziały ufności są właściwie stronnicze. Ich zasięg nie jest tym, czego można się spodziewać, zwłaszcza na granicach rozkładu dwumianowego. Jednak jako ilustracja pokazują, dlaczego logiczne jest, że rozkład dwumianowy ma asymetryczne przedziały ufności.

Przykład w R:

logit <- function(x){ log(x/(1-x)) }
inv.logit <- function(x){ exp(x)/(1+exp(x)) }
x <- c(0.2, 0.5, 0.8)
lx <- logit(x)
upper <- lx + 2
lower <- lx - 2

logxtab <- cbind(lx, upper, lower)
logxtab # the confidence intervals are symmetric by construction
xtab <- inv.logit(logxtab)
xtab # back transformation gives asymmetric confidence intervals

Uwaga : W rzeczywistości R używa rozkładu beta, ale jest to całkowicie równoważne i obliczeniowo nieco bardziej wydajne. Implementacja w R różni się zatem od tego, co tutaj pokazuję, ale daje dokładnie ten sam rezultat.

— Joris Meys
źródło

2

Czy naprawdę chciałeś powiedzieć, że logit „przekształca rozkład dwumianowy w rozkład normalny”?

— whuber

@ whuber: niezły połów formuły i niezły połów formuły. Prawie nie. Zapewnia to, że błędy w regresji logistycznej są zgodne z rozkładem normalnym. Dziękuję za korektę.

— Joris Meys,

Krótka uwaga techniczna, transformacja „łukowa” to taka, która ma szybszą konwergencję do normalności niż transformacja logistyczna. Ustaw

(gdzie

jest liczbą „sukcesów”, a

liczbą prób), a przy pomocy tak zwanej „metody delta” można wykazać, że wariancja

jest w przybliżeniu stała (i niezależna od

, jak powinna być w normalna dystrybucja).

Y = \frac{2}{π} \arcsin \sqrt{\frac{X}{N}}

$Y=\frac{2}{\pi}\arcsin{\sqrt{\frac{X}{N}}}$

X

$X$

N

$N$

Y

$Y$

Y

$Y$

— prawdopodobieństwo prawdopodobieństwa

Link podany dla „dokładnych prawdopodobieństw” jest zepsuty. Czy masz jeszcze jeden?

— S. Kolassa - Przywróć Monikę

@StephanKolassa Można również znaleźć tutaj formuły Cloppera Pearsona: en.wikipedia.org/wiki/…

— Joris Meys,

24

Aby zobaczyć, dlaczego nie powinno być symetryczne, pomyśl o sytuacji, w której a otrzymasz 9 sukcesów w 10 próbach. Następnie i 95% przedział ufności dla wynosi [0.554, 0.997]. Górna granica nie może być większa niż 1 oczywiście, więc większość niepewność może spaść z lewej . $p=0.9$ $\hat{p}=0.9$ $p$ $\hat{p}$

— Rob Hyndman
źródło

9

@Joris wspomniał o symetrycznym lub „asymptotycznym” przedziale, który jest najprawdopodobniej tym, którego się spodziewasz. @Joris wspomniał również o „dokładnych” odstępach Cloppera-Pearsona i dał ci referencje, które wyglądają bardzo ładnie. Istnieje inny przedział ufności dla proporcji, które prawdopodobnie napotkasz (pamiętaj, że nie jest on również symetryczny), przedział „Wilsona”, który jest rodzajem asymptotycznego przedziału opartego na odwróceniu wyniku testu. Końce przedziału rozwiązania (w ) Równanie $p$

(\hat{p} - p) / \sqrt{p (1 - p)} = \pm z_{α / 2}

$(\hat{p} - p)/\sqrt{p(1-p)}=\pm z_{\alpha/2}$

W każdym razie, możesz dostać wszystkie trzy w R z następującymi:

library(Hmisc)
binconf(29, 38, method = "asymptotic")
binconf(29, 38, method = "exact")
binconf(29, 38, method = "wilson")

Zauważ, że metoda „wilson” jest tym samym przedziałem ufności, który jest używany przez prop.test bez korekcji ciągłości Yatesa:

prop.test(29, 38, correct = FALSE)

Zobacz tutaj darmowy podręcznik Laury Thompson SPLUS + R, który towarzyszy Kategorycznej analizie danych Agresti, w której kwestie te są szczegółowo omówione.

1

(+1) Fajnie, że cytujesz podręcznik Laury i dodajesz ten zestaw informacji o CI Wilsona.

— chl

2

Dzięki. Chciałbym zauważyć, że interwał Wilsona omówiono w artykule, do którego odwołuje się @Joris.

9

Tam są symetryczne przedziały ufności dla rozkładu dwumianowego: asymetria nie jest zmuszony na nas, mimo wszystkich powodów już wymienionych. Przedziały symetryczne są zwykle uważane za gorsze

Chociaż są one symetryczne numerycznie , prawdopodobieństwo nie jest symetryczne : to znaczy, że ich jednostronne pokrycia różnią się od siebie. To - niezbędna konsekwencja możliwej asymetrii rozkładu dwumianowego - jest sednem sprawy.
Często jeden punkt końcowy musi być nierealny (mniejszy niż 0 lub większy niż 1), jak wskazuje @Rob Hyndman.

Powiedziawszy to, podejrzewam, że liczbowo symetryczne CI mogą mieć pewne dobre właściwości, takie jak tendencja do bycia krótszymi niż probabilistycznie symetryczne w niektórych okolicznościach.

— Whuber
źródło

\hat{p} = k / n

$\hat p = k/n$

@cb Nie podążam za tym. Po pierwsze, najkrótszy CI niekoniecznie musi mieć jednakową gęstość na każdym końcu. Po drugie, komentarz na temat „nie istnieje” nie ma dla mnie sensu: co oznacza „nie istnieje”?

— whuber

1

najkrótszy CI. Aby obliczyć najkrótszy CI dla danego zasięgu, zacznę od maksymalnej gęstości i powiększę krótki krok w stronę, gdzie gęstość jest większa. Tam uzyskuję największą pewność siebie (to znaczy krótki krok). Powiększam ci wielokrotnie, aż mam pożądany obszar (zasięg). Jeśli moje kroki są małe (nieskończenie małe), gęstość po obu stronach będzie (w przybliżeniu) taka sama. Czy popełniłem błąd w tej strategii?

— cbeleites obsługuje Monikę

p

$p$

\hat{p} = 4 / 5 = 0.8

$\hat p = 4/5 = 0.8$

p

$p$

p

$p$

< 70 %

$< 70 \%$

— cbeleites wspiera Monikę

1

p = 0.8

$p = 0.8$

k \in {3, 4, 5}

$k \in \{3, 4, 5\}$

n = 5

$n = 5$

p

$p$

n

$n$

k

$k$

p

$p$

k = 4

$k = 4$

n = 5

$n = 5$

P r (p | n = 5, k = 4)

$Pr (p | n = 5, k = 4)$

p \in [0, 1]

$p \in [0, 1]$

P r (k | n, p)

$Pr (k | n, p)$

p

$p$

6

$p$ $0$ $1$ $n$ $p\approx 0.5$

— chl
źródło

2

Wiem, że minęło trochę czasu, ale pomyślałem, że będę tu dzwonił. Biorąc pod uwagę nip, łatwo jest obliczyć prawdopodobieństwo określonej liczby sukcesów bezpośrednio przy użyciu rozkładu dwumianowego. Następnie można sprawdzić rozkład, aby zobaczyć, czy nie jest on symetryczny. Zbliża się do symetrii dla dużych np i dużych n (1-p).

Można zebrać prawdopodobieństwa w ogonach, aby obliczyć konkretny CI. Biorąc pod uwagę dyskretny charakter rozkładu, znalezienie szczególnego prawdopodobieństwa w ogonie (np. 2,5% dla 95% CI) będzie wymagało interpolacji między liczbą sukcesów. Dzięki tej metodzie można obliczać CI bezpośrednio bez przybliżenia (innego niż wymagana interpolacja).

— Dr Eric
źródło