Opis jednowymiarowej zbiorowości statystycznej (analiza struktury)

Rozkładem empirycznym cechy/zmiennej nazywamy przyporządkowanie wartościom zmiennej odpowiednich liczebności.

Rozkład i struktura to w zasadzie to samo. Rozkład empiryczny i szereg rozdzielczy to też to samo

Opis tabelaryczny

Szereg szczegółowy jest zamieniany na szereg rozdzielczy (punktowy bądź przedziałowy)

Szereg rozdzielczy punktowy: tablica zawierająca dwie kolumny: wartości oraz liczebności. Kolumna liczebności może zawierać liczebności lub udziały (procenty)

Może być skonstruuowany dla cech mierzalnych, bądź cech niemierzalnych wówczas, gdy liczba wariantów badanej cechy statystycznej jest niewielka. Przedstawiony w postaci tablicy statystycznej (patrz dalej), może zawierać dodatkowe informacje, jak przykładowo łączną liczbę jednostek (wiersz razem)

Przykład 1 (cecha jakościowa): Absolwenci studiów pielęgniarskich w ośmiu największych krajach UE w roku 2022

kraj liczba absolwentów
BE 5376
DE 37165
ES 11166
IT 9639
NL 10990
PL 7406
RO 1382
Razem 83124

Źródło: Eurostat, tablica Health graduates (hlth_rs_grd2)

Przykład 2 (cecha ilościowa skokowa): Gospodarstwa domowe w województwie Pomorskim wg liczby osób w 2021r

liczba osób liczba gospodarstw odsetek
o1 172351 22.34
o2 198084 25.67
o3 147032 19.06
o4 126459 16.39
o5 127676 16.55
Razem 771602 100.00

Źródło: BDL/GUS, Podgrupa P4287 (Gospodarstwa domowe według liczby osób)

Szereg rozdzielczy przedziałowy: tablica zawierająca dwie kolumny: klasy wartości oraz liczebności

Kolumna liczebności może zawierać liczebności lub udziały (procenty)

Na zasadzie wyjątku dopuszcza się aby pierwszy i ostatni przedział były otwarte, tj. nie miały dolnej (pierwszy) lub górnej (ostatni) granicy (co demonstruje przykład)

Przykład 3 (cecha ilościowa ciągła lub skokowa o dużej liczbie wartości):

Podregiony UE wg wartości współczynnika dzietności w 2022 r.

Wsp. dzietności liczba krajów
(0.8,1.05] 28
(1.05,1.3] 254
(1.3,1.55] 372
(1.55,1.8] 432
(1.8,2.05] 66
(2.05,2.3] 9
(2.3,2.55] 2
(2.55,2.8] 1
(2.8,3.05] 0
(3.05,3.3] 0
(3.3,3.55] 0
(3.55,3.8] 1
(3.8,4.05] 0
(4.05,4.3] 0
(4.3,4.55] 1
(4.55,4.8] 0

Uwagi: klasy o zerowej liczebności, w klasach 1,05–1,8 znajduje się 90% całej populacji

Przykład 4 (cecha ilościowa ciągła lub skokowa o dużej liczbie wartości, kontynuuacja):

Zwiększamy liczbę klas, zmniejszamy rozpiętość klasy, ostatni przedział jest otwarty dzięki czemu unikniemy klas z zerową liczebnością.

Wsp. dzietności liczba krajów %
(0.8,0.9] 6 0.51
(0.9,1] 11 0.94
(1,1.1] 34 2.92
(1.1,1.2] 99 8.49
(1.2,1.3] 132 11.32
(1.3,1.4] 119 10.21
(1.4,1.5] 144 12.35
(1.5,1.6] 231 19.81
(1.6,1.7] 210 18.01
(1.7,1.8] 100 8.58
(1.8,1.9] 49 4.20
(1.9,2] 14 1.20
(2,2.1] 5 0.43
(2.1,2.2] 6 0.51
(2.2,2.3] 1 0.09
(2.3,2.4] 2 0.17
(2.4 i więcej 3 0.26
Razem 1166 100.00

Podregiony czyli jednostki administracyjne na poziomie NUTS3 (takich jednostek jest 1166 (2022 rok; uwaga to się nieustannie zmienia!).

Można ewenualnie zastosować następujące wzory na zalecaną liczbę klas szeregu rozdzielczego (log oznacza logarytm dziesiętny):

\(k \leq 5 log(N)\) (15,33) albo \(k = \sqrt{n}\) (34,15)

Następnie mając ustaloną liczbę klas ustalamy rozpiętość klasy jako \(h = (x_{\max} - x_{\min})/k\)

Należy tak zbudować szereg aby pogrupowany szereg był jednomodalny oraz aby nie było przedziału z zerową liczebnością.

Na zasadzie wyjątku dopuszcza się aby pierwszy i ostatni przedział były otwarte, tj. nie miały dolnej (pierwszy) lub górnej (ostatni) granicy

Tablica statystyczna: uporządkowe dane liczbowe, podzielone na wiersze lub kolumny. Mogą być to szeregi rozdzielcze lub coś innego.

Każda tablica statystyczna musi mieć: Część liczbową (kolumny i wiersze); żadna rubryka w części liczbowej nie może być pusta (żelazna zasada); w szczególności brak danych należy explicite zaznaczyć umownym symbolem

Część opisową: tytuł tablicy; nazwy (opisy zawartości) wierszy (opcjonalne); nazwy (opisy zawartości) kolumn; wskazanie źródła danych; ewentualne uwagi odnoszące się do danych.

Wykresy statystyczne

Wykresy statystyczne są graficzną formą prezentacji materiału statystycznego, są mniej precyzyjne i szczegółowe niż tablice, natomiast bardziej sugestywne.

Wykres słupkowy

Wykres kołowy (niezalecany)

Histogram

Poniższy wykres to graficzny odpowiednik przykładu 3 (nie widać/słabo widać klasy z liczebnością równą jeden)

Jeżeli występują przedziały otwarte (przykład 4) to sprawy się zwykle komplikują (od strony technicznej). Poniżej to co wygląda na histogram jest w rzeczywistości wykresem słupkowym. Oś OX zawiera etykiety dla każdego słupka, ponieważ słupków jest względnie dużo, etykiety zachodzą na siebie…

Wykres pudełkowy

Uwaga: poniższy opis wykresu pudełkowego zawiera pojęcia jeszcze nie wyjaśnione (kwartyle, mediana, rozproszenie i IQR), które opisano dokładanie dalej, w punktach miary położenia/miary zmienności.

Konstrukcja pudełka na wykresie: górny/dolny bok równy kwartylom, a linia pozioma w środku pudełka równa medianie; linie pionowe (zwane wąsami) mają długość równą Q1 minus 1,5 IQR oraz Q3 plus IQR (Q1, Q3 to kwartyle, IQR to odstęp między kwartlowy; co to jest kwartyl, mediana i odstęp międzykwartylowy jest wyjaśnione poniżej); Linia pozioma w połowie pudełka określa przeciętny poziom zjawiska; wysokość pudełka/wąsów określa zmienność (im większe wąsy/wysokość tym większa zmienność). Obserwacje nietypowe (czyli takie których wartość jest albo mniejsza od \(Q1 - 1,5IQR\) albo większa od \(Q3 + 1,5IQR\))są zaznaczana indywidualnie jako kropki nad/pod wąsami.

Przykład:

Analiza parametryczna

Miary średnie

Miary przeciętne charakteryzują średni lub typowy poziom wartości zmiennej. Są to takie wartości, wokół których skupiają się pozostałe wartości zmiennej.

Na rysunku po lewej mamy dwa rozkłady zmiennej różniące się poziomem przeciętnym (czerwony ma przeciętnie mniejsze wartości niż turkusowy). Rozkład turkusowy jest bardziej wysmukły, co oznacza

Są to rozkłady jednomodalne, tj. wartości skupiają się wokół jednej wartości. Dla takich rozkładów ma sens obliczanie średnich.

Średnia arytmetyczna

\[\bar x = \frac{1}{N} \sum_{i=1}^N x_i\] W arkuszu kalkulacyjnym policzmy średnią za pomocą funkcji śednia. W języku R będzie to funkcja mean (albo summary)

Przykład:

Plik rugbysci.csv zawiera dane dotyczące wagi i wzrostu zawodników uczestniczących pucharze świata w rugby. Dla zorientowania się w strukturze danych zacznijmy od histogramu

r <- read.csv("rwc2015.csv", sep = ';',  header=T, na.string="NA");

p05 <- ggplot(r, aes(x =weight)) +
  ggtitle('Waga rugbystów') +
  ylab(label="liczba") +
  xlab("waga") + 
  ##geom_histogram(binwidth =bwd,color = "white", fill = "#9FE367",boundary = 0.5) +
  geom_histogram(binwidth =4, color = "forestgreen", fill = "#9FE367")
p05
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

r <- read.csv("rwc2015.csv", sep = ';',  header=T, na.string="NA");

## kolumna weight zawiera wartość brakujące (NA); przed obliczeniem
## średniej trzeba je usunąć, do czego służy opcja na.rm = TRUE
m1 <- mean(r$weight, na.rm = TRUE)
m1

Interpretacja: średnia waga zawodników uczestniczących w turnieju wynosiła 102,8 kilograma (zawsze podajemy jednostkę)

Obliczenie średniej dla danych w postaci szeregów rozdzielczych wymaga zastosowania stosownych wzorów. Po szczegóły odsyłamy do literatury.

Mediana

(Median, kwartyl drugi) dzieli zbiorowość na dwie równe części; połowa jednostek ma wartości cechy mniejsze lub równe medianie, a połowa wartości cechy równe lub większe od Mediany. Stąd też mediana bywa nazywana wartością środkową.

W arkuszu kalkulacyjnym policzymy medianę za pomocą funkcji mediana. W języku R będzie to funkcja median

me1 <- median(r$weight, na.rm = TRUE)
me1

Interpretacja: waga połowy zawodników uczestniczących w turnieju wynosiła 103,0 kilograma i mniej (zawsze podajemy jednostkę)

Obliczenie mediany dla danych w postaci szeregów rozdzielczych wymaga zastosowania stosownych wzorów. Po szczegóły odsyłamy do literatury.

Dominanta

(Mode, Moda, wartość modalna, wartość najczęstsza) jest to wartość cechy statystycznej, która w szeregu empirycznym występuje najczęściej. W szeregach prostych i rozdzielczych jest to wartość cechy, której odpowiada największa liczebność (częstość). Jeżeli liczebność populacji jest niewielka, a cecha jest ciągła to wartość dominanty może być cokolwiek przypadkowa, w szczególności jeżeli istnieje więcej niż jedna wartość naczęstsza to wtedy zamiast jednej mody mamy dwie (lub więcej); co komplikuje…

W arkuszu kalkulacyjnym policzmy dominantę za pomocą funkcji WYST.NAJCZĘŚCIEJ.WART. W języku R będzie to funkcja DescTools::Mode (funkcja Mode z pakietu DescTools)

Obliczenie mediany dla danych w postaci szeregów rozdzielczych wymaga zastosowania stosownych wzorów. Po szczegóły odsyłamy do literatury.

mode <- DescTools::Mode(r$weight, na.rm=TRUE)
mode

Interpretacja: najwięcej zawodników uczestniczących w turnieju ważyło 118 kilogramów (w tym sensie dominanta jest przypadkowa–różni się znacznie od średniej co by wskazywało na znaczną asymetrię, której jednak nie ma)

Kwartyle

(Q, quartille, \(Q_1\)/\(Q_3\)) – wartości zmiennej dla jednostek dzielących populację na cztery równe części. Kwartyl pierwszy dzieli populację w proporcji 25/75%, kwartyl drugi w proporcji 50/50%, a kwartyl trzeci w proporcji 75/25%.

kwartyle <- quantile(r$weight, na.rm=TRUE)
kwartyle
## albo
q1 <- kwartyle['25%']
q3 <- kwartyle['75%']

Interpretacja: waga 1/4 zawodników uczestniczących w turnieju wynosiła 93,0 kilograma i mniej (zawsze podajemy jednostkę); waga 3/4 zawodników uczestniczących w turnieju wynosiła 113,0 kilograma i mniej LUB waga 1/4 zawodników uczestniczących w turnieju wynosiła 113,0 kilograma i więcej.

Miary zmienności

Jak bardzo wartości zmiennej różnią się.

FYI: analiza zmienności to istota statystyki. Gdyby wszystkie jednostki w populacji były takie same, nie byłaby potrzebna ich analiza (statystyczna)

Wariancja

Wariancja (variance) jest to średnia arytmetyczna kwadratów odchyleń poszczególnych wartości cechy od średniej arytmetycznej zbiorowości.

Dlaczego kwadratów? Bo suma/średnia odchyleń wynosi zawsze zero.

Oblicznie wariancji dla szeregu prostego:

\[ s^2 = \frac{1}{N} \sum_{i=1}^N (x_i - \bar x)^2 \]

często zamiast dzielenie przez \(N\) dzielimy przez \(N-1\) (mówiąc precyzyjnie obliczając wariancję dla populacji dzielimy przez \(N\) a dla próby przez \(N-1\); w arkuszu Excel są do tego celu oddzielne funkcje).

W arkuszu kalkulacyjnym policzymy wariancję za pomocą funkcji var.S (próba/sample) lub var.P (populacja/population). W języku R będzie to funkcja var (odpowiednik var.S)

v1 <- var(r$weight, na.rm = TRUE)
v1

Dla przykładu z rugbystami, wartość wariancji wynosi 167,02 kilogramy² (kilogramy do kwadratu). Z uwagi na wariacką jednostkę, wariancji się nie interpretuje

Obliczenie wariancji dla danych w postaci szeregów rozdzielczych wymaga zastosowania stosownych wzorów. Po szczegóły odsyłamy do literatury.

Odchylenie standardowe

Pierwiastek kwadratowy z wariancji. Miara ta ma normalną jednostkę i można ją interpretować.

W arkuszu kalkulacyjnym policzymy wariancję za pomocą funkcji stdev.S (próba/sample) lub stdev.P (populacja/population). W języku R będzie to funkcja sd (odpowiednik stdev.S)

sd1 <- sd(r$weight, na.rm = TRUE)
sd1

Średnie odchylenie od średniej arytmetycznej wynosi 18,1 kilograma

Odchylenie ćwiartkowe/rozstęp ćwiartkowy

Odchylenie ćwiartkowe (Q, midhinge): \[ Q = \frac{Q_3 - Q_1}{2} \] i rozstęp ćwiartkowy (interquartile range, IQR): \[ R_Q = Q_3 - Q_1 \]

q <- (q3 - q1)/2

Interpretacja (odchylenia ćwiartkowego): przeciętne odchylenie 50% środkowych jednostek zbiorowości od mediany. W przykładzie z rugbystami: przeciętne odchylenie 50% środkowych jednostek zbiorowości od mediany wynosi 10 kilogramów.

Współczynnik zmienności

Iloraz bezwzględnej miary zmienności cechy i średniej wartości tej cechy. W analizie struktury korzysta się z różnych miar położenia i zmienności, stąd są różnie definiowane współczynniki zmienności:

\[ V_s = \frac{s}{\bar x}\\ V_Q = \frac{Q_3 - Q_1}{\textrm{Me}}\\ V_Q = \frac{Q_3 - Q_1}{Q_3 + Q_1} \]

Współczynnik zmienności jest wartością niemianowaną. Wartości liczbowe współczynników zmienności najczęściej są podawane w procentach. Przyjmuje się, że jeżeli współczynnik zmienności jest poniżej 10%, to cechy wykazują zróżnicowanie statystycznie nieistotne. Duże wartości tego współczynnika świadczą o dużym zróżnicowaniu, a więc niejednorodności zbiorowości.

Współczynnik zmienności stosuje się zwykle w porównaniach, gdy chce się ocenić zróżnicowanie: kilku zbiorowości pod względem tej samej cechy, tej samej zbiorowości pod względem kilku różnych cech.

UWAGA: ten współczynnik może dawać dziwne rezulataty jeżeli średnia wynosi zero (niezdefiniowany), jest ujemna lub jest bliska zera; zwłaszcza jeżeli użyjemy skali przedziałowej.

Miary asymetrii

Asymetria (skewness), to odwrotność symetrii. Szereg jest symetryczny jeżeli jednostki są rozłożone ,,równomiernie’’ wokół wartości średniej: \[ \bar x = \textrm{Me} = D \] Asymetria prawostronna, lewostronna; wskaźnik asymetrii (skośności), współczynniki asymetrii (skośności).

Współczynnik asymetrii oparty o moment trzeci centralny

Moment trzeci centralny to średnia arytmetyczna z podniesionych do potęgi trzeciej odchyleń wartości cechy od średniej arytmetycznej

\[ \mu_3 = \frac{1}{N} \sum_{i=1}^N (x_i - \bar x)^3 \]

Interpretacja: jeżeli \(\mu_3 = 0\) szereg symetryczny, \(\mu_3 > 0\) asymetria dodatnia (prawostronna), \(\mu_3 < 0\) asymetria ujemna (lewostronna)

Moment trzeci względny określa siłę i kierunek asymetrii: \[ g_1 = \frac{\mu_3}{s^3} \] Na podstawie badań empirycznych: \(-2 < g_1 < 2\), w skrajnych przypadkach może przekraczać ten przedział.

W arkuszu kalkulacyjnym współczynnik \(g_1\) obliczamy za pomocą funkcji skośność. W języku R będzie to funkcja DescTools::Skew (funkcja Skew z pakietu DescTools)

Współczynnik asymetrii oparty na odległościach między średnimi

\[ W_s = \frac{\bar x - D}{s} \] rzadziej używa się: \[ W_s = \frac{3(\bar x - \textrm{Me})}{s} \] Współczynnik asymetrii (skośności) oparty na odległościach między kwartylami lub decylami: \[ W_{sq} = \frac{(Q_3 - Q_2) - (Q_2 - Q_1)}{Q_3 - Q_1} \]

Przykład (rugbyści):

DescTools::Skew(r$weight, na.rm=TRUE)
ws <- (m1 - mode1) /sd1
## Albo
ws.me <- 3*(m1 - me1) / sd1

Wartość współczynnika asymetrii wynosi 0,05 (rozkład symetryczny); albo -1,18 (oparty o różnicę pomiędzy średnią a dominantą). Ponieważ dominanta jest cośkolwiek przypadkowa wartość tego współczynnika wskazuje na umiarkowaną asymetrię lewostronną. Wreszcie wartość współczynnika opartego o różnicę pomiędzy średnią a medianę wynosi 0,05 (co by wskazywało, że rozkład jest jednak symetryczny :-) )

Miary koncentracji

Koncentracja – nierównomierny podział wartości cechy w zbiorowości ALBO rozkład zbiorowości wokół wartości średniej. Ta druga koncentracja jest naszym zdaniem cokolwiek mało istotna więc ją pomijamy odsłyłając zainteresowanych do literatury.

W przypadku braku koncentracji na każdą jednostkę zbiorowości przypada taka sama część ogólnej sumy wartości (wszyscy mają tyle samo); w przypadku maksymalnej koncentracji cała suma przypada na jedną jednostkę (jeden ma wszystko reszta nic).

Współczynnik Giniego

Najczęściej stosowana miara koncentracji; dana wzorem (\(x_i\) oznacza szereg szczegółowy):

\[ G = {\frac{1}{{n^{2}}{\overline {x}}}} {\sum _{i=1}^{n}(2i-n-1)x_i } \] Jeżeli \(G=0\) to brak koncentracji; jeżeli \(G=1\) koncentracja maksymalna.

W języku R współczynnik Giniego obliczymy używając funkcji ineq::ineq (funkcja ineq z pakietu ineq)

gini <- ineq::ineq(r$weight, type = "Gini")

Dla 624 rugbystów, współczynnik koncentracji masy ciała wyniósł 0,09.

Krzywa Lorenza

Krzywa Lorenza jest funkcją określoną na zbiorze liczb dodatnich \(L(C_n)\), gdzie \(C_n\) jest kumulowaną liczebnością zaś \(L\) kumulowaną wartość cechy. Zwykle wartości kumulowane są przedstawione w procentach co pozwala na zgrabną interpretację w stylu: 20% jednostek ma 5% łącznej wartości cechy albo 50% rolników posiada 15% łacznych areałów, itp.

przekątna łącząca lewy-dolny, prawy-górny wierzchołek jest nazywana linią równomiernego rozkładu (koncentracji równomiernej, line of equality). Współczynnik Giniego to iloraz pola A do sumy pół A+B.

Im większa koncentracji, tym krzywa Lorenza jest bardziej wygięta, a wartość współczynnika Ginego większa (maksymalną wartością jest 1)

Przykład:

Wykreślenie krzywej Lorenza za pomocą funkcji ineq::Lc

plot(ineq::Lc(r$weight), main = "Krzywa Lorenza", col = "blue", lwd = 2)

Opis vs wnioskowanie

Zdarzenie losowe: zdarzenie, którego wynik nie da się przewidzieć.

Prawdopodobieństwo to miara naszej pewności co do wyniku zdarzenia losowego. Prawdopodobieństwo to liczba z przedziału [0-1]; zero oznacza. że zdarzenie nie ma szans na realizację; jeden oznacza, że zdarzenie na pewno się wydarzy.

Zmienna losowa: zmienna, której wartości realizują się z określonymi prawdopodobieństwami.

Rozkład zmiennej losowej: określenie prawdopodobieństwa dla wszystkich możliwych wartości zmiennej.

Jeżeli zbiorem wartości zmiennej jest zbiór/podzbiór liczb rzeczywistych zmienna losowa jest ciągła. Dalej zajmować się będziemy tylko zmiennymi losowymi ciągłymi.

Statystyka – dowolna funkcja \(f(x_1, x_2, … , x_n)\) próby losowej \(x_1, x_2, … ,x_n\). Ponieważ \(x_1, x_2, … ,x_n\) są zmiennymi losowymi, zatem \(f(x_1, x_2, … , x_n)\) też jest zmienną losową. Można poczynić pewne ustalenia co do rozkładu tejże statystyki

Estymator to statystyka (funkcja próby losowej), którą stosuje się, aby oszacować wartość nieznanego parametru populacji. Na przykład do oszacowania nieznanej wartości średniej w populacji możemy zastosować następującą statystykę (średnia z próby):

\[\bar x = \frac{1}{n} \sum_{i=1}^n x_i\]

Ocena parametru \(θ\), to wartość estymatora jaką przyjmuje on dla konkretnej próby.

Estymator, którego wariancja dąży do zera dla liczebności próby dążącej do nieskończoności jest estymatorem zgodnym (ze wzrostem próby rośnie dokładność)

Estymator, którego wartość oczekiwana (średnia) jest równa \(θ\), nazywamy estymatorem nieobciążonym

Estymator nieobciążony \(θ^1\) jest efektywniejszy od nieobciążonego estymatora \(θ^2\), jeżeli wariancja estymatora \(θ^1\) jest mniejsza od wariancji estymatora \(θ^2\). Estymator nieobciążony parametru \(θ\) nazywamy efektywnym (najefektywniejszym), gdy jest efektywniejszy od wszystkich estymatorów nieobciążonych tego parametru.

W ramach wnioskowania stosowane są trzy metody (podejścia):

Estymacja punktowa

Szacujemy średnią (albo inny parametr) i tę wartość uznajemy za wartość prawdziwą; dokładność szacunku jest nieokreślona. Inaczej mówiąc wartość estymatora dla konkretnej próby przyjmujemy za ocenę parametru.

Estymatorem punktowym średniej jest średnia z próby a estymatorem punktowym proporcji/ryzyka jest proporcja/ryzyko z próby.

Estymacja przedziałowa

Nie można ustalić prawdopodobieństwa popełnienia błędu dla dokładnej wartości parametru (co wynika z właściwości matematycznych modelu), ale można dla dowolnego przedziału od–do.

Czyli nie można ustalić, że z prawdopodobieństwem 95% oszacujemy wartość średnią czegoś jako 5,000000, ale można z prawdopodobieństwem 95% oszacować przedział, w którym znajdzie się średnia (przykładowo, że będzie to przedział 4,9–5,1).

Estymacja przedziałowa to oszacowanie przedziału wartości od–do, który z zadanym z góry prawdopodobieństwem zawiera prawdziwą wartość parametru.

Z góry wyznaczone prawdopodobieństwo nazywa się poziomem ufności (określa jak często mamy się NIE rąbnąć).

Testowanie hipotez

Większość analiz statystycznych polega na porównaniu. W wyniku tego porównania otrzymujemy liczbę. Załóżmy, że mamy dwie próby dotyczące wieku kandydatów na radnych do sejmików wojewódzkich z roku 2018 (średnia 46,1) oraz z roku 2014 (47,2). Różnica wynosi 1,1 lat i może być spowodowana błędem przypadkowym (tj. gdybyśmy wylosowali jeszcze raz dwie próby, to wynik byłby zupełnie odmienny np 46,9 vs 46,5) i/lub wynikać z tego, że faktycznie w roku 2014 kandydaci byli starsi.

Formalnie stawiamy hipotezę, że różnica średnich wynosi zero. Jest to tzw. hipoteza zerowa. Niezbędne jest także postawienie hipotezy alternatywnej, którą może być proste zaprzeczenie zerowej. Zapisuje się to następująco (\(m_{14}\)/\(m_{18}\) oznacza odpowiednio średnie w latach 2014/2018):

\(H_0\): różnica średnich wieku wynosi zero (\(m_{14} = m_{18}\))

\(H_1\): różnica średnich wieku jest różna od zera (\(m_{14} \not= m_{18}\))

Hipotezy sprawdzamy wykorzystując test statystyczny czyli zmienną losową której rozkład prawdopodobieństwa zależy (jest funkcją powiedziałby matematyk) od wartości testowanych parametrów (w tym przypadku \(m_{14}\) oraz \(m_{18}\)). Tę zmienną losową nazywa się statystyką testu.

Nie jest chyba wielkim zaskoczeniem, że statystyką testu w teście różnicy średnich jest różnica średnich w próbie (poprawnie mówiąc różnica uwzględniająca liczebność próby oraz zmienność obu populacji). Całkiem zdroworozsądkowo możemy przyjąć, że duże różnice statystyki testu świadczą na rzecz hipotezy alternatywnej, natomiast małe na rzecz hipotezy zerowej.

Duża różnica pomiędzy hipotezą a wynikiem z próby może wynikać z tego, że

  1. Pechowo trafiła nam się nietypowa próba, który zdarza się rzadko (rozkład normalny).

  2. Hipoteza jest fałszywa, średnie mają inną wartość niż zakładamy w hipotezie zerowej.

Statystyk zawsze wybierze drugą wersję. Pozostaje tylko ustalić (dla statystyka) co to jest rzadko?

Rzadko to z prawdopodobieństwem mniejszym niż z góry ustalone prawdopodobieństwo otrzymania różnicy (zakładając, że hipoteza zerowa jest prawdziwa), którą otrzymaliśmy w próbie lub większej (coś jak założenie, że zrealizował się najlepszy z najgorszych scenariuszy).

Przyjmijmy przykładowo, że prawdopodobieństwo wystąpienia różnicy 1,1 lat (i większej) oszacowane na podstawie odpowiedniego modelu matematycznego (rozkład normalny) wynosi 0,3 co znaczy że coś takiego zdarza się względnie często – trzy razy na 10 pobranych prób.

Załóżmy z kolei że, ta różnica wyniosła 3,2 lata. Prawdopodobieństwo wystąpienia takiej różnicy (i większej) wynosi 0,009 co znaczy że coś takiego zdarza się względnie rzadko – 9 razy na tysiąc prób.

Przyjmując, że możemy się mylić 5 razy na 100 w pierwszym przypadku statystyk powie, że nie ma podstaw do odrzucenia hipotezy \(H_0\). Różnica 1,1 lat wynika z przypadku. W drugim wypadku statystyk powie, że hipoteza jest fałszywa, bo zdarzyło się coś co nie powinno się zdarzyć.

Prawdopodobieństwo „graniczne” ustalamy z góry i nazywa się ono poziomem istotności. Określa ono jak często możemy się rąbnąć odrzucając hipotezę zerową, która jest prawdziwa.

Ale jest jeszcze drugi przypadek popełnienia błędu: przyjmujemy hipotezę zerową, która jest fałszywa. W testach statystycznych nie określa się prawdopodobieństwa popełnienia tego błędu, a w związku z tym nie można przyjąć hipotezy zerowej (bo nie znamy ryzyka popełnienia błędu).

W konsekwencji hipotezę zerową albo się odrzuca albo nie ma podstaw do odrzucenia. Wniosek cokolwiek niekonkluzywny, ale tak jest.

Dlatego też często „opłaca się” tak postawić hipotezę zerową aby ją następnie odrzucić, bo taki rezultat jest bardziej konkretny.

Testy nieparametryczne

Można testować hipotezy na temat wartości parametrów, ale można też testować przypuszczenia o charakterze mniej konkretnym. Na przykład, że dwie zmienne są niezależne (co to znaczy wyjaśniono w następnym rozdziale), albo że dwa rozkłady są podobne do siebie (rozkłady nie średnie). Takie hipotezy/testy określa się jako nieparametryczne. Przykładami są testy niezależności chi-kwadrat albo normalności Shapiro-Wilka (opisane w następnym rozdziale)

Oczywiste, ale powtórzmy: przypuszczenia o charakterze nieparametrycznym możemy tylko testować (sprawdzać hipotezy); nie obliczamy wtedy ani ocen ani nie wyznaczamy przedziałów ufności.

Analiza współzależności

Problem

Związki między cechami mogą być: funkcyjne (nauki przyrodnicze) – wartościom jednej zmiennej odpowiada tylko jedna wartość drugiej zmiennej lub stochastyczne – wartościom jednej zmiennej odpowiadają z pewnym przybliżeniem wartości innej zmiennej.

Problem: czy istnieje związek (zależność) pomiędzy zmiennymi? Jaki jest charakter zależności? Jaka jest siła zależności? Przykładowo czy istnieje związek pomiędzy wielkością dochodu (przyczyna) a wielkością spożycia mięsa (skutek), albo jako jest zależność pomiędzy wielkością produkcji, nakładami kapitałowymi i wielkości nakładów pracy?

Zależność stochastyczna występuje wtedy, gdy wraz ze zmianą wartości jednej zmiennej zmienia się rozkład prawdopodobieństwa drugiej zmiennej.

Zależność korelacyjna występuje wtedy gdy gdy wraz ze zmianą wartości jednej zmiennej zmienia się wartość średnia drugiej zmiennej.

Wykresy: wykres rozrzutu

W układzie kartezjańskim każdej obserwacji odpowiada kropka o współrzędnych XY.

O występowaniu związku świadczy układanie się kropek według jakiegoś kształtu (krzywej). O braku związku świadczy chmura punktów niepodobna do żadnej krzywej.

Punkty układające się według prostej świadczą o zależności liniowej (wyjątek: linia pozioma lub pionowa) Punkty układające się według krzywej świadczą o zależności nieliniowej.

Przykład: Zależność pomiędzy zamożnością (gdp per capita) a współczynnkiem dzietności (fert) na poziomie 229 makroregionów (NUTS2) w Unii Europejskiej.

Przykład: Zależność pomiędzy odsetkiem kobiet w wieku 15-64 lata posiadających wyższe wykształcenie (edu) a współczynnkiem dzietności (fert) na poziomie 229 makroregionów (NUTS2) w Unii Europejskiej.

Tablice wielodzielne

Test niezależności chi-kwadrat

Współczynnik korelacji liniowej Pearsona

Kowariancja to średnia arytmetyczna iloczynów odchyleń wartości zmiennych \(X\) oraz \(Y\) od ich wartości średnich. Dla \(n\) obserwacji na zmiennych \(X\) oraz \(Y\) można to zapisać w postaci następującej formuły:

\[\mathrm{cov} (xy) = \frac{\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{n}\]

Kowariancja zależy od rozproszenia (im większe tym większa), ma też dziwną jednostkę (jednostkaX · jednostkaY) oraz zależy od wybranych skal (tony vs gramy na przykład.)

Z powyższych powodów do pomiaru związku pomiędzy cechami używa się standaryzowanego współczynnika kowariancji, zwanego współczynnikiem korelacji liniowej, (Pearson linear correlation coefficient). Standaryzacja polega na podzieleniu wartości kowariacji przez iloczyn odchyleń standardowych \(s_x\) oraz \(s_y\).

\[r_{xy} = \frac{\mathrm{cov}(xy) }{s_x \cdot s_y}\]

Współczynnik jest miarą niemianowaną, przyjmującą wartości ze zbioru \([-1;1]\); Skrajne wartości \(\pm 1\) świadczą o związku funkcyjnym (wszystkie punkty układają się na linii prostej); wartość zero świadczy o braku związku (linia pozioma/pionowa).

Interpretacja opisowa: wartości powyżej 0,9 świadczą o silnej zależności.

Kwadrat współczynnika korelacji liniowej nazywamy współczynnikiem determinacji.

W arkuszu kalkulcyjnym Excel współczynnik korelacji liniowej obliczmy za pomocą funkcji PEARSON(x; y). W języku R służy do tego funkcja cor.

Przykład: Zależność pomiędzy zamożnością a dzietnościa na poziomie makroregionów 229 (NUTS2) w Unii Europejskiej

r1 <- cor(fgwe$fert, fgwe$gdp, method = "pearson")

wartość współczynnika korelacji wynosi -0.001 (brak korelacji)

Współczynnik korelacji liniowej Spearmana

Współczynnik korelacji rang (Spearmana vel Spearman’s Rank-Order Correlation) może być stosowany w przypadku gdy cechy są mierzone w skali porządkowej (lub lepszej)

Obliczenie współczynnika Spearmana dla \(N\) obserwacji na zmiennych \(XY\) polega na zamianie wartości zmiennych \(X\) oraz \(Y\) na rangi (numery porządkowe od \(1...N\)). Następnie stosowana jest formuła współczynnika korelacji liniowej Pearsona (\(\tau_x\) oraz \(\tau_y\) oznaczają rangi):

\[\rho_{xy} = \frac{\textrm{cov}(\tau_x, \tau_y)}{s_{\tau_x} s_{\tau_y}}\]

Współczynnik \(\rho_{xy}\) to – podobnie jak oryginalny współczynnik korelacji liniowej Pearsona – miara niemianowana, o wartościach ze zbioru [-1;1];

Przykład

Regresja liniowa

Regresja liniowa zakłada, że istnieje związek przyczyna-skutek i ten związek można opisać linią prostą (stąd liniowa). Skutek jest jeden i nazywa się go zmienną zależną zaś przyczyn może być wiele i noszą nazwę zmiennych niezależnych (albo predyktorów). W przypadku gdy związek dotyczy dwóch zmiennych mówi się o regresji prostej.

Jeżeli zmiennych niezależnych jest więcej niż jedna, to mówimy o regresji wielorakiej.

Regresja prosta

Równanie regresji dla zmiennych \(Y\) (skutek) oraz \(X\) (przyczyna) można zapisać następująco:

\[Y = b_0 + b_1 \cdot X + e \]

\(Y = b_0 + b_1 \cdot X\) to część deterministyczna, a \(e\) oznacza składnik losowy.

O tym składniku zakładamy, że średnia jego wartość wynosi zero. Można to sobie wyobrazić, że w populacji jest jakaś prawdziwa zależność \(Y = b_0 + b_1 \cdot X\) pomiędzy \(X\) a \(Y\), która w próbie ujawnia się z błędem o charakterze losowym. Ten błąd może wynikać z pominięcia jakiejś ważnej zmiennej (model to zawsze uproszczenie rzeczywistości), przybliżonego charakteru linii prostej jako zależności pomiędzy \(X\) a \(Y\) (prosta ale nie do końca prosta) albo błędu pomiaru.

Współczynnik \(b_1\) (nachylenia prostej) określa wielkość efektu w przypadku regresji, tj. siły zależności pomiędzy zmiennymi.

Współczynnik \(b_1\) ma prostą interpretację: jeżeli wartość zmiennej \(X\) rośnie o jednostkę to wartość zmiennej \(Y\) zmienia się przeciętnie o \(b_1\) jednostek zmiennej \(Y\).

Wyraz wolny zwykle nie ma sensownej interpretacji (formalnie jest to wartość zmiennej \(Y\) dla \(X=0\))

Oznaczmy przez \(y_i\) wartości obserwowane (zwane też empirycznymi) a przez \(\hat y_i\) wartości teoretyczne (leżące na prostej linii regresji).

Wartości \(b_0\) oraz \(b_1\) wyznacza się minimalizując sumę kwadratów odchyleń wartości teoretycznych od wartości empirycznych, tj.:

\[\sum_{i=1}^n (y_i - \hat y_i)^2 = \sum_{i=1}^n e_i^2 \to \min\]

Rozwiązując powyższy problem minimalizacyjny otrzymujemy wzory definiujące parametry \(b_0\) oraz \(b_1\). Metoda wyznaczania parametrów linii prostej w oparciu o minimalizację sumy kwadratów odchyleń nosi nazwę metoda największych kwadratów.

Przypominamy, że estymatorem nazywamy metodę oszacowania parametru na podstawie próby. Ponieważ traktujemy \(b_0\) oraz \(b_1\) jako parametry jakieś populacji generalnej to wzory na \(b_0\) oraz \(b_1\) statystyk nazwie estymatorami parametrów \(b_0\) oraz \(b_1\). W konsekwencji tego \(b_0\)/\(b_1\) posiadają jakąś wartość średnią oraz wariancję.

Przypominamy że wartość średnia dobrego estymatora powinna wynosić zero (bo wtedy nie ma błędu systematycznego) oraz że wariancja estymatora powinna maleć wraz ze wzrostem liczebności próby. Można udowodnić że estymatory parametrów \(b_0\)/\(b_1\) uzyskane metodą najmniejszych kwadratów posiadają obie właściwości.

Graficznie kryterium minimalizacyjne przedstawia rysunek

Suma podniesionych do kwadratu odległości pomiędzy czerwonymi i niebieskimi kropkami ma być minimalna. Kropki niebieskie to wartości empiryczne; kropki czerwone to wartości teoretyczne. Zadanie wyznaczenie parametrów takiej prostej sprowadza się do (dla dwóch zmiennych):

\[\sum (y_i - b_0 - b_1 x)^2 \to \min\]

Rozwązując powyższe ze względu na \(b_0\) i \(b_1\) otrzymujemy oceny parametrów regresji:

\[b_1 = \frac{\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^n (x_i - \bar x)^2}\]

\[b_0 = \bar y - b\bar x\]

Można udowodnić, że bez względu czy punkty na wykresie układają się w przybliżeniu wzdłuż prostej czy nie, zawsze jakaś prosta zostanie dopasowana (jeżeli tylko punktów jest więcej niż jeden.) Jak to ocenić w sposób bardziej konkretny a nie tylko na oko dopasowanie prostej do wartości empirycznych?

Ocena dopasowania: wariancja resztowa oraz średni błąd szacunku

Oznaczając resztę jako: \(e_i = y_i - \hat y_i\), definiujemy wariancję resztową jako:

\[s_e^2 = \frac{1}{n-k} \cdot \sum e_i^2\].

Gdzie \(n\) oznacza liczbę obserwacji (liczebność próby), a \(k\) liczbę szacowanych parametrów bez wyrazu wolnego czyli jeden w regresji prostej (a więcej niż jeden w regresji wielorakiej o czym dalej.)

Pierwiastek kwadratowy z wariancji resztowej. nazywamy średnim błędem szacunku (mean square error, MSE)

Ocena dopasowania: współczynniki zbieżności i determinacji

Suma kwadratów reszt (albo odchyleń wartości teoretycznych od wartości empirycznych, albo suma kwadratów błędów vel resztowa suma kwadratów):

\[\mathrm{RSK} = (y_1 - \hat y_1)^2 + (y_2 - \hat y_2)^2 + ... + (y_n - \hat y_n)^2\].

Suma kwadratów odchyleń wartości empirycznych od średniej (ogólna suma kwadratów):

\[\mathrm{OSK} = (y_1 - \bar y)^2 + (y_2 - \bar y)^2 + ... + (y_n - \bar y)^2\]

Suma kwadratów odchyleń wartości teoretycznych od średniej (wyjaśniona suma kwadratów):

\[\mathrm{WSK} = (\hat y_1 - \bar y)^2 + (\hat y_2 - \bar y)^2 + ... + (\hat y_n - \bar y)^2\]

Można wykazać, że \(\mathrm{OSK} = \mathrm{WSK} + \mathrm{RSK}\) zatem (po podzieleniu obu stron równania przez \(\mathrm{OSK}\) otrzymujemy:

\[ 1 = \mathrm{WSK}/\mathrm{OSK} + \mathrm{RSK}/\mathrm{OSK}\]

Współczynnik determinacji oznaczany jako \(R^2\) to \(\mathrm{WSK}/\mathrm{OSK}\).

Współczynnik zbieżności oznaczany jako \(\Phi^2\) (duża grecka litera Fi) to \(\mathrm{RSK}/\mathrm{OSK}\).

Współczynnik zbieżności $^2$100% informuje ile % zmienności zmiennej \(Y\) nie została wyjaśniona przez zmienne regresji.

Współczynniki przyjmują wartość z przedziału \([0,1]\) lub \([0, 100]\)% jeżeli ich wartości zostaną pomnożone przez 100.

Interpretacja współczynnika zbieżności: udział (procent) zmienność wyjaśnianej przez linię regresji. Im \(R^2\) jest bliższe jedności (lub 100% jeżeli jest współczynnik zbieżności jest wyrażony w procentach) tym lepiej.

Odchylenie standardowe składnika resztowego \(S_e\)

\[S_e = \sqrt{\frac{1}{n-2} \cdot \sum e_i^2 }\]

O ile średnio (\(\pm\)) odchylają się wartości rzeczywiste zmiennej objaśnianej od teoretycznych obliczonych na podstawie regresji.

Współczynnik zmienności resztowej \(V_e\)

\[V_e = \frac{S_e}{\bar y} \cdot 100\]

Jaką część średniej wartości zmiennej objaśnianej stanowi jej odchylenie standardowe reszt.

Ocena dopasowania: istotność parametru \(b_1\)

Jeżeli: \(Y= 0 \cdot X + b_0\), to \(Y = b_0\) czyli nie ma zależności pomiędzy \(X\) oraz \(Y\). Wartości \(b_1\) bliskie zero wskazują na słabą zależność pomiędzy cechami.

Przypominamy, że estymator parametru \(b_1\) ma średnią równą prawdziej wartości \(b_1\). Dodatkowo zakładamy, że rozkład tego estymatora jest normalny. To założenie pozwala wiarygodnie oszacować wariancję; w konsekwencji znamy dokładny rozkład (bo przypominamy, że rozkład normalny jest określony przez dwa parametry: średnią oraz właśnie wariancję)

Można teraz zadać pytanie jeżeli faktycznie \(b_1=0\), to jakie jest prawdopodobieństwo, że współczynnik \(\hat b_1\) oszacowany na podstawie \(n\) obserwacji będzie (co do wartości bezwzględnej) większy niż \(b_e\). Albo inaczej: otrzymaliśmy \(b_e\), jakie jest prawdopodobieństwo otrzymania takiej wartości (lub większej co do wartości bezwzględnej) przy założeniu, że istotnie \(b_1=0\).

Jeżeli takie prawdopodobieństwo jest duże, to uznajemy, że być może \(b_1 = 0\), a jeżeli małe to będziemy skłonni uznać, że \(b_1 \not= 0\). Duże/małe przyjmujemy arbitralnie, zwykle jest to \(0,1\), \(0,05\) lub \(0,01\). Tak zgadza się, to prawdopodobieństwo to poziom istotności

Testowanie istotności współczynnika regresji jest ważnym kryterium oceny jakości dopasowania. Regresja z nieistotnym współczynnikiem nie może być podstawą do interpretowania zależności pomiędzy \(X\) oraz \(Y\).

Przykład: Zależność pomiędzy odsetkiem kobiet w wieku 15–64 lata posiadających wyższe wykształcenie (edu) a współczynnkiem dzietności (fert) na poziomie 229 makroregionów (NUTS2) w Unii Europejskiej

lm1 <- lm(fert ~ edu, data=fgwe)
summary(lm1)
## 
## Call:
## lm(formula = fert ~ edu, data = fgwe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59914 -0.16875  0.01085  0.13729  2.07026 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  1.597315   0.057328  27.863 <0.0000000000000002 ***
## edu         -0.003775   0.001728  -2.185              0.0299 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2575 on 227 degrees of freedom
## Multiple R-squared:  0.0206, Adjusted R-squared:  0.01629 
## F-statistic: 4.775 on 1 and 227 DF,  p-value: 0.0299

Przypadek ogólny: regresja wieloraka

Uogólnieniem regresji prostej jest regresja wieloraka. W modelu regresji wielorakiej po lewej stronie równania występuje zmienna liczbowa oznaczona jako \(Y\), a po prawej zmienne liczbowe lub nominalne, \(X_1, \ldots, X_k\):

\[Y = b_0 + b_1 \cdot X_1 + b_2 \cdot X_2 + ... + b_k \cdot X_k \]

Wpływ każdej zmiennej \(X_i\) na zmienną zależną \(Y\) jest określony przez odpowiedni współczynnik \(b_i\).

Podobnie jak w przypadku regresji prostej do oceny stopnia dopasowania modelu do danych wykorzystuje się: średni błąd szacunku, współczynnik zbieżności \(R^2\) oraz weryfikuje się istotność współczynników \(b_i\).

Przykład

Przykład: Zależność pomiędzy odsetkiem kobiet w wieku 15-64 lata posiadających wyższe wykształcenie (edu) oraz wskaźnikiem deprywacji materialnej i społecznej (sd) a współczynnkiem dzietności (fert) na poziomie 229 makroregionów (NUTS2) w Unii Europejskiej

lm2 <- lm(fert ~ edu + sd, data=fgwe)
summary(lm2)
## 
## Call:
## lm(formula = fert ~ edu + sd, data = fgwe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67767 -0.16876  0.00823  0.14872  1.97401 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  1.425408   0.068082  20.937 < 0.0000000000000002 ***
## edu         -0.001169   0.001770  -0.661                 0.51    
## sd           0.013030   0.003018   4.317            0.0000237 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.248 on 226 degrees of freedom
## Multiple R-squared:  0.09522,    Adjusted R-squared:  0.08721 
## F-statistic: 11.89 on 2 and 226 DF,  p-value: 0.00001229

Predykcja

Uwagi wstępne

Przewidywanie co będzie. Anglicy odróżniają forecasting (prognozowanie) or predicting (przewidywanie); w sumie obie rzeczy dotyczą przeszłości ale istnieje subtelna różnica:

jaka będzie wielkość GDP w przyszłym roku (forecasting)

jaka będzie wielkość produkcji jeżeli zwiększymy nakłady pracy o \(x\) jednostek (predicting)

Prognozowanie w modelu regresji

Załóżmy, że:

\[y_i = b_1 + b_2 x_{2} + b_3 x_{3} + \ldots + b_k x_{k} + u_i\]

prognozowana wartość:

\[y^* = b_1 + b_2 x_{2}^* + b_3 x_{3}^* + \ldots + b_k x_{k}^* + u^*\]

gdzie: \(x_{2}^*\ldots x_{k}^*\) jest wektorem wartości przyjmowanych przez zmienne objaśniające, a \(u^*\) jest składnikiem losowym o takim samym rozkładzie jak składniki losowe z próby (\(N(0, \sigma)\))

Prognozy dokonujemy obliczając

\[\hat y^* = \hat b_1 + \hat b_2 x_{2}^* + \hat b_3 x_{3}^* + \ldots + \hat b_k x_{k}^*\]

gdzie wartości \(\hat b_1, \ldots, \hat b_k\) oznaczają oceny parametrów

\(\hat y^*\) nazywamy prognozą punktową (w odróżnieniu od przedziałowej o czym za chwilę)

Wówczas błąd prognozy wynosi:

\[e^* = y^* - \hat y^*\]

Oczywiście wartość \(y^*\) jest nieznana, ale można udowodnić, że

\[E(e^*) = 0\]

(średnia wartość błędu wynosi zero)

oraz, że wariancja jest równa:

\[{\rm var} (e^*) = {\rm var} (\hat y^*) + {\rm var}(u) = {\rm var} (\hat y^*) + \sigma\]

gdzie \(\sigma = {\rm var}(u)\) (dla zwięzłości)

Wariancja błędu prognozy jest więc większa od wariancji składnika losowego. Prognozujemy z błędem \({\rm var} (\hat y^*)\) plus zjawisko ma charakter losowy także w przyszłości \(\sigma = \sqrt{{\rm var}(u)}\). Te dwa rodzaje błędów się sumują…

Ponieważ z założenia \(u^*\) jest składnikiem losowym o takim samym rozkładzie jak składniki losowe z próby \(u^*\), zatem \(\sigma\)jest jedynym nieznanym parametrem (inaczej mówiąc \({\rm var} (\hat y^*)\) jest jakąś funkcją \(\sigma\); jaką można to analitycznie ustalić ale nie będziemy tego robić)

Wstawiając zamiast \(\sigma\) jego estymator \(s_e^2\) (zwany błędem standardowy reszt albo odchyleniem standardowym składnika losowego – patrz wyżej) otrzymujemy nieobciążony estymator wariancji błędu prognozy \({\rm var} (e^*) = s^{*2}\)

Pierwiastek kwadratowy z wariancji błędu prognozy \(s^*\) nosi nazwę średniego błędu predykcji.

Iterpretacja: ile średnio wyznaczona prognoza może się odchylać od wartości rzeczywistych zmiennej prognozowanej.

Ponieważ \(u\) ma rozkład \(N(0, \sigma)\), to m.in. duże odchylenia są mniej prawdopodobne niż małe można zatem wyznaczyć dwie wartości:

\[\hat y^*_l < y < \hat y^*_u\]

pomiędzy którymi – z zadaną dokładnością – znajduje się prognozowana wartość. Ta zadana dokładość jest określona prawdopodobieństwem zwykle 0,95 lub 0,99.

Taka konstrukcja nazywa się prognozą przedziałową.

Wniosek: jeżeli model jest słabo dopasowany to prognoza jest jeszcze gorsza a prognoza przedziałowa może być mało przydatna, typu 20 plus/minus 30.

Miary dokładności predycji

Niewątpliwie średniego błędu predykcji jest miarą jej dokładności. Im mniejszy tym lepiej.

Średni błąd predykcji należy do miar dokładności predykcji ex-ante, tj. wyliczanych bez znajomości prawdziwych wartości (prognoza się jeszcze nie zrealizowała)

Oprócz miary ex-ante są miary ex-post tj. takie w których prognozy porównuje się z wartościami zrealizowanymi.

Żeby te wartości zrealizowane mieć to albo trzeba poczekać :-) Albo dokonać następującego tricku: dzielimy dane na dwie części zwane zwykle zbiorem uczącym (training set) i zbiorem testowym (test set)

Na podstawie zbioru uczącego szacujemy model.

Na podstawie zbioru testowego sprawdzamy właściwości predyktywne modelu.

Ponieważ dysponujemy prognozami oraz realizacjami to możliwe jest oszacowanie jakości prognoz. Są to tego stosowane następujące miary:

Średniokwadratowy błąd prognozy (albo średni błąd kwadratowy; MSE z angielska tj mean square error)

\[{\rm MSE} = \frac{1}{n^*}\cdot \sum_{i=1}^n (y_i - \hat y_i)^2\]

Pierwiastek MSE zwany RMSE (root MSE):

\[{\rm RMSE} = \sqrt{ {\rm MSE }}\]

Średni błąd bezwzględny (Mean Absolute Error; MAD):

\[{\rm MAD} = \frac{1}{n^*}\cdot \sum_{i=1}^n |y_i - \hat y_i|\]

Interpretacja: ile średnio wyznaczona prognoza odchylała się od wartości rzeczywistych zmiennej prognozowanej (było może się odchylać dla miar ex-ante)

Analiza szeregów czasowych

Pojęcia wstępne

Szereg czasowy to ciąg pomiarów z kolejnych momentów/okresów czasu.

Szereg czasowy jest szeregiem okresów jeżeli dotyczy strumieni (zasobów). Przykładowo liczba urodzeń żywych w powiecie kwidzyńskim w latach 2010–2019, to szereg okresów. W szczególności można dodać liczbę urodzeń w poszczególnych latach otrzymując łączną liczbę urodzonych w tym okresie.

Urodzenia żywe w powiecie kwidzyńskim w latach 2010–2019

Rok 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
L. urodzeń 965 937 907 821 831 785 813 861 834 833

Źródło: https://bdl.stat.gov.pl/

W latach 2010–2019 w powiecie kwidzyńskim urodziło się 965 + 937 + … + 833 = 8587 dzieci (średnio 858,7 rocznie).

Albo szereg czasowy jest szeregiem momentów jeżeli dotyczy stanów. Przykładowo liczba ludności powiatu kwidzyńskiego w latach 2010–2019 (stan na 31.12), to szereg czasowy momentów.

Ludność ogółem w powiecie kwidzyńskim w latach 2010–2019

Rok 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
L.ludności 83562 83734 83783 83611 83627 83464 83495 83423 83291 83174

Źródło: https://bdl.stat.gov.pl/

(Zwróćmy uwagę że taki szereg zwykle w tytule ma magiczną formułę stan na dzień-miesiąc – w przykładzie stan na ostatni dzień roku) Nie można sumować liczby mieszkańców z kolejnych lat, bo taka suma nie ma merytorycznie sensu. Można podać średnią (przeciętny stan) ludności jako wartość średniej chronologicznej:

\[\bar x_{ch} = \frac{(\frac{1}{2} x_1 + x_2 + ... + x_{n-1} + \frac{1}{2} x_{n} )}{n-1} = (83562/2 + 83734 + ... + 83291 + 83174/2)/9 = 83532,9\]

Średnia roczna liczba ludności w powiecie kwidzyńskim w latach 2010–2019 wyniosła zatem 83532,9 mieszkańców.

Szeregi czasowe mogą różnić się częstotliwością (dzienna, tygodniowa, miesięczna, kwartalna, roczna).

Wykresy: liniowe i słupkowe

Jeżeli oś OY nie zaczyna się od zera to krzywa będzie bardziej stroma bo efektywnie wykres będzie miał większe proporcje wysokość/szerokość:

Przyrosty i indeksy

Oznaczmy wartość zmiennej \(Y\) w okresie \(t\) jako \(y_t\) (mówi się okres badany) oraz wartość zmiennej \(Y\) w okresie \(0\) jako \(y_0\) (okres podstawowy).

Przyrost absolutny to po prostu \(y_t - y_0\). Jeżeli obliczamy przyrosty dla wielu okresów, to można je porównywać do jednego okresu/momentu (\(y_1\) w zapisie poniżej):

\[y_2 - y_1, y_3 - y_1, \ldots, y_{n-1} - y_1, y_{n} - y_1\]

lub dla poprzedniego okresu/momentu

\[y_2 - y_1, y_3 - y_2, \ldots, y_{n-1} - y_{n-2}, y_{n} - y_{n-1}\]

W pierwszym przypadku mówimy o przyrostach jednopodstawowych, w drugim zaś o przyrostach łańcuchowych. Zwróćmy uwagę, że podstawa nie musi być równa pierwszej wartości. Może to być np. wartość ostatnia lub środkowa:

\[y_1 - y_n, y_2 - y_n, \ldots, y_{n-2} - y_n, y_{n-1} - y_n\]

Interpretacja: o ile wzrosła/spadła wielkość zjawiska w okresie badanym w porównaniu do podstawowego w jednostkach zmiennej \(Y\) oczywiście.

Przyrost względny to \((y_t - y_0)/y_0\) Jeżeli obliczamy przyrosty dla wielu okresów, to podobnie jak miało to miejsce w przypadku przyrostów absolutnych możemy się posługiwać
indeksami jednopodstawowymi lub łańcuchowymi.

Przyrost względny zwykle wyrażone są w procentach jeżeli \(y_t = 115\) a \(y_0 = 100\), to \((y_t - y_0)/y_0 \cdot 100 =115\). Zjawisko ma o 15% wyższy poziom w roku \(t\).

Wskaźniki dynamiki (indeksy) określają względną zmianę wartości szeregów czasowych (popularny procent zmiany), czyli \(y_t/y_0\). Jeżeli podstawą porównań jest okres poprzedni \(y_t/y_{t-1}\) to indeks nazywany jest łańcuchowym, jeżeli podstawa porównań jest stała, tj. \(y_t/y_{c}\), to indeks nazywany jest jednopodstawowy.

Indeks zwykle wyrażone są w procentach jeżeli \(y_t = 115\) a \(y_0 = 100\), to $y_t/y_0 =15 $. Zjawisko ma o 15% wyższy poziom w roku \(t\).

Używamy średniej geometrycznej do obliczenia średniego tempa zmian zjawiska w czasie:

\[\bar y_g = \sqrt[n-1]{y_n/y_1}\]

Przykład: Zwiedzający Muzeum Zamkowe w Malborku w latach 2015–2016

rok-miesiąc zwiedzający (1) (2) (3) (4) (5) (6)
2015-01 3278 NA 0 NA 0.00 NA 100.00
2015-02 4149 871 871 26.57 26.57 126.57 126.57
2015-03 4623 474 1345 11.42 41.03 111.42 141.03
2015-04 13083 8460 9805 183.00 299.12 283.00 399.12
2015-05 68485 55402 65207 423.47 1989.23 523.47 2089.23
2015-06 66199 -2286 62921 -3.34 1919.49 96.66 2019.49
2015-07 105843 39644 102565 59.89 3128.89 159.89 3228.89
2015-08 95255 -10588 91977 -10.00 2805.89 90.00 2905.89
2015-09 28021 -67234 24743 -70.58 754.82 29.42 854.82
2015-10 13532 -14489 10254 -51.71 312.81 48.29 412.81
2015-11 4382 -9150 1104 -67.62 33.68 32.38 133.68
2015-12 3204 -1178 -74 -26.88 -2.26 73.12 97.74
2016-01 3382 178 104 5.56 3.17 105.56 103.17
2016-02 6772 3390 3494 100.24 106.59 200.24 206.59
2016-03 5361 -1411 2083 -20.84 63.54 79.16 163.54
2016-04 18435 13074 15157 243.87 462.39 343.87 562.39
2016-05 73453 55018 70175 298.44 2140.79 398.44 2240.79
2016-06 70743 -2710 67465 -3.69 2058.11 96.31 2158.11
2016-07 109008 38265 105730 54.09 3225.44 154.09 3325.44
2016-08 116078 7070 112800 6.49 3441.12 106.49 3541.12
2016-09 37274 -78804 33996 -67.89 1037.10 32.11 1137.10
2016-10 14841 -22433 11563 -60.18 352.75 39.82 452.75
2016-11 5416 -9425 2138 -63.51 65.22 36.49 165.22
2016-12 3369 -2047 91 -37.80 2.78 62.20 102.78

Zawartość kolumn: (1) przyrosty łańcuchowe, (2) przyrosty jednopodstawowe (2015-01=100%), (3) przyrosty względne łańcuchowe, (4) przyrosty względne jednopodstawowe (2015-01=100%), (5) indeksy łańcuchowe, (6) indeksy jednopodstawowe (2015-01=100%).

Indeksy agregatowe

Model wahań w czasie

W szeregu czasowym można zwykle wyróżnić długookresową tendencję (trend); powtarzalne wahania (sezonowość); resztę traktuje się jako wartości przypadkowe. Reasumując:

\[TS = T + S + E\]

lub

\[TS = T \cdot S \cdot E\]

Pierwszy wariant nazywa się addytywny drugi multiplikatywny. W wariancie addytywnym zmiany (trendu/sezonowości) okres/okres są stałe; w wariancie multiplikatywnym tempo zmiany jest stałe, tj. zjawisko okres/okres rośnie/spada o x%. W jednostkach bezwzględnych oznacza to, że rośnie/spada coraz szybciej.

Problem: oszacowanie \(T\) oraz \(S\)

Szacowanie trendu

Metoda mechaniczna MA

Średnia ruchoma (moving average MA). Idea tego wygładzania jest prosta: sumujemy kolejne wartości szeregu i dzielimy przez liczbę elementów sumy (średnia \(k\)-okresowa); Ile elementów sumujemy jest dobierane metodą prób/błędów…

średnia trzy okresowa

\[\bar y_{n-1} = (y_{n-2} + y_{n-1} + y_n)/3\]

\[\bar y_{n-2} = (y_{n-4} + y_{n-1} + y_n)/5\] scentrowana czterookresowa:

\[\bar y_{n-2} = \frac{\frac{1}{2} y_{n-4} + y_{n-3} + y_{n-2} +y_{n-1} + \frac{1}{2} y_n }{4}\]

Przykład dzienne dane nt. liczby zgonów z powodu COVID (w okresie 1.10.2020–5.2.2021; źródło komunikaty MZ via Twitter a od 28.01.2021 https://www.gov.pl/web/koronawirus/wykaz-zarazen-koronawirusem-sars-cov-2; strona reklamowana przez Google!):

Na powyższym wykresie mamy 4 średnie ruchome 3, 7, 14, 28 okresową. Średnia trzy okresowa jest za mało wygładzona. Średnie 7, 14, 28 są podobne ale każda kolejna jest krótsza.

Najlepsza średnia ruchoma: 7 okresów. (Dostatecznie wygładza trend i jest najdłuższa)

Metoda analityczna

Polega na dopasowaniu określonej funkcji matematycznej; w najprostszym przypadku prostej przy użyciu metody najmniejszych kwadratów, czyli zakłada się że trend jest postaci:

\[Y = a + b \cdot t + e\]

gdzie: \(e\) oznacza składnik losowy;

Parametry \(a\) i \(b\) są wyznaczane w taki sposób aby suma kwadratów różnic pomiędzy punktami na prostej, a odpowiadającymi im obserwacjami empirycznymi była jak najmniejsza.

Wielkość tej sumy (lub pierwiastek kwadratowy z sumy) jest miarą dokładności dopasowania (wariancja składnika losowego; albo średni błąd składnika losowego/resztowego – dla pierwiastka kwadratowego):

\[S_y = \sqrt{\frac{1}{n-2} \sum_{i=1}^n (y_t - \hat y_t)^2}\]

o ile średnio wartości empiryczne odchylają się od wartości teoretycznych wyznaczonych na podstawie funkcji trendu.

Udział wariancji składnika losowego w całości wariancji zmiennej \(Y\) jest inną miarą dopasowania (znaną jako współczynnik determinacji \(\Phi^2\):

\[\Phi^2 = \frac{\sum_{t=1}^n}{(y_t - y_t)^2}{\sum_{t=1}^n (y_t - \bar y_t)^2}\]

im ten udział jest mniejszy tym lepiej (lub jeżeli współczynnik zdefiniujemy jako 1 minus \(\Phi^2\) to im większy tym lepiej – współczynnik zbieżności czyli \(R^2\))

współczynnik zmienności resztowej:

\[V_{sy}=\frac{S_y}{\bar y_t}\cdot 100\]

Przykład Dane nt zgonów z powodu COVID w okresie 1.10.2020–5.2.2021 (linia niebieska trend dopasowany metodą Najmniejszych Kwadratów)

Współczynnik kierunkowy trendu liniowego wynoszący 1.5583885 jest interpretowany jako przeciętna zmiana z okresu na okres. Równanie prostej można zapisać jako:

zgony = 1.5583885 czas + 182.7417569

Interpretacja: w omawianym okresie przeciętnie umierało 1,5 osoby więcej dziennie. Ale dopasowanie linii prostej do danych jest słabe co widać oraz o czym świadczą wartości \(R^2\) (9.2613784%, tj. 9.2613784% zmienności jest objaśniane przez model) oraz średni błąd składnika losowego \(S_e\) (181.6602058). Ten błąd warto porówać do średniej wartości zmiennej objaśmnianej (liczby zgonów), która w omawianym okresie wynosi 283.2578125. Zatem błąd jaki popełniamy stanowi 64.1324609% średniej. Dużo (coś jakby średnio 1 \(\pm\) 0,6)

Inny przykład (nadwaga w PL):

Współczynnik kierunkowy trendu liniowego wynoszący 0.3502796 jest interpretowany jako przeciętna zmiana z okresu na okres. Równanie prostej można zapisać jako:

nadwaga = 0.3502796 czas + 9.8785134

Interpretacja: w omawianym okresie przeciętnie przybywało 0.3502796% osób z nadwagą rocznie. Dopasowanie linii prostej do danych jest bardzo dobre co widać oraz o czym świadczą wartości \(R^2\) (99.4245443%, tj. 99.4245443% zmienności jest objaśniane przez model) oraz średni błąd składnika losowego \(S_e\) (0.3309818). Ten błąd warto porówać do średniej wartości zmiennej objaśmnianej (liczby zgonów), która w omawianym okresie wynosi 17.4095238. Zatem błąd jaki popełniamy stanowi 1.9011537% średniej.