niedziela, 26 listopada 2017

Parametr μ jako stała (liczba) lub jako zmienna (losowa).

Zabiegi statystyczne służą do tego, aby móc na podstawie dobrze dobranej próby wnioskować o parametrach w populacji. Parametrach - czyli co? Zwykle uczymy się, że parametr to pewna stała liczba, której nie znamy i poznać dokładnie nie możemy, bo całej populacji nie przebadamy. A, i że oznacza się je greckimi  literami, np μ. Z tym drugim nie ma co dyskutować, skoro chcemy upamiętnić gracki alfabet, to niech będzie. Lepszy grecki alfabet niż chińskie znaczki. Natomiast myślenie o parametrze w kategoriach populacji to scheda po Ronaldzie A. Fisherze - genialnym matematyku i statystyku XX w., niemniej jednak był to tylko człowiek, który miał swoje przyzwyczajenia myślowe - postrzeganie parametru jako charakterystyce populacji jest jedną z nich. Z jednej strony można właśnie tak myśleć o parametrze: jako czymś, co charakteryzuje populację, np. średni poziom wzrostu wszystkich Polaków. Z drugiej strony łatwiej ogarnąć dalsze techniki statystyczne, gdy myśli się o parametrze jako o czymś, co bardziej charakteryzuje jakiś rozkład cechy niż samą populację. Dlaczego? Przecież wprawdzie badamy populację, ale dokładniej to celem badań jest rozkład cechy. To rozkład ma parametr (np. wykładniczy) lub parametry (normalny: średnia i wariancja), a nie sama populacja - populacja ma wiele cech, wzrost, przeżywalność, poziom depresji. Przy założeniu, że wzrost ma rozkład normalny to jest sens pytać o jego parametr μ (i średnią z próby wykorzystać jako przybliżenie tego μ). Niby jest to niewielka sprawa, ale w pierwszym przypadku parametr ‘fruwa’ w powietrzu jako charakterystyka populacji, a w drugim jest przytwierdzony do rozkładu cechy.

Tak czy siak, w szkole czy na studiach uczymy się, że parametr ma być liczbą. Wbija się to do głów tak mocno, że pomyślenie o parametrze w populacji w innych kategoriach jest trudne. Dobrze. To może wobec tego wyobraźmy sobie, że parametr nie jest jedną konkretną liczbą, tylko zmienia się, np. w czasie. Ludzie w XXI wieku są średnio wyżsi niż ludzie w XI w. Nie wiem, do końca czy to prawda, bo nie ryję w ziemi w poszukiwaniu średniowiecznych szkieletów, ale coś tam słyszałam, że byli krótsi. W każdym razie chodzi o zmienność wartości parametru μ w czasie. Dziś jest tak, kiedyś było inaczej. Skoro tak, to w takim wypadku μ będzie zmienną losową. Niekoniecznie czas jest tym czynnikiem modyfikującym wartość μ. Bywa, że przedział liczbowy (lub zbiór) jest lepszym opisem parametru. Trzeba byłoby tylko podać prawdopodobieństwa kolejnych wartości μ.
W ten sposób płynnie przeszłam od rozumienia parametru μ jako stałej liczby do pojmowania parametru μ jako zmiennej losowej.
Można? Można. Skąd wziąć taki rozkład? O tym później. Gdzie przydaje takie pojęcie parametru w kategoriach zmiennej losowej, a więc zbioru liczb z odpowiadającymi im szansami pojawienia się? ... W statystyce bayesowskiej (i paru innych miejscach też).

Tym krótkim i cichym postem, bez szumnego tytułu, zapowiadam wejście mojego bloga w kosmiczny pas postów o statystyce bayesowskiej. Robi się ona modna od lat jakoś 90-tych w naukach społecznych, choć pierwszy artykuł w psychologii o statystyce bayesowskiej jest z 1963 roku. Niniejszy post jest cichy, bo cichość ta oddaje mój brak ekscytacji metodami bayesowskimi: jak każda szkoła myśli statystycznej, bayesianizm ma swoje założenia i ograniczenia, plusy i minusy, które przystępnie opiszę w kolejnych postach.
A zatem - wstęga przecięta! :)

wtorek, 24 października 2017

Moc testu

Dokonując wnioskowania statystycznego, badacz może pomylić się na dwa uzupełniające się sposoby: (1) może stwierdzić, że efekt eksperymentalny (leczenia, manipulacji, ogólnie: związek między zmiennymi) istnieje, kiedy w rzeczywistości go nie ma, np. istnieje związek domknięcia poznawczego z samooceną; (2) badacz może nie wykryć efektu eksperymentalnego (leczenia, manipulacji, ogólnie: związku między cechami psychologicznymi), podczas gdy w populacji ten efekt lub związek istnieje. Przykładowo: może uznać, że nie ma związku palenia z zapadalnością na nowotwór płuc (to na szczęście inne wieloletnie badania dowodzą). Te dwa błędy są w zasadzie komplementarne i mają swoje odpowiednie nazwy. Podsumowując, jeśli mamy dwie hipotezy: zerową H0 i alternatywną H1, to decyzje odnośnie ich statusu
logicznego (‘prawdziwa’, ‘fałszywa’) mogą prowadzić do dwóch rodzajów błędu we wnioskowaniu badacza. Badacz może stwierdzić, że hipoteza zerowa H 0 jest fałszywa (a w rzeczywistości jest ona prawdziwa), podjąć decyzję o jej odrzuceniu na korzyść hipotezy alternatywnej H1, która to jest nieprawdziwa. Badacz może również uznać, że hipoteza zerowa H0 jest prawdziwa, podczas gdy w rzeczywistości jest ona fałszywa (za to prawdziwą jest hipoteza alternatywna H1).
Trzeba odróżnić dwie rzeczy: decyzję badacza odnośnie rzeczywistości i stan rzeczywisty rzeczy. Czasami wydaje się nam, że istnieje coś, czego nie ma i nazywamy to np. złudzeniem wzrokowym (sylwetka oprawcy w ciemnościach pokoju lub urojenia paranoiczne). Z kolei przy pomijaniu stronnym obserwuje się, że pacjent nie dostrzega jednej ze stron, mimo tego, że ona przecież istnieje.

Błąd I-go rodzaju, zapisywany grecką literą α, to błąd polegający na odrzuceniu prawdziwej hipotezy zerowej, uznaniu jej za fałszywą i w konsekwencji przyjęciu hipotezy alternatywnej (nieprawdziwej).
Błąd II-go rodzaju, zapisywany jako β, to błąd, który polega na przyjęciu fałszywej hipotezy zerowej, co powoduje zignorowanie hipotezy alternatywej, kiedy to właśnie ona jest w rzeczywistości tą prawdziwą (np. ‘’nie ma związku między paleniem papierosów a ryzykiem nowotworu płuc”).
Tab. 1. Błędy I i II rodzaju (LJK).


Jeśli hipoteza zerowa jest hipotezą o braku różnic lub związku, to badacz popełniający błąd II rodzaju stwierdza, że pomiędzy dwoma cechami psychologicznymi nie ma współzależności, a tak naprawdę one istnieją. Dobrze to ilustruje przykład medyczny. Ustalmy tutaj takie dwie hipotezy:

hipoteza zerowa: brak choroby versus hipoteza alternatywna: obecność choroby.

Badacz popełniający błąd I-go rodzaju stwiedziłby, że pacjent jest chory i zatrzymałby go w szpitalu. Błąd II-go rodzaju polegałby na stwierdzeniu, że chory pacjent jest zdrowy i pozbawienie pacjenta opieki.

Kolizja oznaczeń. Często słyszę, że ‘’poziom istotności to błąd pierwszego rodzaju’‘. Jest to całkowite poplątanie z pomieszaniem, choć przyczynili się do tego sami autorzy tych pojęć (R.A. Fisher oraz J.S. Neyman). α co innego oznacza w podejściu Fishera, a co innego w podejściu Neymana. Rozpatrywanie błędów I, α, oraz II rodzaju, β, jest elementem szkoły Neymana-Pearsona, w której statystyk zestawiał dwie konkurujące hipotezy i przeprowadzał badanie mające na celu wybrać którąś z nich. Natomiast, w szkole Fishera była tylko jedna hipoteza: hipoteza zerowa, dlatego on nie rozważał błędów II-go rodzaju (i jak się później okaże: mocy testu). Dla Fishera, α to poziom istotności (significance level). W książce pisze on tak: ’zwyczajowo i z wygodną dla badacza dobrze byłoby przyjąć 5 procent jako poziom istotności. W tym sensie, że badacz ignorowałby wszystkie te wyniki, które nie przekroczą tego progu (Design of Experiments, s.13)’. Dalej zachęca, aby ignorować wyniki mające prawdopodobieństwo wyższe niż 1 na 20 (s.13). Neyman również zaadoptował oznaczenie α oraz regułę 1-na-20, ale interpretacja tej α w szkole Neymana-Pearsona jest zupełnie inna niż w szkole Fishera. Tutaj, eksperyment (czy ogólnie: badanie) jest teoretycznie wykonywany nieskończenie wiele razy (mimo tego, że badacz raz robi eksperyment). Za każdym razem podejmowana jest decyzja dotycząca prawdziwości hipotezy zerowej bądź alternatywnej.
α i β są odsetkami błędnych decyzji. Ponieważ to są dwie różne szkoły statystyczne, uznanie, że poziom istotności α jest tożsamy z błędem I-go rodzaju α jest samo w sobie błędem. Ten sam znaczek α, a dwie różne interpretacje.

Poziom błędu I-go rodzaju. Zwyczajowo przyjmuje się, że poziom błędu I-go rodzaju wynosi 0.05 (lub: 5%). Jest to zwyczaj, który nie ma swojego naukowego dowodu, zresztą tak jak każda reguła kciuka. Można spotkać badania, gdzie poziom błędu I-go rodzaju ma być niższy i wynosić 0.01%. Możliwe jest również podwyższenie tego błędu do 10% - to jest kwestia oceny badacza w danej dziedzinie nauki. Test, który ma tendencję do nieodrzucania hipotezy zerowej i w związku z tym, jego rzeczywisty poziom błędu I-go rodzaju nie przewyższa ustalonego przez badacza (nominalnego, np. 5%) jest testem konserwatywnym. Żeby popełniać mało pomyłek polegających na odrzuceniu prawdziwej hipotezy zerowej, trzeba mieć w ogóle niechęć do odrzucacania jakichkolwiek hipotez. Ponieważ hipotezy zerowe mówią zwykle o istnieniu jakiegoś efektu, a konserwatywny jest słowem o etymologii łacińskiej (conservare znaczy przechowywać, zachowywać) i ma znaczenie zachowawczy (‘lepiej nic nie zmieniać’), więc test, który ma tendencję do nieodrzucania hipotezy o zerowym efekcie, jest testem konserwatywnym. Test, który ma rzeczywisty poziom alfa wyższy niż 5%, czyli mimo tego, że teoretycznie powinien popełniać nie więcej niż 5% pomyłek w rzeczywistości popełnia ich więcej, bo częściej odrzuca hipotezę zerową na korzyść alternatywnej, zwany jest testem liberalnym (antykonserwatywnym).

Poziom błędu II rodzaju. Zwykle żąda się, aby błąd II-go rodzaju w teście statystycznym pozostawał na poziomie 20%. Oznacza to pomyłkowe orzeczenie zdrowia u co piątego pacjenta (nawiązując do medycznej ilustracji pojęcia). Rysunek 1 przedstawia rozkłady statystyk testowych (nie: cech 1 ). Czarnym kolorem jest zaznaczony rozkład statystyki testowej gdy hipoteza zerowa H 0 jest prawdziwa. Zielonym kolorem jest zaznaczony rozkład zgodny z hipotezą alternatywną H1. Różowy punkt na dole to wartość krytycznej statystyki testowej - wartości statystyki testowej większe od niej spowodują odrzucenie hipotezy zerowej. Chcielibyśmy, aby w sumie ich prawdopodobieństwa nie przekroczyły pewnego pułapu, czyli nie przekroczyły poziomu błędu I-go rodzaju. Stąd też, niebieskie pole pod czarną krzywą na prawo od różowego punktu odpowiada błędowi I rodzaju, α i wedle naszego życzenia, chcielibyśmy aby jego powierzchnia wynosiła α = 0.05.

Rys. 1. Poziom błędu I i II  rodzaju (LJK).
Zielona krzywa na prawo od różowego punktu opisuje prawdopodobieństwa wartości tych statystyk testowych, które są mniejsze od niego. Tutaj też chcielibyśmy, aby sumaryczne prawdopodobieństwo nie przekraczało progu. Próg ten odpowiada prawdopodobieństwo popełnienia błędu II-go rodzaju, β.
Zgodnie z regułą kciuka, chcielibyśmy, aby czerwone pole wynosiło 20%.

Rzeczywisty a nominalny poziom błędu II-go rodzaju. Sprawdźmy symulacyjnie za pomocą programu R, jaki jest rzeczywisty poziom błędu II rodzaju dla testu t-Studenta, porównującego dwie grupy (klik), który ma zbadać takie hipotezę:
H0 : μ1 − μ2 = 0
przeciwko jawnie określonej hipotezie alternatywnej
H1 : μ1 − μ2 = 0.1

Opis symulacji. Badanie wykonujemy w następujący sposób: generujemy dwie próby z rozkładu normalnego o ustalonej liczebności. U nas będzie to 30 w każdej z grup. Pierwsza z nich będzie miała średnią μ1 równą zero i wariancją σ2 = 1 (zatem cecha w tej grupie ma rozkład normalny, standardowy), druga będzie miała średnią μ2 równą 0.1 i wariancją σ2 = 1. Stąd na pewno wiemy, że hipoteza zerowa jest nieprawdziwa a hipoteza alternatywna - prawdziwa (0 = 0.1). Sprawdzamy, czy p-wartość będzie mniejsza od 0.05, czyli czy testowi uda się odrzucić (błędną) hipotezę zerową. Jeśli tak, to test nie popełnia błędu II-go rodzaju, jeśli nie - test popełnia błąd II-go rodzaju.
Tutaj trudno jest badać rzeczywisty poziom błędu I-go rodzaju, ponieważ z góry ustaliliśmy, że hipoteza zerowa jest fałszywa (a błąd I-go bada sytuację, gdy hipoteza zerowa jest prawdziwa). Całą tę procedurę - od wygenerowania nowych prób do sprawdzenia decyzji testu - powtarzamy co najmniej kilkaset razy (np. tysiąc razy). Ręcznie zajęłoby to wieki, na szczęście są od tego komputery i ich moce obliczeniowe.
set.seed(20062017)
errII<-sum(replicate(1000,t.test(rnorm(30,0,1), rnorm(30,0.1,1))$p.value)>0.05)
paste("błąd II-go rodzaju wynosi ", errII/1000, sep="")
Błąd II-go rodzaju wynosi 93,5%. Jest to wartość o wiele większa niż zakładane 20%. Zatem dla testu  t-Studenta dwóch 30-elementowych grup niezależnych, gdzie cecha (w jednej i w drugiej grupie) ma rozkład normalny o tej samej wariancji σ1 = σ2 = 1 z tą różnicą międzygrupową, że średnia w pierwszej grupie μ 1 wynosi 0, a średnia w drugiej grupie μ2 wynosi 0.1, błąd II-go rodzaju wynosi około 93%. Oznacza to, że jeśli różnice w średnich między dwoma grupami są takie małe, to dla prób o liczności 30 elementów, test t-Studenta ma problem z wykryciem tej różnicy (będziemy później o tym mówić, że moc tego testu w tych specyficznych warunkach jest słaba, niska).

Moc testu statystycznego. Dzięki symulacjom dowiedzieliśmy się, że rzeczywisty błąd II-go rodzaju dla testu t-Studenta, gdzie obie grupy mają po 30 obserwacji a średnie w dwóch grupach różnią o 0.1, wynosi 93.5%. Oznacza to,
że 6% przypadków test jest zdolny prawidłowo odrzucić fałszywą hipotezę zerową i przyjąć hipotezę alternatywną, ponieważ to ona jest prawdziwa. Zazwyczaj mówi się wtedy, że moc takiego testu wynosi 6% (ponieważ 100% − 93.6% ≈ 6%).
Do oceny zdolności testu do nie-popełnienia błędu II-go rodzaju posługujemy się pojęciem mocy testu statystycznego. Oblicza się z prostego wzoru:
moc = 1 − β gdzie β to błąd II-go rodzaju.
Moc nie jest liczbą odwrotną do β (liczbą odwrotną do liczby a nazywamy liczbę a 1 ), ani przeciwną (liczbą przeciwną do liczby a nazywamy liczbę −a). Po prostu nie ma nazwy na związek liczb a i b,które sumują się do jedynki, a + b = 1.
Moc testu to prawdopodobieństwo odrzucenia fałszywej hipotezy zerowej, gdyby eksperyment wykonywano wiele razy (long-run probability, w interpretacji częstościowej).

Ogólnie można rozumieć moc testu statystycznego jako zdolność tego testu do prawidłowego odrzucenia błędnej hipotezy zerowej i przyjęcia prawdziwej hipotezy alternatywnej. Oczekiwanie tego, aby błąd II-go rodzaju wynosił 20%, pociąga za sobą moc testu na poziomie 80% (ponieważ 1 − 20% = 80%).

Poniższy rysunek przedstawia relacje między błędem I - go rodzaju (niebieskie pole), II-go rodzaju (czerwone pole) a mocą (pomarańczowe pole). Skoro moc to jest jeden minus błąd drugiego rodzaju, to czerwone pole i pomarańczowe pole powinny składać się na pole pod krzywą po prawej stronie.

Rys. 2. Relacje między błędem I, II rodzaju oraz mocą (LJK).


O mocy testu statystycznego można myśleć tylko i wyłącznie w warunkach gdy jest jasno określona hipoteza alternatywna - bo to na jej podstawie jest określany błąd II-go rodzaju. Bez jasno zdefiniowanej postaci hipotezy alternatywnej (np. μ 1 − μ 2 = 0.5) nie ma sensu mówić o błędzie II-go rodzaju, ani o mocy testu statystycznego. Powyższe przykłady są dość specyficzne, bo bezpośrednio podają postać hipotezy alternatywnej, ale badacze zwykle radzą sobie za pomocą wielkości efektu (effect size), która to opisuje wielkość zjawiska w populacji, np. d Cohena albo r Pearsona. Przykłady zostały tak dobrane, że
wielkość efektu jest równa hipotezie alternatywnej.

Ile potrzeba obserwacji, aby ten błąd spadł do 20%? Weźmy liczebność obserwacji n 1 = n 2 = 1571 w obu grupach i dla takich sprawdźmy błąd II-go rodzaju (a zarazem moc testu).
set.seed(20062017)
errII<-sum(replicate(1000, t.test(rnorm(1571,0,1), rnorm(1571,0.1,1))$p.value)>0.05)
print(paste("Błąd II-go rodzaju, beta, wynosi",errII/1000), sep="")
## [1] "Błąd II-go rodzaju, beta, wynosi 0.178"
print(paste("Moc wynosi 1-beta =",1-errII/1000), sep="")
## [1] "Moc wynosi 1-beta = 0.822"

Wychodzi na to, że dopiero, gdy badacz przebada 2 · 1571 osób badanych, może być spokojny o to, że moc testu statystycznego wyniesie 80% a test będzie zdolny wykryć różnicę między średnimi w populacji równą 0.1.

Dobrze, ale skąd wzięłam akurat tę liczbę: 1571? Czy są jakieś sposoby wyznaczania ilości osób badanych? Odpowiedź brzmi: tak - dzięki wielkości efektu. Ustalając przed rozpoczęciem zbierania danych, moc wybranego testu statystycznego (β = 80%), poziom błędu I-go rodzaju (α = 0.05) oraz wielkość efektu, da się zaplanować liczbę obserwacji. Akurat w moim przykładzie, gdzie wariancje w obu grupach były identyczne i równe 1, wartość wielkość efektu dla dwóch grup, d-Cohena, uwzględniająca wariancję w obu grupach, jest tożsama z formułą zawartą w hipotezie alternatywnej μ1 − μ2 = 0.1. Zwykle jednak
wygodniej jest odwoływać się do wielkości efektu niż bezpośrednio do postaci hipotezy alternatywnej. Aby test był w stanie rzetelnie wykryć mały efekt, należy przebadać 2 · 1571 osób (czyli łącznie 3142.). Więcej o wielkości efektu Klik

Moc, błąd I-go rodzaju, wielkość efektu i liczebność próby są ze sobą związane w taki sposób, że wystarczy znać wartości dowolnych trzech z nich, aby wyliczyć wartość czwartej. Jeśli wiemy, że a = 1, b = 2, c = 3 oraz że a + b + c + d = 10 łatwo wyliczyć wartość d, podstawiając w tym równaniu wartości poszczególnych danych 1 + 2 + 3 + d = 10. Stąd d = 10 − 6 i d = 4.

Analiza mocy post-hoc.
Obliczanie mocy po wykonaniu badań (post-hoc). Analiza mocy przydaje się do dwóch rzeczy: planowanie wielkości próby oraz do meta-analiz (które nie
są tu omawiane). Można też zastanawiać się jaką moc będą miały testy przy oczekiwanej różnicy między zmiennymi, związku między nimi (ogólnie: przy znanej lub spodziewanej wielkości efektu). Jeśli badamy moc przed zebraniem danych, to wówczas taka analiza mocy jest prospektywną analizą mocy. Wówczas ma ona sens - badacz dowiaduje się, jakiej wielkości próby potrzebuje lub jaką moc będą miały jego testy. Część badaczy natomiast oblicza moc testu po wykonaniu analiz, mając dane: ilość obserwacji, które udało się zebrać oraz otrzymaną wielkość efektu (w próbie). W ten sposób obliczają moc testu post
hoc, czyli wykonują retrospektywną analizę mocy. Interpretacja takich wyników nastręcza trudności, co w konsekwencji prowadzi do pytania o sens retroanalizy. Dlaczego? Jeśli do wyznaczenia poziomu mocy korzystamy z wielkości efektu uzyskanej z próby, to zakładamy, że wielkość efektu w próbie jest tą
7wielkością efektu, która występuje w populacji, co samo w sobie jest kontrowersyjne, ponieważ prowadzi do udowadniania czegoś, co chcemy udowodnić. Niech posłuży ten przykład. Badacz chciałby dowiedzieć się jakiej mocy będzie jego test, jeśli pobierze dwie próby o tej samej liczebności 30 oraz będzie spodziewał się różnicy między grupami równej 0.1
(wielkość efektu będzie tu również równa 0.5).
set.seed(5072017)
proba1<-rnorm(30,0,1) #dla 1072017 są takie same
proba2<-rnorm(30, 0.5,1)
power.t.test(n=30,delta=0.5, sd=1) #moc wynosi 0.48
Okazuje się, że moc jego testu wynosi zaledwie 48%.
Drugi badacz, który również dokonuje porównań międzygrupowych najpierw przeprowadza badanie, potem wykonuje analizę testem t-Studenta, bada różnicę między średnimi w grupach a następnie oblicza moc testu.
t.test(proba1, proba2)
(ds=abs(mean(proba1)-mean(proba2))) #różnica między średnimi w próbach
power.t.test(n=30,delta=ds) #moc = 0.8
Jego moc testu jest wysoka, wynosi 80%.
Z tego punktu widzenia, trochę bezsensowna jest analiza mocy post-hoc na podstawie wielkości efektu z próby, ale to jest to, co robi obecny SPSS.

Zależność między błędami.
Często mówi się o tym, że błędy α i β są ze sobą związane - zmieniając jeden, zmieniamy też i drugi. Zależność tę można wykazać analitycznie (tj. za pomocą formuł matematycznych), można też zilustrować dwoma skrajnymi przykładami: wówczas gdy badacz minimalizuje jeden z nich. Poniższe dwa przykłady pokazują działanie testu jednostronnego, aby uniknąć zaciemnienia rysunków.
Konsekwencje minimalizacji błędu I-go rodzaju alfa = 0. Jeśli chcemy, aby poziom błędu I rodzaju był równy 0, to automatycznie maksymalizujemy błąd II-go rodzaju. Jak to należy rozumieć? Aby nie popełniać błędu I rodzaju, czyli uniknąć odrzucenia prawdziwej hipotezy zerowej, trzeba na wszelki wypadek za każdym razem przyjmować, że hipoteza zerowa jest prawdziwa. Wprawdzie będziemy się mylić i uznawać brak efektu za prawdę, ale przynajmniej nie pomylimy się i nie stwierdzimy efektu w sytuacji, gdzie go nie ma.

Przykład prawniczy. Mamy 100 podejrzanych o popełnienie przestępstwa. Hipoteza zerowa brzmi: osoba jest niewinna (brak winy). 
Hipoteza alternatywna: osoba jest winna (wina). 
Błąd I-go rodzaju oznaczałby uznanie osoby niewinnej za przestępcę (błędne odrzucenie hipotezy zerowej o niewinności). Aby nie popełnić błędu I-go rodzaju w 100 przypadkach należy uznać, że wszyscy są niewinni. Przyjmując hipotezę zerową jako prawdziwą we wszystkich tych przypadkach, wyzerujemy błąd I-go rodzaju. Osoby naprawdę niewinne nie otrzymują kary. Natomiast wypuszczamy na wolność również przestępców. Błąd I-go rodzaju jest zminimalizowany (do zera), ale za to błąd II-go rodzaju ma najwyższą wartość (jest równy 1).  

Konsekwencje minimalizacji błędu II-go rodzaju beta = 0. W przykładzie prawniczym błąd II-go rodzaju oznacza błędne odrzucenie prawdziwej hipotezy alternatywnej (dana osoba jest przestępcą) i przyjęcie fałszywej hipotezy zerowej (dana osoba jest niewinna). Chcemy w żadnym ze 100 przypadków nie popełnić tego błędu. Zatem uznajemy, że wszyscy są winni. Przestępcy słusznie siedzą za kratkami, ale osoby niewinne również wędrują do więzienia. Błąd II-go rodzaju jest równy zero, ale zmaksymalizowaliśmy błąd I-go rodzaju. Chcielibyśmy zminimalizować błąd I-go i błąd II-go rodzaju. Jednocześnie zdajemy sobie sprawę z tego, że wyzerowanie go spowoduje zmaksymalizowanie błędu I-go rodzaju (patrz sekcja 1.2.3). Jaki byłby odpowiedni kompromis między poziom błędu I-go a II-go rodzaju? Propozycja 20% dla błędu II-go rodzaju wydaje się być sensowna. Oczekuje się, dany test dla ustalonej hipotezy alternatywnej popełnia błąd II-gow 1 z 5 przypadków.
Interpretacja moc = 0.5 Jak zinterpretować wartość mocy równą 0.5? Jeśli wiemy, że moc = 0.5, to zakładając, że hipoteza alternatywna jest prawdziwa, błąd II-go rodzaju wynosi β = 1 − moc = 1 − 0.5 = 0.5. Oznacza to, że w połowie przypadków nie stwierdzamy obecności związku (różnicy). Nie, nie jest to dobra wiadomość.

Ile osób mam przebadać? To jest pytanie, na które każdy badacz chciałby znać odpowiedź, bo nawet na chłopski rozum wychodzi, że im więcej, tym lepiej, ale to myślenie ma swoje ograniczenia, np. niektóre zjawiska są bardzo, bardzo rzadkie albo kosztowne.
To zależy od szkoły: Fisher powiada, że aby uzyskać istotny statystycznie wynik, im więcej, tym lepiej. Dla Neymana: skorzystaj z wielkości efektu (najpierw zaplanuj ją!), w podejściu bayesowskim wielkość próby nie jest aż tak ważna - 'dobadanie' próby powoduje takiego problemu jak z p-wartością.

niedziela, 8 października 2017

Centralne Twierdzenie Graniczne CTG

Centralne Twierdzenia Graniczne to grupa praw mówiących o zachowaniu się sum zmiennych losowych.

Najbardziej rozpowszechnione w podręcznikach do statystyki jest poniższe:


Niech X1, X2, … Xn będą niezależnymi zmiennymi losowymi o tym samym rozkładzie, każda o tej samej wartości oczekiwanej EX = 0 oraz wariancji σ2 = 1. Wtedy:
Wtedy rozkład średnich tych zmiennych losowych jest coraz bliższa rozkładowi normalnemu.

I po krzyku. Pozostałe twierdzenia to przypadki, kiedy średnia każdej z tych zmiennych jest różna od zera, a wariancja jest różna od jedności, a także kiedy zmienne nie mają tego samego rozkładu.
Wyjaśnienia wymaga punkty: rozkład, rozkład cechy a rozkład średnich cechy oraz skoro rozkład średnich to może jest jeszcze rozkład wariancji? I o tym jest poniższy wpis.


Rozkład.
Pojęcie rozkładu zostało już wprowadzone tu (klik). Przypomnę, że w dużym uproszczeniu rozkład zmiennej można utożsamiać z częstotliwością występowania poszczególnych wyników danej cechy.

Rozkład średnich a rozkład cechy.
Przypuśćmy, że wybrana zmienna psychologiczna ma rozkład beta (klik). Już może mniejsza o to, co można robić z rozkładem beta, poza faktem, że jest to rozkład dobry dla tych cech, które mogą przyjmować wartości od zera do jednego. Specjalnie wybrałam taki dziwny rozkład. Próba złożona z 30 osób mogłaby mieć taki rozkład jak po lewej stronie rysunku. Pionową, czerwoną kreską jest zaznaczona średnia. Widać, że rozkład jest dwugarbny (dwumodalny).
Rozkład cechy (po lewej) a rozkład średnich (po prawej). Ryc:LJK.


Teraz wyobraźmy sobie, że robimy wykres średnich. Mamy już nie jedną trzydziestoelementową próbę z takiego beta rozkładu, ale tysiąc. Wtedy wykres wygląda tak, jak na rysunku po prawej stronie. Jest jednogarbny (jednomodalny), w miarę symetryczny. Zwiększając liczbę obserwacji w każdej z tych tysiąca prób wykres średnich coraz bardziej przypominałby rozkład normalny.
Taka jest różnica między rozkładem cechy (rysunek po lewej), który może być bardzo daleko od normalnego, a rozkładem średnich cechy, który coraz bardziej jest symetryczny i gaussowski.

Rozkład wariancji.
Tak, popuśćmy wodze fantazji kawalerskiej, aby poczuć różnicę między rozkładem cechy a rozkładem średnich. Wymyślmy sobie, że chcemy zbadać rozkłady wariancji. I znowu: potrzeba mi więcej niż jedną próbę, aby móc wyrysować wykres i zauważyć jakiś kształt. Pozwolisz, drogi Czytelniku, że wykorzystam powyższy tysiąc prób z rozkładu beta - na pewno mają jakieś wariancje.

Rozkład zestawu tysiąca wariancji w n=30-elementowych próbach. Ryc:LJK.



W tym miejscu pewnie już  Czytelnik odczuwa trzewiami różnicę między rozkładem cechy, a rozkładem (pewnej ilości) średnich cechy, stąd dobrze byłoby wspomnieć o błędzie w niektórych podręcznikach.

Ważne.
Ponieważ CTG opisuje zachowanie się średnich a nie cechy, to niekiedy w podręcznikach można znaleźć taki błąd, pewnie wynikający ze skrótowości, że: 'CTG oznacza, że zwiększając liczebność próby zwiększamy prawdopodobieństwo uzyskania rozkładu normalnego danych'. Otóż, jest to nieprawda. Zwiększając liczebność próby zwiększamy prawdopodobieństwo uzyskania rozkładu normalnego średnich.
Byłoby to bardzo dziwne, gdyby gromadzenie danych zmieniałoby jego naturę.























piątek, 7 października 2016

Bootstrap i bootstrapowanie średniej.

Spróbuję upiec dwie pieczenie na jednym ogniu – opowiedzieć jak działa bootstrap* i pokazać, co to znaczy: rozkład średnich wybranej cechy w populacji. Dlaczego to ostatnie? Otóż niekiedy spotyka się założenie o normalności cechy, np. w teście t-Studenta lub w ANOVA, lub w analizie regresji i bywa, że jest to zbyt mocne założenie. Zależy nam, aby owszem rozkład był normalny, ale nie wyjściowej cechy tylko ... pewnego przekształcenia wartości tej zmiennej. Takim przekształceniem jest średnia. Ale co to znaczy rozkład średnich w próbie? Dlaczego średnia miałaby mieć jakiś rozkład? O tym będzie niżej. A co ma do tego bootstrap? O tym będzie jeszcze niżej...

Spróbujmy zrozumieć na konkretnym przykładzie.
Mamy próbę, X jest zmienną psychologiczną – wyniki dwudziestu osób w skali ABC, która przyjmuje wartości od 1 do 10. Te dwadzieścia osób otrzymało taki oto wynik: 

0 2 1 7 6 3 5 1 10 3 2 6 5 10 5 1 9 8 7 5 

Przy okazji moglibyśmy dowiedzieć się coś o tym, ile jest wystąpień konkretnych wartości w próbie, np. ile jest osób o wyniku trzy, ile siedem, itd. Rozłożenie występowania wyników skali ABC przedstawia poniższy histogram.



Rys. Rozkład wartości skali ABC dla dwudziestu osób (LJK).
Rozkład ten można spróbować zinterpretować. Najczęściej występującą wartością jest 5 (aż cztery razy). Średnia w tej próbie wynosi 4,80 (użyłam wzoru na średnią arytmetyczną). Średnia w próbie (ta obliczana na podstawie wzoru) jest estymatorem średniej w populacji (czyli średniego poziomu cechy). Skoro jest estymatorem, to tak jak każdy estymator jest zmienną losową. Gdybyśmy przebadali inną grupę dwudziestu osób, to średnia mogłaby wyjść nie 4,80, a np. 4,52. Widać tutaj, że wartości średniej w próbie mogą zmieniać się w zależności od przebadanej próby. 
Zatem chcielibyśmy dowiedzieć się, w jakich granicach można spodziewać się wartości średniej w próbie. Innymi słowami, pytamy jak precyzyjny jest pomiar tej średniej. Zwykle korzysta się z przedziałów ufności. Interpretacja przedziałów ufności jest podana tutaj (klik). W przypadku naszej próby przedział ufności dla średniej wynosi [3.35,6.25]. Obliczyłam go na podstawie wzorów.
 
Centralne Twierdzenie Graniczne (CTG). 
Problem z przedziałami ufności obliczanymi w przy pomocy wzorów jest taki, że podskórnie zakładają one korzystanie z Centralnych Twierdzeń Granicznych. 
Centralne Twierdzenia Graniczne to taka grupa twierdzeń, które określają warunki, dzięki którym grupy rozkładów mogą zbliżać się do pewnego rozkładu wzorcowego. Na przykład coś takiego: jeśli będziesz mieć odpowiednio dużą liczbę elementów w próbie (czyt. przebadanych osób), to rozkład średniej w próbie będzie coraz bliżej rozkładowi normalnemu. 

Rozkład średnich w próbie. 
Jak to działa? Przykład. Mam 1000 grup po 20 osób każda. Każda z tych osób została zbadana pod kątem interesującej psychologa cechy. Jaki typ rozkładu mają średnie tej cechy? Jeśli ta cecha ma w populacji rozkład normalny, to nie dość, że wykres będzie wykresem rozkładu normalnego (to trochę stylistyczne powtórzenie), to rozkład średnich również jest rozkładem normalnym.


Rys. Rozkład cechy w populacji i rozkład średnich w próbie (LJK).


Dla cechy która ma rozkład normalny, rozkład średnich w różnych próbach będzie rozkładem normalnym (na to jest odpowiedni dowód matematyczny). Natomiast dla cechy, której rozkład nie jest rozkładem normalnym, korzysta się z CTG. 
Poniżej mamy cechę z rozkładu wykładniczego. Jej rozkład jest po lewej stronie. A po prawej stronie mamy rozkład średnich tej cechy w tysiącu dwudziestoelementowych próbek.
 
Rys. Rozkład cechy w populacji i rozkład średniej cechy w próbie (LJK).
Rozkład średniej w próbie jest nieco skośny z prawej strony. Rozkład normalny tymczasem jest symetryczny. 
Skoro w obliczaniu końców przedziału ufności korzysta się z kwantyli jakiegoś rozkładu, to tym sposobem zaocznie wskazuje się na ten rozkład. My używamy kwantyli rozkładu normalnego. Może ktoś pamięta wartość 1.96? 1.96 jest kwantylem rozkładu normalnego. 
Więc... spróbujmy postąpić inaczej. Gdybyśmy skądś wiedzieli, jaki jest rozkład wszystkich dwudziestoelementowych średnich tejże cechy, to łatwo byłoby skonstruować taki przedział (niżej). 
Ale tego nie mamy, możemy za to próbować odgadnąć jak wygląda taki rozkład. W tym celu przyda się bootstrap. Jak to wygląda? 

Boostrap.
Na bootstrap składają się dwa pojęcia – repróbkowanie za pomocą losowania ze zwracaniem oraz zasada wtyczki (plug in principle) – czyli metaforyczna nazwa na użycie formuły (po ludzku: wzoru) pozwalającej obliczyć interesującą nas statystykę opisową. 

Losowanie ze zwracaniem.
Najpierw poznamy pojęcie losowania ze zwracaniem. Wyjaśnię to w obrazowy sposób. Mamy naszą oryginalną próbę dwudziestu wartości cechy X. Repróbkowanie ogólnie będzie polegało na utworzeniu nowej próby z elementów tej oryginalnej. Re-próbka ma mieć tyle samo elementów, czyli też dwadzieścia. Powstanie za pomocą mechanizmu losowania ze zwracaniem (nie jest to jedyny sposób repróbkowania, ale w bootstrapie robimy właśnie ten).
Losowanie ze zwracaniem polega na tym, że po wylosowaniu elementów z worka i odnotowaniu jego wartości z powrotem wkładamy go do wora. Można go wylosować jeszcze raz, i jeszcze jeden i w ogóle może zdarzyć się tak, że tylko ten jeden wylosujemy te dwadzieścia razy (ale naprawdę rzadko tak się dzieje :).

Przykład losowania ze zwracaniem. Oryginalna próba ma trzy elementy (np. kulki z etykietami): {a,b,c}. Losujemy trzy razy ze zwracaniem, za każdym razem zapisując etykietę elementu i odkładając ją do puli). Próba jaka mogłaby powstać to: {b,b,c}. Inna to mogłaby być {a,b,a}. Jeszcze ktoś inny dostanie {c,c,a}.

Tak utworzona próbę wtórną {b,b,c} nazywa się próbą bootstrapową. Spróbowałam otrzymać jedną próbę bootstrapową dla naszych danych i mam taki wynik: 
3 1 1 1 3 6 2 1 5 1 5 2 5 6 1 3 1 6 3 0 
Co się stało? Bootstrapowa próba różni się od oryginalnej. Mamy siedem razy wartość 1 – w oryginalnej próbie były trzy. Na poniższym rysunku widać rozłożenie poszczególnych wartości.
Rys. Rozkład wartości pojedynczej próby bootstrapowej (LJK).

Czy można policzyć średnią w tej bootstrapowej próbie? Na to pytanie odpowiada zasada wtyczki.
Zasada wtyczki (plug-in principle).
Jest to tak intuicyjne, że w pierwszej chwili można nie wiedzieć o co chodzi, więc po kolei.

My bootstrapujemy teraz średnią. Wzór na średnią (zsumować wszystkie wartości i podzielić przez ich ilość) działa dla jakichkolwiek liczb.

Tak więc można zastosować wzór na średnią do naszej bootstrapowej próby numer 1. W wyniku zastosowania wzoru otrzymamy liczbę: 2.8. Cóż, dużo mniej niż ta oryginalna.
Skoro mamy jedną próbę bootstrapową i jedną bootstrapową wartość średniej, to tę procedurę można wykonywać nieskończenie wiele razy. Okazuje się, że rozsądnie jest zrobić to co najmniej 1000 razy.


Będziemy mieć tysiąc prób boostrapowych i tysiąc obliczonych dla każdej z nich wartości średniej.
Można byłoby zobaczyć, jak wygląda histogram takich średnich.
 

Rys. Histogram wartości średnich bootstrapowych (LJK).

Histogram to nic innego jak rozkład. Pionowa kreskowana linia to bootstrapowa wartość średniej.

Mając rozkład możemy policzyć przedział ufności bez używania wzorów, które to z kolei odwołują się do CTG. Jest kilka sposób obliczania przedziałów, my weźmiemy najprostszy i mamy wynik: [3,45; 6,20]. W porównaniu do przedziału ufności obliczanego tradycyjnie przez wzory jest on trochę węższy (tamten wynosił [3,35; 6,25]). Bootstrapowy przedział ufności tym się różni od tradycyjnego, że nie odwołuje się do Centralnych Twierdzeń Granicznych (CTG). Obliczając go, opieramy się na empirycznym rozkładzie estymatora średniego poziomu cechy w populacji, czyli średniej w próbie.

Tak wygląda najbardziej podstawowe (i popularne też) wykorzystanie bootstrapu i najprostszy bootstrapowy przedział ufności dla średniej. 

Posłowie
*Wiem, że jest propozycja ze strony pewnego uznanego psychologa, aby w języku polskim bootstrap nazywał się „samowspornością”. Cóż, podczas moich studiów na matematyce, ani też późniejszych rozmów z matematykami nikt nie znał tego pojęcia. Bootstrap to bootstrap – nie ma dobrego polskiego odpowiednika, tak jak na słowo 'computer' (choć Chińczycy zrobili z tego „elektryczny mózg”, 电脑). Mnie nie przekonuje propozycja polskiego odpowiednika, po prostu nie brzmi mi dobrze. No, ale jak się przyjmie, to się przyjmie. Uważam, że się nie przyjmie, bo trzeba byłoby przekonać grona innych naukowców nie tylko ze środowiska psychologicznego.