Ocena mocy testu normalności (w R)

Chcę ocenić dokładność testów normalności dla różnych wielkości próbek w R (zdaję sobie sprawę, że testy normalności mogą być mylące ). Na przykład, aby spojrzeć na test Shapiro-Wilka, przeprowadzam następującą symulację (a także sporządzam wyniki) i oczekuję, że wraz ze wzrostem wielkości próby maleje prawdopodobieństwo odrzucenia wartości zerowej:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

Myślałem, że wraz ze wzrostem wielkości próby powinna być mniejsza częstotliwość odrzucania, jednak wydaje się ona dość jednolita. Myślę, że nie rozumiem tego - wszelkie myśli są mile widziane.

Symulujesz w ramach hipotezy zerowej (rozkład normalny), dlatego wskaźnik odrzucenia będzie dążył do poziomu istotności zgodnie z oczekiwaniami. Aby ocenić moc, musisz przeprowadzić symulację w dowolnym niestandardowym rozkładzie. Istnieją nieskończone możliwości / scenariusze (np. Rozkłady gamma ze wzrostem skośności, rozkład t ze spadkiem df itp.) Do wyboru, w zależności od zakresu badania.

Zrozumienie analizy mocy testów hipotez statystycznych można poprawić, przeprowadzając je i uważnie przyglądając się wynikom.

Z założenia test wielkości $\alpha$ma na celu odrzucenie hipotezy zerowej z co najmniej szansą$\alpha$gdy wartość null jest prawdziwa (oczekiwany współczynnik fałszywie dodatnich wyników ). Gdy mamy możliwość (lub luksusu) wyboru spośród alternatywnych procedur z tą właściwością, wolelibyśmy te, które (a) faktycznie zbliżają się do nominalnego wskaźnika fałszywie dodatnich i (b) mają stosunkowo większe szanse na odrzucenie hipotezy zerowej, gdy jest ona nie prawda.

Drugie kryterium wymaga od nas ustalenia, w jaki sposób (y) i o ile zero nie jest prawdziwe. W przypadkach podręczników jest to łatwe, ponieważ alternatywy mają ograniczony zakres i są wyraźnie określone. W przypadku testów dystrybucji, takich jak Shapiro-Wilk, alternatywy są znacznie bardziej niejasne: są „nienormalne”. Wybierając spośród testów dystrybucji, analityk najprawdopodobniej będzie musiał przeprowadzić własne jednorazowe badanie mocy, aby ocenić, jak dobrze testy działają w odniesieniu do bardziej szczegółowych alternatywnych hipotez, które budzą obawy w omawianym problemie.

Przykład motywowany odpowiedzią Michaela Mayera zakłada, że ​​rozkład alternatywny może mieć cechy podobne do cech rodziny rozkładów Studenta. Ta rodzina, sparametryzowana liczbą$\nu\ge 1$ (a także według lokalizacji i skali) obejmuje limit dużych $\nu$ rozkład normalny.

W obu sytuacjach - niezależnie od tego, czy oceniamy rzeczywisty rozmiar testu, czy jego moc - musimy wygenerować niezależne próbki z określonego rozkładu, uruchomić test na każdej próbce i znaleźć szybkość, z jaką odrzuca hipotezę zerową. Jednak w każdym wyniku testu dostępnych jest więcej informacji: jego wartość P. Zachowując zestaw wartości P wytworzonych podczas takiej symulacji, możemy później ocenić szybkość, z jaką test odrzuciłby wartość zerową dla dowolnej wartości$\alpha$możemy się tym przejmować. Sercem analizy mocy jest zatem podprogram, który generuje taki rozkład wartości P (albo przez symulację, jak właśnie opisano, albo - czasami - z formułą teoretyczną). Oto zakodowany przykład R. Jego argumenty obejmują

• rdist, nazwa funkcji do generowania losowej próbki z pewnego rozkładu

• n, wielkość próbek na żądanie rdist

• n.iter, liczba takich próbek do uzyskania

• ..., wszelkie opcjonalne parametry, które należy przekazać rdist(takie jak stopnie swobody$\nu$).

Pozostałe parametry kontrolują wyświetlanie wyników; zostały one uwzględnione głównie w celu ułatwienia generowania liczb w tej odpowiedzi.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

Widać, że obliczenia zajmują tylko jedną linię; reszta kodu wykreśla histogram. Aby to zilustrować, użyjmy go do obliczenia oczekiwanych fałszywie dodatnich wskaźników. „Stawki” są w liczbie mnogiej, ponieważ właściwości testu zwykle różnią się w zależności od wielkości próbki. Ponieważ dobrze wiadomo, że testy dystrybucyjne mają dużą moc w porównaniu z jakościowo małymi alternatywami, gdy wielkości próbek są duże, niniejsze badanie koncentruje się na szeregu małych rozmiarów próbek, w których takie testy często są stosowane w praktyce: zazwyczaj około$5$ do $100.$ Aby zaoszczędzić czas obliczeń, raportuję tylko o wartościach $n$ od $5$ do $20.$

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

Po określeniu parametrów ten kod jest również tylko jedną linią. Daje następujące dane wyjściowe:

Jest to oczekiwany wygląd: histogramy pokazują prawie jednolite rozkłady wartości P w pełnym zakresie od$0$ do $1$. Przy nominalnym rozmiarze ustawionym na$\alpha=0.05,$ raport z symulacji między $.0481$ i $0.0499$wartości P były w rzeczywistości mniejsze niż ten próg: są to wyniki wyróżnione na czerwono. Bliskość tych częstotliwości do wartości nominalnej potwierdza, że ​​test Shapiro-Wilka wykonuje się zgodnie z reklamą.

(Wydaje się, że istnieje tendencja do niezwykle wysokiej częstotliwości wartości P w pobliżu $1$. Nie ma to większego znaczenia, ponieważ w prawie wszystkich aplikacjach jedyne wartości P, na które się patrzy, to$0.2$ lub mniej.)

Przejdźmy teraz do oceny mocy. Pełny zakres wartości$\nu$ dla Studenta rozkład t można odpowiednio zbadać, oceniając kilka przykładów z całego $\nu=100$ aż do $\nu=1$. Skąd to wiem? Wykonałem kilka wstępnych testów, używając bardzo małej liczby iteracji (z$100$ do $1000$), co wcale nie zajmuje czasu. Kod wymaga teraz podwójnej pętli (aw bardziej skomplikowanych sytuacjach często potrzebujemy potrójnych lub poczwórnych pętli, aby uwzględnić wszystkie aspekty, które musimy zmieniać): jeden do badania, w jaki sposób moc zmienia się w zależności od wielkości próbki, a drugi do badania, w jaki sposób zmienia się w zależności od stopnie swobody. Jednak po raz kolejny wszystko odbywa się w jednym wierszu kodu (trzeci i ostatni):

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Krótkie studium tego obrazu zapewnia dobrą intuicję na temat mocy. Chciałbym zwrócić uwagę na jego najbardziej istotne i przydatne aspekty:

• W miarę zmniejszania się stopni swobody z $\nu=64$ po lewej stronie do $\nu=1$po prawej stronie, coraz więcej wartości P jest małych, co pokazuje, że wzrasta moc odróżniania tych rozkładów od rozkładu Normalnego. (Moc jest określana ilościowo w każdym tytule wykresu: jest równa proporcji czerwonego pola histogramu).

• W miarę wzrostu wielkości próbki od $n=5$ w górnym rzędzie do $n=20$ na dole również rośnie moc.

• Zwróć uwagę, że ponieważ rozkład alternatywny różni się bardziej od rozkładu zerowego, a wielkość próbki rośnie, wartości P zaczynają się zbierać w lewo, ale nadal istnieje ich „ogon” rozciągający się aż do $1$. Jest to charakterystyczne dla badań mocy. Pokazuje, że testowanie jest ryzykowne : nawet gdy hipoteza zerowa zostanie rażąco naruszona, a nawet gdy nasza próbka jest dość duża, nasz test formalny może nie dać znaczącego wyniku.

• Nawet w skrajnym przypadku w prawym dolnym rogu, gdzie próbka $20$ pochodzi z rozkładu t Studenta z $1$ stopień swobody (rozkład Cauchy'ego), moc nie jest $1$: tam jest $100 - 86.57 = 13\%$ szansa, że ​​próbka $20$ Warianty iid Cauchy nie będą uważane za znacząco różne od normalnych na poziomie $5\%$ (to znaczy z $95\%$ pewność siebie).

• Możemy ocenić moc przy dowolnej wartości $\alpha$wybieramy kolorując więcej lub mniej pasków na tych histogramach. Na przykład, aby ocenić moc przy$\alpha=0.10$, pokoloruj dwa lewe słupki na każdym histogramie i oszacuj jego powierzchnię jako ułamek całości.

  (To nie zadziała zbyt dobrze dla wartości $\alpha$ mniejszy niż $0.05$z tą postacią. W praktyce histogramy ograniczałyby się do wartości P tylko w zakresie, który zostałby zastosowany, być może od$0$ do $20\%$i pokaż je na tyle szczegółowo, aby umożliwić wizualną ocenę mocy do $\alpha=0.01$ lub nawet $\alpha=0.005$. (Do tego służy breaksopcja sim.) Przetwarzanie wyników symulacji może dostarczyć jeszcze więcej szczegółów.)

Zabawne jest to, że tyle można uzyskać z tego, co w rzeczywistości składa się z trzech wierszy kodu: jeden do symulacji próbek iid z określonego rozkładu, jeden do zastosowania do tablicy rozkładów zerowych, a trzeci do zastosowania do tablica alternatywnych dystrybucji. Są to trzy etapy analizy mocy: reszta to tylko podsumowanie i interpretacja wyników.

(Więcej niż komentarz, być może nie pełna odpowiedź)

[I] spodziewałby się, że wraz ze wzrostem wielkości próby maleje prawdopodobieństwo odrzucenia wartości zerowej

Pomijając rozważania stronnicze testy (które nie są rzadkie pod względem dopasowania, więc warto wspomnieć), istnieją trzy sytuacje związane ze współczynnikiem odrzucenia, które warto rozważyć:

1) współczynnik odrzucenia podczas symulacji od zera (jak wydaje się, że robisz w swoim pytaniu)

Tutaj współczynnik odrzucania powinien być na poziomie istotności lub blisko niego, więc nie, jeśli utrzymasz stały poziom istotności, współczynnik odrzucania nie zmniejszy się, gdy n wzrośnie, ale pozostanie na poziomie bliskim / bliskim$\alpha$.

2) współczynnik odrzucenia przy symulacji z jakiejś alternatywy

Tutaj współczynnik odrzucania powinien rosnąć wraz ze wzrostem n .

3) wskaźnik odrzucenia dla niektórych zbiorów rzeczywistych danych

Praktycznie wartość zerowa nigdy nie jest prawdziwa, a rzeczywiste dane będą zawierały pewną mieszankę ilości nienormalności (mierzonych przez statystyki testowe). Jeżeli stopień nienormalności nie jest związany z wielkością próby, wskaźnik odrzucenia powinien rosnąć wraz ze wzrostem n .

Tak więc w rzeczywistości w żadnej z tych sytuacji nie powinniśmy obserwować spadku wskaźnika odrzucania wraz z wielkością próby.