wtorek, 2 stycznia 2018

Dystrybuanta rozkładu prawdopodobieństwa.

Zaczynamy zajęcia ze statystyki, prowadzący lub prowadząca przedstawia podstawy rachunku prawdopodobieństwa, zdarzenie elementarne, losowe i przechodzimy do zagadnienia przedstawienia jakoś szans pojawienia się wyników zmiennej losowej - rzutu monetą, kostką, okiem i beretem. I widzimy taki jakiś wzorek:
F(t) = P( -∞< t) 
i koniec. Język matematyczny to bardzo ładny i zwięzły język, którym można wiele myśli wyrazić i wiele idei przekonać. Wiele nie znaczy: wszystko i dlatego mogę pisać tego bloga. W tym poście przedstawię pojęcie matematyczne - dystrybuantę - które służą do opisu szans wystąpienia możliwych wartości zmiennej losowej. Jakiejś zmiennej losowej.


Prawdopodobieństwo trzeba jakoś przedstawić. Można to zrobić w postaci: (a) funkcji prawdopodobieństwa, (b) gęstości, (c) zbiorczo: w postaci skumulowanej, tj. pokazać, jakie jest prawdopodobieństwo spotkania obserwacji do wybranego poziomu, np. jakie jest odsetek dzieci, który osiągnęły wzrost do danego, od najmniejszego możliwego.
Wychodzi na to, że nie wystarczy podać odsetka dzieci, które osiągnęły dany wzrost, tylko wszystkie mniejsze też. Trochę to nieintuicyjne na pierwszy rzut oka, bo na pytanie ile masz centymetrów wzrostu, odpowiadamy konkretnie: "Mam 172 cm", a nie "Do 172", ale uwierzcie mi, że to ma sens przy korzystaniu ze statystyki. W tym miejscu chodzi o to, aby rozmówcy zadać pytanie: 'jaka jest szansa, że Twój wzrost wynosi do 172 cm'. To na razie załatwia sprawę skumulowania. Ale dystrybuanta to nie synonim na skumulowane prawdaopodobieństwo, choć mają wiele wspólnego.Dalej okaże się, co jest różne.



Aby przekształcić rozkład prawdopodobieństwa (podany w powyższej postaci) do dystrybuanty należy w kolejno dodawać do siebie wartości tego prawdopodobieństwa. Na pierwszym miejscu jest liczba oczek, na drugim miejscu skumulowane prawdopodobieństwo uzyskania liczby oczek do podanej:



W przypadku rozkładów dyskretnych (rozkład dyskretny?klik), takich jak rzut kostką (monetą również) punkty, których następuje dodanie prawdopodobieństwa, nazywamy punktami skoku. Więc punktami skokowymi (skoku) są wszystkie
możliwości wyrzutu oczek.

Rys. Skumulowany rozkład zmiennej o tytule Rzut kostką (LJK).

Na powyższym rysunku wygląda to tak, jakby wykres miał wartości jedynie w czarnych punktach a poza nimi nic nie było. Tak może być tylko i wyłącznie gdy rozmawiamy o skumulowanym prawdopodobieństwie w potocznym tego słowa znaczeniu, gdzie nie ma sensu rozmawiać o tym, ile ono wynosi w punkcie 3.5, bo nie można wyrzucić trzy-i-pół oczka. Okazuje się, że różnica między skumulowanym prawdopodobieństwie a dystrybuantą jest taka, że dystrybuanta muruje dziury w wykresie skumulowanego prawdopodobieństwa: zobaczcie sami (niżej).

Rys. Dystrybuanta rozkładu zmiennej o tytule Rzut kostką (LJK).


Dystrybuanta jako pojęcie matematyczne to nieco szersza sprawa niż skumulowane prawdopodobieństwo. Żąda się bowiem, aby dystrybuanta, oznacza zwykle: F (x), była określona dla wszystkich liczb rzeczywistych, czyli dla wszystkch x ∈ R, nawet jeśli cecha nie pokrywa wszystkich liczb rzeczywistych (np. nie można wyrzucić trzy i pół oczka, nikt nie ma ujemnej samooceny).

Z tej przyczyny powstaje pozorny problem - co z takimi wartościami dystrybuanty, których cecha nie przyjmuje? Nie ma problemu dla wartości wyższych niż maksymalny poziom cechy - dla nich skumulowane prawdopodobieństwo będzie równe 1. Ponadto, nawet jeśli cecha nie przyjmuje jakiejś wartości, to może mieć wartości z jakiegoś mniejszego zakresu, np, P (X < 3.5) = P (X < 3) z tego powodu, że cecha nie przyjmuje
wartości z przedziału (3, 3.5) - ale przyjmuje wartości do 3 - w związku zachowujemy się tak, jakby pytano o liczbę oczek do trzech, zamiast trzy i pół. Wygląda na to, że prawdopodobieństwo przyjęcia wartości z przedziału od ponad 3 do 3.5 jest równe zero;P (3 < X < 3.5) = 0. Stawiamy zerowe prawdopodobieństwo tam, gdzie cecha nie przyjmuje takiej wartości.
Dodatkowo, zgodnie z intuicją skoro żadne prawdopodobieństwo nie jest ujemne, to cokolwiek dodawane do całości kumulacji wywoła powiększenie tej sumy - nigdy jej nie zmniejszy. Innymi słowami, kumulacja można pozostawać na stałym poziomie, ale nigdy nie może się zmniejszać. W ten sposób otrzymaliśmy trzy własności dystrybuanty - zapiszemy je formalnie.

Kryteria dystrybuanty - która funkcja jest dystrybuantą, a która nie jest?
Aby funkcja była dystrybuantą musi spełniać poniższe trzy kryteria:

1. dla x biegnących do minus ∞ F (x) = 0 oraz dla x biegnących do ∞ F (x) = 1, czyli:

2. prawostronnie ciągła
3. niemalejąca - czyli nie ma ujemnych prawdopodobieństw

To są formalne warunki dystrybuanty. Najciekawsze jest to, że jeśli narysujemy dowolną funkcję, lecz spełniającą powyższe kryteria to na pewno opisuje rozkład jakiejś zmiennej. Nie wiadomo, jaka i czy jest odkryta, ale już wiemy, że będzie to rozkład.

Dystrybuanta teoretyczna a empiryczna.
W matematycznym depozycie teoretycznym znajdują się idealne rozkłady - między innymi normalny modelujący, wykładniczy, jednostajny i inne.

Dystrybuanta teoretyczna to dystrybuanta wynikająca z teorii matematycznej, można ją potraktować jako dystrybuantę cechy w (niedostępnej) populacji. Dystrybuanta empiryczna to dystrybuanta uzyskana na podstawie próby.
Niech x będzie dowolnym wynikiem zmiennej (cechy) - możliwym bądź nie. Po prostu: obserwacją. Wartość gdzie n to liczebność próby.

Wzór gwiazdka *


Na przykład: symbol Fˆ5(3) [ef z daszkiem] oznacza ile elementów jest mniejszych lub równych 3 w pięcioelementowej próbce.


Dystrybuanta a gęstość.
Porównajmy teraz dwa widoki: dystrybuantę rozkładu standardowego normalnego N(0, 1) oraz gęstość rozkładu standardowego normalnego N(0, 1): 

Rys. Dystrybuanta (po lewej) a gęstość (po prawej) rozkładu normalnego standardowego (LJK).
 Należy pamiętać, że obie te krzywe opisują rozkład tej samej cechy. Ktoś, kto umie czytać dystrybuantę i gęstość potrafi poradzić sobie z odpowiedzią na dowolne pytanie o rozkład posługując się jedną i drugą zmienną.


Przykład.
Rzucaliśmy osiem razy kostką do gry - więc n = 8. Otrzymaliśmy takie wyniki: 3, 5, 2, 3, 2, 1, 4, 3. Ani razu szóstki, jedna jedynka, trzy trójki, jedna czwórka i jedna piątka. Najpierw policzymy dystrybuantę dla tej sytuacji, to będzie dystrybuanta empiryczna. Korzystam ze wzoru oznaczonego gwiazdką (skrolnij wyżej).
 

Jak narysować dystrybuantę w tej sytuacji? Możemy mieć dystrybuantę empiryczną i teoretyczną. Wiemy, że skoro teoretycznie wszystkie wyniki są jednakowo prawdopodobne, to otrzymalibyśmy dystrybuanta teoretyczna. Z kolei dystrybuanta empiryczna odpowiada konkretnej sytuacji z zadania.

Rys. Dystrybuanta empiryczna i teoretyczna dla rzutu kostką (LJK).

Widać, że się różnią. Mają różne 'progi'. Wyniki to z tego, że dystrybuanta teoretyczna odpowiada sytuacji idealnej, wszystkie rzuty są jednakowo prawdopodobne, stąd te punkty skoku są takie same. Po prawej mamy dystrybuantę empiryczną, która pokazuje, co się zdarzyło. A że bywa różnie, i mimo tego, że szanse na jakiekolwiek wynik są równe, to i tak w konkretnej sytuacji otrzymujemy różne wyniki.

Uwaga.
W matematyce istnieje coś takiego jak 'dystrybucja' (nie, nie dóbr jak w ekonomii) i nie ma to nic wspólnego z 'dystrybuantą'.

piątek, 29 grudnia 2017

Sprawa Kołmogorov-Smirnov vs. Anderson-Darling

Poznamy teraz dwa testy badające, czy rozkład empiryczny jest pożądanym rozkładem. Oba te testy należą do tej samej grupy - pracują na dystrybuantach empirycznych. A dokładniej - mierzą stopień rozbieżności, między dystrybuantą empiryczną a teoretyczną. Oba jednak robią to na dwa różne sposoby.

TEST KOŁMOGOROVA-SMIRNOFFA

Zacznijmy od przykładu, a najlepsze przykłady to te życiowe. Mamy chłopców i dziewczynki - pytamy o to, w jakim momencie życiowym różnica między wzrostem chłopców a dziewczynek jest największa. W zasadzie to nawet mniej nas interesuje sam moment życiowy, jak ta maksymalna możliwa różnica.


Rys. Dystrybuanty skumulowane dla chłopców i dziewczynek (LJK).


Kołmogorova Smirnoffa właśnie bada największą odległość między krzywymi - ta różnica wynosi 0.6 (w moim przykładzie).

Badanie odległości między wykresami
W teście Kołmogorowa-Smirnoffa tak naprawdę chodzi o zbadanie odległości między dwoma wykresami funkcji: jedną jest rozkład cechy w naszej próbie, drugą jest rozkład teoretyczny, od którego podobieństwo pytamy test. Ale na moment zapomnijmy od rozkładach, testach i dystrybuantach i wróćmy do liceum lub technikum.

Przykład. Mamy dwie funkcje: jedna (liniowa)
y1 = f1(x) = x 
To jest tzw. tożsamość (identity): jedynce odpowiada jedynka, 1/2 odpowiada 1/2, pi odpowiada pi, dwójce dwójka, a zombi zombi itd.
Druga funkcja jest kwadratowa
y2 = f2(x) = x2
Dowolnej liczbie x przyporządkowuje jej kwadrat x2, np. 1/2 odpowiada 1/2 do kwadratu, czyli 1/4. 1/3 odpowiada 1/9, a 1/4 odpowiada 1/16. Zajmijmy się tylko liczbami od zera do jednego, czyli ułamkami - zaraz się okaże, dlaczego takie ograniczenie.

Pytanie o odległość między tymi funkcjami jest w zasadzie pytaniem o maksymalną różnicę między wartościami obu funkcji y1 oraz y2, jaka tylko może się zdarzyć.

Zróbmy to najpierw po kolei.
Dla x = 0 wartość pierwszej funkcji wynosi: f1(0) = 0 a wartość drugiej funkcji: f2(0) = 02=0, więc tu nie ma różnicy 0 − 0 = 0.
Dla x = 1/2 mamy w pierwszej funkcji f1(1/2) = 1/2, a w drugiej f2(1/2) = 1/4, więc różnica między nimi wynosi 1/2-1/4=1/4.
Dla x = 1/3 mamy w pierwszej funkcji f1(1/3) = 1/3, a w drugiej f2(1/3) = 1/9, więc różnica między nimi wynosi 1/3-1/9=2/9
I tak dalej, i tak dalej, i tak nieskończenie wiele razy, bo tyle jest liczb w przedziale od zera do jednego.
Gdybyśmy przejrzeli wszystkie punkty z przedziału od zera do jednego w końcu znaleźlibyśmy taki ułamek x, dla którego różnica między obiema funkcjami byłaby największa. Ale po co robić to ręcznie, skoro jest rachunek różniczkowy - istnieją matematyczne sposoby zamiast ręcznego szukania ułamek po ułamku :-).

Jeśli ktoś z liceum/technikum pamięta pochodną i jakieś porównywanie do zera, to właśnie w tym momencie zostało to wykorzystane, a my zobaczmy efekty:

Rys. Istota testu K-S, różnica między krzywymi (LJK).

Maksymalna odległość między tymi dwoma funkcjami f1(x) oraz f2(x) znajduje się dla x = 1/2 i wynosi 1/2 . Niby to samo 1/2, ale różne znaczenia. Czy to jest dużo czy mało? ‘To zależy’. Test K-S właśnie bada, czy maksymalna różnica między krzywymi (dokładniej te krzywe tu mają piękną nazwę: dystrybuantami) jest duża.
Akurat w przypadku dystrybuant można spodziewać się pewnego rodzaju wzorca takich różnic i jedne są uznawane za duże, a inne różnice za małe - dlatego to jest test statystyczny (może sprawdzić czy mamy do czynienia z typową różnicą, czy z ekstremalną <- to tak bardzo skrótowo).

Tak, jak już wspomniałam istotą testu K-S jest maksymalna różnica między dwie funkcjami reprezentującymi rozkład cechy i widać to w postaci analitycznej (czyli we wzorze) statystyki testowej:
T = sup |Fn(x) − F(x)|
gdzie T - statystyka testowa, sup oznacza, że chodzi o możliwie największą wartość różnicy między dwoma zwierzakami: Fn(x) to dystrybuanta empiryczna a F(x) - to dystrybuanta teoretyczna. Te pionowe kreski oznaczają wartość bezwzględną - interesuje nas wartość różnicy, a nie jej znak - wszystko jedno, czy minus 10 czy plus 10 stopni na dworze, mi i tak jest zimno :-)

TEST ANDERSONA-DARLINGA

Drugi test, który pomaga zbadać czy rozkład empiryczny jest rozkładem normalnym i który również w jakiś sposób wykorzystuje dystrybuanty empiryczne, to test Andersona-Darlinga. W zasadzie omówię nieco prostszą wersję testu Andersona-Darlinga, noszą kolejną nazwę pochodzącą od nazwisk twórców, tj. test Cramera-von Misesa. Zrobię to dlatego, że oba te testy: AD oraz CVM różnią się bardzo, bardzo niewiele, główna idea między nimi jest taka sama, przy czym łatwiej objaśnić CVM.
Wróćmy na chwilę do przykładu z funkcjami f1(x) (czyli liniowa) oraz f2(x) (czyli kwadratowa). Test K-S badał maksymalną różnicę między nimi. Kiedy tak patrzymy na rysunek, to przychodzi do głowy pytanie, dlaczego szukać maksymalnej różnicę - może lepiej byłoby scalić wszystkie różnice w jedną? Popatrzeć na pole między jedną a drugą krzywą?

Rys. Istota testu A-D - pole zostawione między dwiema krzywymi (LJK).

Jak obliczyć to pole zieleni? Na papierze rysujemy przerywaną linią kwadrat i rysujemy jego przekątną (zielona linia). Umieszczamy w kwadracie okrąg. Wycinamy z papieru wzdłuż przekątnej i przerywanych linii trójkąt. A potem wycinamy kawałek koła i mamy pole. Tak to mniej-więcej wygląda. Na szczęście
matematycy opracowali narzędzia, które pozwalają liczyć pole między dwoma krzywymi, bo to jest to, co w zasadzie robimy, bez wycinanek. Jeśli ktoś pamięta z edukacji szkolnej (zależnie kto, do jakiej chodził, zgodnie z reformą, ja chodziłam do prastarego liceum), do liczenia pól przydawały się całki i to jest to, co stanowi podstawę statystyki testowej w teście CVM.

gdzie T - statystyka testowa, Fn(x) to dystrybuanta empiryczna a F(x) to dystrybuanta teoretyczna. Tak, ja wiem, że to jest całka, ale powyższy wzór da radę rozumieć praktycznie następująco - oblicz pole, jakie znajduje się między wykresem jednej dystrybuanty, teoretycznej F(x) a wykresem drugiej dystrybuanty empirycznej Fn. "Oblicz powierzchnię błony jak tworzy się na dwóch dystrybuantach" - takie porównanie przychodzi mi do głowy.

Ściśle rzecz ujmując, test Andersona-Darlinga mniej więcej wygląda podobnie, poza tym wprowadza pewne ulepszenie do testu CVM - nadaje wagi tym różnicom. Dlatego też uprzedzam, że powyżej jest statystyka testowa dla CVM, dla AD nie chciałam podać, aby już nie mieszać sprawy z wagami.
Można powiedzieć, że test Andersona-Darlinga to taki upgrade testu Cramera-von-Misesa

Różnica między testem A-D a K-S
Różnicę między testem Andersona-Darlinga a testem Kołmogorova-Smirnoffa dobrze obrazuje różnica między dwie charakterystykami jeziora Loch Ness (niżej).

Rys. Jezioro Loch Ness jako ilustracja do różnicy między KS-AD (LJK).

Test Andersona Darlinga podaje powierzchnię jeziora (56.4km2) natomiast może też interesować nas maksymalna szerokość. Maksymalna szerokość Loch Ness to 1.5km (szerokość podajemy w zwykłych kilometrach, a powierzchnię w kilometrach kwadratowych). Jeden i drugi sposób daje pogląd na jezioro, i oba również tracą informację (pytanie o szerokość pomija powierzchnię jeziora; pytanie o powierzchnię jeziora pomija jego głębokość). Czerwona kreska oznacza maksymalną szerokość jeziora (oczywiście w tym ułożeniu, względem południków).

Jak to wykonać w SPSS-ie?
W SPSS-ie są dostępne tylko Kołmogorov-Smirnoff (oraz nieomówiony powyżej Shapiro-Wilk). W zasadzie to jest K-S z poprawką Lillieforsa i należy po kolei przeplikać:
Analiza→ Opis Częstości → Eksploruj 
a potem należy zaznaczyć okienko z poleceniem ‘Testy normalności wraz z wykresami’. I poszukać tabelki między wykresami (tak, ona może być w gąszczu wykresów).

czwartek, 28 grudnia 2017

Wesołego - Sylwestra i Szczęśliwego Nowego Roku!

Drodzy Czytelnicy mojego bloga, ... blożka w zasadzie, bo niektóre są strasznie obfite w treści, a ten w sumie nieduży :-)

Fot. LJK.


Z okazji Nowego Roku, który zdarzy się już w najbliższy poniedziałek, składam Wam najserdeczniejsze życzenia noworoczne. A żeby nie było, że tylko w Nowy Rok, to też i w pozostałe 365 dni również - zdrowia, szczęścia i ogólnie, żeby Wam się statystycznie lepiej wiodło niż w poprzednim. I czekolady! - Lili.

A powyżej to moja bombka choinkowa i przednia elewacja aparatu fotograficznego. Hej!


Ps. Jeśli zdarza mi się nie odpowiedzieć na wiadomość przez formularz, to trzeba wysłać drugi, przypominający :( Za wszelkie pominięcia przeproszę :-)

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 (niekoniecznie niezależnych, ale my poznajemy je w tym warunku).

Po ludzku idzie to tak: mamy zmienną X1, ma ona jakiś tam sobie rozkład, nie obchodzi nas jaki, ale grunt, że ma. Dodajemy do niej kolejną zmienną X2. Ona ma dokładnie taki sam rozkład, co pierwsza X1, ale nie wiemy, jaki. Musi pewnie być tak, że X1 + X2 ma jakiś rozkład. No, ma. To się nawet da policzyć. A jeśli dodamy X3, i dodamy X4, i ... i Xn. Z tym n to aż do nieskończoności. Wiem, że to abstrakcyjnie brzmi. To teraz podzielcie tę sumę przez n, będzie średnia. Więc matematycy zastanawiali się, jakiż to rozkład miałaby taka średnia. I wymyślili :-)
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.