Przedział ufności dla mediany

40

Muszę znaleźć 95% CI na medianie i innych percentylach. Nie wiem jak do tego podejść. Głównie używam R jako narzędzia programistycznego.

r confidence-interval median Dominic Comtois
źródło

Odpowiedzi:

31

Oto ilustracja klasycznego zestawu danych R.

> x       = faithful$waiting
> bootmed = apply(matrix(sample(x, rep=TRUE, 10^4*length(x)), nrow=10^4), 1, median)
> quantile(bootmed, c(.025, 0.975))
2.5% 97.5% 
 73.5    77

co daje przedział ufności (73,5; 77) dla mediany.

( Uwaga: Poprawiona wersja, dzięki Jana . Kiedyś w wcześniej, co doprowadziło do zamieszania!)103nrow

Xi'an
źródło
7
Wydaje mi się podejrzanie wąski. Użycie funkcji z library(boot)wydaje się to potwierdzić:> boot.ci (boot (x, funkcja (x, i)) mediana (x [i]), R = 1000)) Interwały: Poziom Normalny Podstawowy 95% (74,42; 78,22) (75,00 , 78,49) Poziom Percentile BCa 95% (73,51, 77,00) (73,00, 77,00)
onestop
2
nie ma za co Xi'an ... Nawiasem mówiąc, zawsze wolę ustawić oryginalną wartość N w macierzy, ponieważ jest to stała dla różnych rozmiarów bootstrap, które mogę zrobić. Więc zazwyczaj powiedziałbym, że ncol = długość (x). Uważam, że w ten sposób jest mniej szansa na błąd.
John
6
Jest to po prostu nieefektywny sposób obliczenia kwantyli dwumianowych jak w odpowiedzi onestop .
whuber
30

Inne podejście opiera się na kwantylach rozkładu dwumianowego.
na przykład:

> x=faithful$waiting
> sort(x)[qbinom(c(.025,.975), length(x), 0.5)]
[1] 73 77
jeden przystanek
źródło
4
Podoba mi się prostota tego ... Wyniki są zbliżone do metody bootstrap.
Dominic Comtois
1
Jest to oczywiście o wiele bardziej wydajne niż ładowanie początkowe dla ciągłego przypadku, ale jedną wadą jest to, że nie uwzględnia powiązanych szeregów. Czy zdarza ci się wiedzieć o obejściu tego problemu?
ali_m
15

Sprawdź resampling bootstrap. Wyszukaj w pomocy R funkcję rozruchu. W zależności od danych z ponownym próbkowaniem możesz oszacować przedziały ufności dla prawie wszystkiego.

Tharen
źródło
Zgodzić się. To jest najlepsze podejście. Moim zdaniem niedostatecznie wykorzystywany w naukach biomedycznych.
pmgjones
10
Rozważ spojrzenie na wygładzony pasek startowy do szacowania kwantyli populacji, ponieważ wydaje się, że w tym przypadku występują problemy z konwencjonalnym boostrap - odniesienia można znaleźć w tym pliku pdf . Jeśli interesowała Cię tylko teoretyczna mediana, można zastosować estymator Hodgesa-Lehmana - podany np. Przez wilcox.test(..., conf.int=TRUE)funkcję R.
caracal
4

Są też inne podejścia: jedno oparte jest na teście sumy rang Wilcoxona zastosowanym dla jednej próbki z korektą ciągłości. W R można to podać jako:

wilcox.test(x,conf.level=0.95,alternative="two.sided",correct=TRUE)

I tutaj jest CI CI dla mediany Davida Olive'a:

CI dla mediany

Germaniawerks
źródło
1

Wynik oparty na metodzie qbinom jest nieprawidłowy dla małych próbek. Załóżmy, że x ma 10 składników. Następnie qbinom (c (.025; .975), 10, .5) daje 2 i 8. Wynikowy interwał nie traktuje statystyk rzędu dolnego ogona symetrycznie z tymi z górnego ogona; powinieneś dostać 2 i 9 lub 3 i 8. Prawidłowa odpowiedź to 2 i 9. Możesz sprawdzić w proc univariate w SAS. Złap tutaj, nie potrzebujesz więcej niż 0,025 prawdopodobieństwa poniżej i powyżej; dolny kwantyl tego nie robi, ponieważ daje co najmniej 0,025 co najmniej na poziomie. Zostajesz zapisany na dole, ponieważ liczba, która powinna wynosić 1, powinna zostać odwzorowana na statystykę drugiego rzędu, licząc 0, więc „wyłączanie o jeden” anuluje się. To przypadkowe anulowanie nie zdarza się na górze, więc tutaj otrzymujesz złą odpowiedź. Kod sort (x) [qbinom (c (.025, .975), długość (x) ,. 5) + c (0,1)] prawie działa, a .5 można zastąpić innymi wartościami kwantyli, aby uzyskać przedziały ufności dla innych kwantyli, ale nie będzie dobrze, gdy istnieje takie, że P [X <= a ] =. 025. Patrz np. Higgins, statystyki nieparametryczne.

John Kolassa
źródło