Laboratorium 8 -Regresja liniowa

1 Regresja

W analizie regresji celem jest wyjaśnienie lub modelowanie relacji pomiędzy zmienną \(Y=Y_1, \dots, Y_n\), nazywaną zmienną objaśnianą lub zmienna zależną, a jedną lub wieloma zmiennymi objaśniającymi (lub niezależnymi) \(X_1, \dots, X_p\).

Regresję można zapisać w następujący sposób: \[E(Y|X)=f(X,\beta).\] Interesującym nas zagadnieniem, będzie opisanie wartości oczekiwanej zmiennej objaśnianej \(Y\) za pomocą zmiennych objaśniających \(X_1, \dots, X_p\). A więc, celem jest wyznaczenie wektora parametrów \(\beta=\beta_1, \dots, \beta_p\), przy danej postaci modelu opisanego funkcją \(f\).

Regresja jest popularnym narzędziem statystycznym pozwalającym nie tylko na opisanie zależności pomiędzy zmiennymi objaśniającymi, a zmienną objaśnianą, ale również na oszacowanie średniej wartości zmiennej objaśnianej w zależności od zmiennych objaśniających oraz umożliwiającą interpretację zmiennych istotnie wpływających na zmienną objaśnianą.

# Regresja Liniowa

Model regresji w przypadku gdy funkcja \(f\) jest liniowa można zapisać w następujący sposób: \[y=X\beta+\epsilon,\] gdzie:

• \(y=y_1, \dots, y_n\) jest zmienną o długości \(n\) objaśnianą przez \(X\),
• \(X=\begin{bmatrix}1 &x_{11} & \dots & x_{1p}\\ 1 &x_{21} & \dots & x_{2p}\\ \vdots & &\vdots &\\ 1 & x_{n1}& \dots & x_{np} \end{bmatrix}\) jest macierzą (na ogół nie jest to macierz kwadratowa),
• \(\beta=(\beta_0, \dots, \beta_p)^t\) to wektor współczynników,
• \(\epsilon\) jest szumem o rozkładzie \(N(0, \sigma^2)\).

Model regresji liniowej opisuje oczekiwaną wartość średnią zmiennej \(y\) jako kombinację liniową zmiennych \(X\).

Celem jest estymacja nieznanych współczynników \(\beta\).

Założenia regresji liniowej

Określając model regresji liniowej musimy przyjąć poniższe założenia:

1. Postać modelu liniowego - zależność między zmienną zależną, a zmiennymi niezależnymi jest liniowa.

2. Obserwacje \(X\) są niezależne od siebie. Oznacza to, że macierz \(X\) jest rzędu p. Założenie to jest konieczne, aby istniało jednoznaczne rozwiązanie dla \(\hat{\beta}\).

3. Homoskedastyczność - wariancja błędów (zmiennych \(\epsilon\)) jest stała.

4. Zmienne \(\epsilon\) są niezależne i mają rozkład normalny o średniej \(0\).

Przykład: Rozważmy następujący model: \[y=f(X_1, X_2, X_3)+\epsilon.\] Zakładając, że funkcja \(f\) jest liniowa możemy go zapisać jako \[y=\beta_0+\beta_1 X_1 +\beta_2 X_2 +\beta_3 X_3+\epsilon,\] gdzie \(\beta_0\) jest stałą. Celem w tym przypadku jest estymacja czterech nieznanych współczynników \(\beta_0, \beta_1, \beta_2, \beta_3\). Dla ułatwienia zapisu posłużymy się zapisem macierzowym, tj: \[y=X\beta+\epsilon,\] gdzie \(y=(y_1, \dots, y_n)^t\), \(\epsilon=(\epsilon_1, \dots, \epsilon_n)\), \(\beta=(\beta_0, \dots, \beta_3)^t\) oraz \[X=\begin{bmatrix} 1 &x_{11} & x_{12} & x_{13}\\ 1 &x_{21} & x_{22} & x_{23}\\ \vdots & &\vdots &\\ 1 & x_{n1}& x_{n2}& x_{n3} \end{bmatrix}.\] Pierwsza kolumna jedynek w macierzy \(X\) odpowiada wyrazowi wolnemu \(\beta_0\).

Uwaga: do modelu liniowego zmiennymi objaśniającymi mogą być zarówno dane ilościowe jak i jakościowe, natomiast zmiennymi objaśnianymi powinny być tylko zmienne ilościowe.

1.1 Estymacja wektora współczynników \(\beta\)

Estymator wektora współczynników można uzyskać korzystając z następującego wzoru:

\[\hat{\beta}=\left(X^t X\right)^{-1}X^t y.\]

1.1.1 Wyprowadzenie

Naszym celem jest znalezienie takich współczynników wektora \(\beta\), aby \(X\beta\) było najbliżej ,jak to możliwe zmiennej objaśnianej \(y\).

Różnice pomiędzy estymowanymi wartościami przez model liniowy a prawdziwymi nazywać będziemy resztami lub residuami, tj \[\hat{\epsilon_i}=y_i-\hat{y_i},\] gdzie \(\hat{y_i}=X_i\beta^t.\)

Można to zrobić za pomocą metody najmniejszych kwadratów (Least Squares Estimation). Polega ona znalezieniu takich \(\beta\), które minimalizują sumę kwadratu błędu popełnianego przy estymacji, tj: \[SSE=\text{argmin}_{\beta}\sum_{i=1}^n\left(\epsilon_i\right)^2 =\text{argmin}_{\beta}\sum_{i}^n\left(\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\beta_3X_{3i}-y_i\right)^2.\]

W celu znalezienia takiego \(\beta\), które minimalizuje powyższą funkcję obliczmy pochodną ze względu na \(\beta_0, \beta_1, \beta_2, \beta_3\) i przyrównujemy do \(0\).

\[\frac{\partial SSE}{\partial \beta_0}=\frac{\partial \sum_{i=1}^n\left(\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\beta_3X_{3i}-y_i\right)^2}{\partial \beta_0}=\sum_{i=1}^n 2\left(\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\beta_3X_{3i}-y_i\right)\cdot \mathbf{1}\] \[\frac{\partial SSE}{\partial \beta_0}=\sum_{i=1}^n2X_{0i}^t\left(\hat{y}_i-y_i\right)=0\] Ze względu na \(\beta_1\) mamy \[\frac{\partial SSE}{\partial \beta_1}=\frac{\partial \sum_{i=1}^n\left(\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\beta_3X_{3i}-y_i\right)^2}{\partial \beta_1}=\sum_{i=1}^n 2\left(\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\beta_3X_{3i}-y_i\right)\cdot X_{1i}\] \[\frac{\partial SSE}{\partial \beta_0}=\sum_{i=1}^n2X_{1i}^t\left(\hat{y}_i-y_i\right)=0\] Postępując analogicznie dla pozostałych współczynników, otrzymujemy że dla dowolnego \(i =1, \dots, n\) zachodzi \[(\ast)\qquad2X_{\cdot i}^t\hat{\epsilon}_i=0\] Używając zapisu macierzowego mamy \[2X^t\left(X\beta-y\right)=0\] \[2X^tX\beta-2X^t y=0\] Ostatecznie, otrzymujemy wzór na estymator \(\hat{\beta}\):

\[\hat{\beta}=\left(X^tX\right)^{-1}X^t y.\] W praktyce stosuje się inną zależność o dużo lepszych właściwościach numerycznych, dajce te same rezultaty.

Uwaga: nie pokazaliśmy, że znaleziony estymator jest to ekstremum (a tym bardziej minimum) dla naszej funkcji. W tym celu musielibyśmy policzyć drugą pochodną (a w przypadku macierzy hesjan).

1.1.2 Uzasadnienie geometryczne

Można również pokazać, że \(\hat{\beta}\) jest minimum dla naszej funkcji podejściem geometrycznym. W tym celu, opiszemy kilka pojęć.

Niech macierz \(N\) rozmiaru \(n\times 1\) (inaczej wektor \(n\) -elementowy) należy do \(\mathbb{R}^n\). Dodatkowo macierz \(N\) należy do przestrzeni Euklidesowej o wymiarze \(n\), którą będziemy oznaczać przez \(E^n\). Dla przestrzeni \(E^n\) zdefiniujemy iloczyn skalarny jako \[\left\langle \mathbf{x,y}\right\rangle=\mathbf{x^ty}.\] Dodatkowo, dowolny zbiór wektorów \(n\)-elementowych \(x_1,\dots, x_k\) rozpościera podprzestrzeń \(k\) wymiarową przestrzeni \(E^n\). Dokładniej, będziemy to oznaczać przez \[\sigma(X)=\sigma(x_1, \dots, x_k)=\left\{z\in E^n|z=\sum_{i=1}^k b_ix_i, b_i\in \mathbb{R} \right\},\] gdzie wektory \(x_1,\dots, x_k\) nazywać będziemy wektorami bazowymi. Innymi słowy \(\sigma(X)\) jest powłoką liniową (span(\(X\))) zbioru \(X\).

Dopełnienie ortogonalne \(\sigma(X)\) oznaczymy przez \(\sigma(X)^{\perp}\). Dokładnie, \[\sigma(X)^{\perp}=\left\{\mathbf{w}\in E^n|\mathbf{w^t z=0}\text{ for every }z\in \sigma(X)\right\}.\]

Dla przykładu, rozważmy dwa dowolne wektory \(\mathbf{x_1, x_2}\) takie, że \(\mathbf{x_1\neq x_2}\) oraz wspólnym początku (patrzy rysunek poniżej). Wtedy \(\sigma(\mathbf{x_1, x_2})\) jest płaszczyzną zawierającą, wszystkie liniowe kombinacje wektorów \(\mathbf{x_1, x_2}\).

Używając tego wprowadzania możemy przedstawić geometrycznie metodę najmniejszych kwadratów (patrz poniży wykres).

W metodzie najmniejszych kwadratów minimalizujemy \(\epsilon^2\). Graficznie więc będzie to spełnione dla \(\epsilon\) ortogonalnego do \(span(X\beta)\) (patrzy rysunek poniżej).

Zatem aby pokazać, że \(\hat{\beta}\) jest ekstremum (minimum) dla naszej funkcji, wystarczy pokazać, że \[\left\langle X\beta, \hat{\epsilon}\right\rangle=0\] \[(X\beta)^t\hat{\epsilon}=0\] \[\beta^tX^t\hat{\epsilon}=0.\] Co zachodzi dla z warunku \((\ast)\).

1.1.3 Model liniowy w R

W modelu liniowym zakładamy liniowość parametrów, natomiast nie musimy mieć liniowości zmiennych objaśniających. Dla przykładu, \[y=\beta_0+\beta_1X_2+\beta_2 \log(X_2)+\epsilon\] jest modelem liniowym. Natomiast modelem liniowym nie jest model: \[y=\beta_0+\beta_1X_1^{\beta_2}+\epsilon.\]

Funkcją służąca do budowy modelu liniowego w języku R jest lm(). Model liniowy opisujemy za pomocą formuły, czyli symbolicznego opisu zależności pomiędzy zmiennymi. Składnia formuły jest następująca: Lewa.strona ~ Prawa.strona. Lewa.strona to najczęściej jedna zmienna (zmienna objaśniana) lub wyrażenie. Prawa.strona to jedna lub więcej zmiennych rozdzielonych najczęściej znakiem +. W celu konstrukcji bardziej skomplikowanych modeli odsyłamy do: tworzenie modelu w R.

Przyjrzyjmy się przykładowi, w którym modelujemy cenę mieszkania na podstawie jego powierzchni oraz liczby pokoi.

library("Przewodnik")

## Warning in check_dep_version(): ABI version mismatch: 
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 2
## Please re-install lme4 from source or restore original 'Matrix' package

model_lin <- lm(mieszkania$cena~mieszkania$powierzchnia+mieszkania$pokoi)

Funkcja lm() nie tylko wyznacza dopasowanie modelu, ale również wyznacza wartości współczynników \(\beta\) ich statystyczną istotność oraz wylicza wartości reszt. Aby usunąć z modelu wyraz wolny należy dodać do formuły -1. Wynikiem jest obiekt klasy lm przechowujący informacje o dopasowaniu modelu liniowego.

Aby z obiektu klasy lm wydobyć obliczony wektor współczynników należy użyć $coeff:

model_lin$coeff

##             (Intercept) mieszkania$powierzchnia        mieszkania$pokoi 
##              82407.0883               2070.8966               -840.1008

Zilustrujmy zależność pomiędzy ceną a powierzchnią za pomocą prostej: \(y=82407+2070x\)

plot(mieszkania$powierzchnia,mieszkania$cena, xlab="Powierzchnia", ylab="Cena")
abline(a=82407, b=2070, col="red")

Więcej informacji o modelu można uzyskać używając funkcji summary(). W wyniku jej wywołania, wypisywane są informacje o dopasowanym modelu, takie jak: wartości estymowanych współczynników, wyniki testów istotności oraz wartości \(R^2\). Wynikom testów istotności i wartości \(R^2\) poświecimy w dalszej części tego laboratorium więcej uwagi.

summary(model_lin)

## 
## Call:
## lm(formula = mieszkania$cena ~ mieszkania$powierzchnia + mieszkania$pokoi)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39705  -9386   -863   9454  35097 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              82407.1     2569.9  32.066   <2e-16 ***
## mieszkania$powierzchnia   2070.9      149.2  13.883   <2e-16 ***
## mieszkania$pokoi          -840.1     2765.1  -0.304    0.762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14110 on 197 degrees of freedom
## Multiple R-squared:  0.8937, Adjusted R-squared:  0.8926 
## F-statistic: 828.3 on 2 and 197 DF,  p-value: < 2.2e-16

W wyniku otrzymujemy następujące informacje:

• Call - formuła modelu, którą użyto do oszacowania,

• Residuals - reszty, czyli różnice między rzeczywistymi wartościami zmiennej zależnej (cena mieszkań) a wartościami przewidywanymi przez model. Wyświetlone mamy tutaj podstawowe statystyki (minumum, pierwszy kwantyl, medianę, trzeci kwantyl oraz wartość maksymalną) dotyczącą reszt dla naszego modelu,

• Coefficients - oszacowane współczynniki (Estimate) \(\beta\) modelu liniowego. Dodatkowo, znajdują się również następujące informacje informację:

  • błąd standardowy (std. error) - niskie wartości wskazują na precyzyjniejsze oszacowania.

  • t value - im wyższa wartość bezwzględna, tym bardziej prawdopodobne, że dana zmienna jest istotna.

  • Pr(>|t|) - p-wartość bardzo istotna kolumna która określa czy dana zmienna objaśniająca jest istotna statystycznie.

    • p<0.001 (***) — bardzo istotny współczynnik,

    • p<0.01p<0.01 (**) — istotny współczynnik,

    • p<0.05p<0.05 (*) — mniej istotny, ale akceptowalny,

    • p<0.1p<0.1 (.) — marginalnie istotny,

    • p>0.1p>0.1 (brak oznaczeń) — zmienna nieistotna.

    W naszym przypadku powierzchnia mieszkania ma bardzo niską p-wartość (<2e−16), co oznacza, że jest to bardzo istotna zmienna, natomiast liczba pokoi ma wysoką p-wartość (0.762), co sugeruje, że nie ma istotnego wpływu na cenę.

• Residual standard error - średni błąd prognozy. Mniejsza wartość resztowego błędu standardowego wskazuje, że model jest bardziej dopasowany do danych.

• F-statistic - jest obliczana przy przeprowadzaniu testu statystycznego sprawdzającego czy istnieje statystycznie istotna liniowa zależność między zmienną objaśnianą a zmiennymi objaśniającymi. Na jej podstawie obliczana jest p-wartość. Niska p-wartość (np. poniżej 0.05) świadczy o tym, że że istnieje istotna zależność między zmienną zależną a przynajmniej jedną zmienną objaśniającą. Z drugiej strony wysoka p-wartość sugeruje, że zmienne objaśniające nie mają istotnego wpływu na zmienną zależną (brak liniowej zależności).

1.2 Współczynnik \(R^2\)

Określa jak dobrze model opisuje dane. Wyraża się on następującym wzorem: \[R^2=1-\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{\sum_{i=1}^n(y_i-\bar{y})^2},\] gdzie \(\hat{y}_i=X_i\beta^{t}\) i \(\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\). Współczynnik \(R^2\) przyjmuje wartości z zakresu od 0 do 1. Im bliżej \(R^2\) jest wartości 1, tym lepsze jest dopasowanie modelu do danych.

Uwaga: współczynnik \(R^2\) ma pewne ograniczenia w szczególności, gdy model liniowy nie posiada wyrazu wolnego. Wtedy pojawia się problem przy obliczaniu mianownika, ponieważ całkowita suma kwadratów jest mierzona wokół zera, a nie średniej ruchomej. Skutkuje to tym, że \(R^2\) może osiągnąć wtedy wartości ujemne co wpływa na to, że model dopasowany jest do danych gorzej niż średnia.

Dlatego też w praktyce często zwraca się uwagę na tak zwany skorygowany \(R^2\) (Adjusted R-squared). Wyraża się go za pomocą poniższej formuły: \[R^2_{adj}=1-\frac{(1-R^2)(n-1)}{n-k-1},\] gdzie \(n\) to liczba obserwacji, a \(k\) to ilość parametrów w modelu (u nas \(p+1\), bo \(\beta_0, \dots,\beta_p\)). Również przyjmuje wartości z przedziału \([0,1]\) i podobnie jak zwykły \(R^2\) im bliżej osiąga wartości 1 tym lepiej. Dodatkowo, uwzględnia on liczbę zmiennych niezależnych w modelu oraz liczbę obserwacji. To znaczy, karze za dodawanie zmiennych niezależnych, które nie przyczyniają się znacząco do wyjaśnienia zmienności zmiennej zależnej. Dzięki temu może służyć do porównywania modeli.

Informacje na temat współczynników

\(R^2\) oraz \(R^2_{ajd}\) można znaleźć przy wyświetleniu summary(model_lin) lub wprost wartość można wyświetlić za pomocą komendy summary(model_lin)$r.squared lub summary(model_lin)$adj.r.squared.

1.2.1 Zadania - regresja liniowa

Załaduj dane gala z biblioteki faraway. Zbiór danych opisuje liczbę gatunków żółwi na różnych wyspach Galapagos. W zbiorze danych znajduje się 30 przypadków (wysp) i 7 zmiennych, takich jaki:

• Species - liczba gatunków żółwi występujących na wyspie,

• Endemics - liczba gatunków endemicznych (czyli żyjących tylko na danym obszarze),

• Elevation - najwyższe wzniesienie wyspy (w metrach),

• Nearest - odległość od najbliższej wyspy (w kilometrach),

• Scruz - odległość od wyspy Santa Cruz (w kilometrach),

• Adjacent - powierzchnia najbliższej wyspy (w kilometrach kwadratowych)

1. Używając funkcji lm() stwórz model liniowej zależności liczby gatunków (Species) od pozostałych zmiennych. Wyświetl szczegóły dotyczące tego modelu.

2. Oblicz samodzielnie estymator współczynników \(\hat{\beta}\) korzystając z podanego wzoru (otrzymane wartości powinny być takie same jak te uzyskane za pomocą funkcji lm()).

3. Oblicz samodzielnie współczynnik dopasowania \(R^2\) oraz \(R^2_{adj}\) korzystając z podanych wzorów (otrzymane wartości powinny być takie same jak te uzyskane za pomocą funkcji lm()). Wskazówka: reszty (residua) dla modelu można otrzymać za pomocą polecenia nazwa_modelu_liniowego$residuals.

4. Stwórz analogiczny model jak w podpunkcie a) (używając funkcji lm()), ale bez wyrazu wolnego. Porównaj otrzymany model z modelem z podpunktu a). Który model jest lepszy i dlaczego?

5. Stwórz model liniowy (używając funkcji lm()) opisujący liniową zależność liczby gatunków (Species) od zmiennych Elevation i Adjacent i również bez wyrazu wolnego. Porównaj go z modelami z podpunktów a) i d). Który model jest najlepszy i dlaczego?

1.3 Wykresy diagnostyczne dla modelu liniowego

Po dopasowaniu modelu, powinniśmy sprawdzić, czy spełnione są przyjęte założenia modelu liniowego. Jeżeli założenia są spełnione, to residua powinny mieć rozkład normalny o jednorodnej wariancji. Dlatego też, będziemy weryfikować założenia modelu badając właściwości reszt przy pomocy wykresów diagnostycznych. Wykorzystamy w tym celu przeciążoną funkcję plot() dla obiektów klasy lm . Otrzymamy wtedy cztery wykresy.

par(mfrow=c(2,2))
plot(model_lin)

1. Wykres o nagłówku “Residuals vs fitted” na osi x przestawia wartości dopasowane przez model \(\hat{y_i}\), a na osi y wartości reszt \(\hat{\epsilon_i}=y_i-\hat{y_i}\). Na podstawie tego wykresu możemy ocenić (czerwona krzywa), czy średnia zależy od wartości \(\hat{y_i}\) (to źle), czy jest bliska 0 (do dobrze). Dodatkowo, możemy zauważyć, czy wariancja jest jednorodna dla reszt, czy zmienia się w zależności od wartości \(\hat{y_i}\).

2. Wykres “Q-Q Residuals” (opisywany dokładniej w labolatorium 7) pozwala ocenić, czy reszty mają rozkład normalny. Przypomnijmy, że dla prawidłowego modelu punkty powinny układać się wzdłuż linii prostej \(y=x\).

3. Wykres “Scale-Location” na osi x przestawia wartości dopasowane przez model \(\hat{y_i}\), a na osi y pierwiastki zestandaryzowanych reszt. Podobnie, jak w przypadku wykresu pierwszego, tutaj również badamy jednorodność wariancji i liniowość modelu. Wykres spełniający założenia modelu liniowego powinien przedstawiać równomiernie rozmieszczone punkty wzdłuż pionowej czerwonej prostej.

4. Wykres “Residuals vs Leverage” jest używany do wykrywania wartości nietypowych. Na osi pionowej mamy zestandaryzowane reszty, a na osi poziomej tzw. siłę dźwigni \(h_i\), czyli miarę wpływu obserwacji na ocenę współczynników (ang. laverage). Jest ona wyznaczona według poniższego wzoru \[h_i=X_i(X^tX)^{-1}X_i^t.\] W języku R siłę dźwigni można obliczyć za pomocą funkcji hat() na macierzy \(X\), którą można wyodrębnić z modelu liniowego używając funkcji model.matrix(model_lin).

Siła dźwigni określa jaki wpływ na wartość współczynnika \(\beta\) ma obserwacja \(X_i\). W prawidłowym modelu liniowym pojedyncza obserwacja nie powinna mieć znacząco silniejszego wpływu na wartości współczynników niż pozostałe. Obserwacje można uznać za nietypową według tak zwanej reguły kciuka, gdy \(h_i\geq \frac{2p}{n}\), gdzie \(p\) jest liczbą parametrów (zmiennych objaśniających) w modelu liniowym.

1.3.1 Zadania - wykresy diagnostyczne

1. Dla najlepszego modelu z poprzedniego zadania wykonaj wykresy diagnostyczne i opisz co na ich podstawie można powiedz czy model spełnia założenia.

2. Znajdź nazwy wysp, które odpowiadają najmniejszej i największej wartości reszt.

3. Oblicz siłę dźwigni i nanieś uzyskane wartości na wykres. Następnie znajdź obserwacje nietypowe korzystając z reguły kciuka (zaznacz próg na wykresie poziomą prostą). Ile jest obserwacji nietypowych? Wypisz je. Spróbuj usunąć obserwacje nietypowe z zbioru danych, a następnie dopasować nowy model liniowy i dla niego wykonać testy diagnostyczne. Czy usuniecie zmiennych nietypowych pomogło ?

2 Testy statystyczne

Testowanie hipotez to kluczowy koncept w statystyce, który pozwala nam podejmować decyzje oparte na danych.

Pierwszym krokiem jest sformułowanie dwóch przeciwstawnych hipotez: hipotezę zerową ( \(H_0\)) i hipotezę alternatywną ( \(H_1\)).

Przed wykonaniem testu określa się również poziom istotności \(\alpha\). Jest to maksymalne prawdopodobieństwo odrzucenia prawdziwej hipotezy zerowej (czyli popełnienie błędu I rodzaju). Najczęściej przyjmujemy \(\alpha=0.05\), ale możemy również \(\alpha=0.01\) lub \(\alpha=0.001\).

Na podstawie posiadanego zestawu danych, wybieramy jaki test statystyczny chcemy przeprowadzić i obliczamy statystykę testową. Następnie, określamy rozkład statystyki testowej przy prawdziwości hipotezy zerowej i obliczamy p-wartość.

P-wartość jest to prawdopodobieństwo uzyskania statystyki testowej (przy prawdziwości hipotezy zerowej) co najmniej tak ekstremalnej jak zaobserwowana.

Jeśli p-wartość jest mniejsza niż założony poziom istotności, odrzucamy hipotezę zerową na korzyść hipotezy alternatywnej.

Jeśli wartość p jest większa niż poziom istotności, nie mamy podstaw do odrzucenia hipotezy zerowej.

Przy testowaniu hipotez możemy popełnić dwa rodzaje błędów:

• Błąd I rodzaju: hipoteza zerowa jest prawdziwa, ale odrzucamy ją na podstawie testu statystycznego.

• Błąd II rodzaju: hipoteza zerowa jest fałszywa, ale proces testowania stwierdza, że powinna zostać zaakceptowana.

Przyjmujemy, że błąd I rodzaju jest poważniejszy niż błąd II rodzaju i chcemy go kontrolować poprzez określenie poziomu istotności testu.

Przy testowaniu mówimy również o mocy testu. Określa on prawdopodobieństwo nie wystąpienia błędu II rodzaju, czyli prawdopodobieństwo odrzucenia hipotezy zerowej gdy jest ona fałszywa. Czyli jeżeli przez \(\beta\) określimy prawdopodobieństwo popełnienia błędu II rodzaju, to moc testu to jest \(1-\beta\).

W praktyce zatem poszukujemy testów na poziomie istotności \(\alpha\) i o maksymalnej mocy testu.

Możemy przeprowadzić trzy rodzaje testów statystycznych:

Rodzaje testów statystycznych

Przykład: Uczciwa moneta Rozważmy przykład testowania hipotezy, czy moneta, którą rzucamy, jest uczciwa (tzn. czy ma równe prawdopodobieństwo wypadnięcia orła i reszki)?.

Załóżmy, że rzuciliśmy kostką 30 razy i uzyskaliśmy 22 razy orła i 8 razy resztkę. Zaczęliśmy się zastanawiać, czy taki wynik jest statystycznie możliwy rzucając symetryczną monetą.

2.1 Test obustronny

Test dwustronny jest używany, gdy chcemy sprawdzić, czy badany parametr różni się od wartości oczekiwanej w dowolnym kierunku (czy jest większy lub mniejszy).

Sformułujmy naszą hipotezę \[H_0: p=1/2\qquad \text{vs}\qquad H_1: p\neq1/2,\] gdzie \(p\) oznacza prawdopodobieństwo wypadnięcia orła.

Przyjmijmy poziom istotności \(\alpha=0.05\).

Do przeprowadzenia testu użyjemy gotowej funkcji binom.test(), która służy do do testowania hipotezy dotyczącej prawdopodobieństwa sukcesu w próbie o rozkładzie dwumianowym. W szczególności jest ona przydatna, gdy mamy do czynienia z eksperymentem typu “sukces-porażka”, gdzie chcemy sprawdzić, czy obserwowany odsetek sukcesów różni się od zakładanego prawdopodobieństwa sukcesu.

k <- 22 # Liczba orłów
n <- 30 # Całkowita liczba prób

# Przeprowadzenie testu dwumianowego
test <- binom.test(k, n, p = 0.5, alternative = "two.sided")
test

## 
##  Exact binomial test
## 
## data:  k and n
## number of successes = 22, number of trials = 30, p-value = 0.01612
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5411063 0.8772052
## sample estimates:
## probability of success 
##              0.7333333

W wyniku otrzymujemy p-wartość mniejszą od założonego poziomu istotności \(0.05\), więc odrzucamy hipotezę zerową na rzecz alternatywy i stwierdzamy, że moneta może nie być uczciwa.

2.2 Test jednostronny

Test jednostronny stosujemy, gdy chcemy sprawdzić, czy badany parametr jest większy lub mniejszy od określonej wartości w jednym kierunku.

Sformułujmy teraz naszą hipotezę \[H_0: p\leq1/2\qquad \text{vs}\qquad H_1: p>1/2,\] gdzie \(p\) oznacza prawdopodobieństwo wypadnięcia orła.

test_jednostronny <- binom.test(k, n, p = 0.5, alternative = "greater")
test_jednostronny

## 
##  Exact binomial test
## 
## data:  k and n
## number of successes = 22, number of trials = 30, p-value = 0.008062
## alternative hypothesis: true probability of success is greater than 0.5
## 95 percent confidence interval:
##  0.5700661 1.0000000
## sample estimates:
## probability of success 
##              0.7333333

Zauważmy, ze otrzymana p-wartość jest ponownie mniejsza niż założony poziom istotności, zatem odrzucamy hipotezę zerową na rzecz alternatywy. Test zatem sugeruje, że prawdopodobieństwo uzyskania orła jest większe niż \(0.5\).

2.3 Istotność współczynników w modelu liniowym

Wyświetlając informacje o modelu liniowym dla zbioru danych mieszkania w ostatniej kolumnie otrzymujemy wyniki testów na istotność współczynników w modelu. Przypomnijmy:

summary(model_lin)

## 
## Call:
## lm(formula = mieszkania$cena ~ mieszkania$powierzchnia + mieszkania$pokoi)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39705  -9386   -863   9454  35097 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              82407.1     2569.9  32.066   <2e-16 ***
## mieszkania$powierzchnia   2070.9      149.2  13.883   <2e-16 ***
## mieszkania$pokoi          -840.1     2765.1  -0.304    0.762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14110 on 197 degrees of freedom
## Multiple R-squared:  0.8937, Adjusted R-squared:  0.8926 
## F-statistic: 828.3 on 2 and 197 DF,  p-value: < 2.2e-16

Wyświetlając summary(model_lin) w ostatniej kolumnie możemy zauważyć p-wartości, które wartości współczynników są istotnie większe od zera, a które nie. Dodatkowo, ilość gwiazdek również sygnalizuje jak bardzo dany współczynnik jest istotny dla naszego modelu.

Możemy zatem powiedzieć że stała (Intercept) oraz powierzchnia są współczynnikami istotnymi, ponieważ p-wartość jest mniejsza od założonego poziomy istotności.

Test statystyczny w tym przypadku jest postaci: \[H_0: \hat{\beta_i}=0\qquad \text{vs}\qquad H_1: \hat{\beta_i}\neq 0.\] Statystyka testowa dla tego testu, ma rozkład t-Studenta o \(n-p\) stopniach swobody. Można ją wyznaczyć korzystając z poniższego wzoru: \[t_i=\frac{\hat{\beta_i}}{se(\hat{\beta_i})},\] gdzie \(se(\hat{\beta_i})\) jest błędem standardowym. Dokładniej, błąd standardowy w tym przypadku wyraża się jako \[se(\hat{\beta_i})=\sqrt{\hat{\sigma}^2\left[(X^t X)^{-1}\right]_{ii}},\] gdzie \(\hat{\sigma}^2=\frac{1}{n-k-p-1}\epsilon^t\epsilon\) oraz \(\left[(X^t X)^{-1}\right]_{ii}\) oznacza elementy na diagonali. Następnie przeprowadzamy test obustronny. P-wartością jest prawdopodobieństwo osiągnięcia wartości równej statystyki testowej \(t\) lub bardziej ekstremalnej, tj \(2P(T\geq |t|)=2\left(1-P(T\leq |t|)\right)\). Jeżeli otrzymana wartość jest mniejsza od przyjętego poziomu istotności to odrzucamy hipotezę za rzecz alternatywy.

2.3.1 Zadanie - istotność współczynników

Sprawdź, które współczynniki są istotne statystycznie dla modelu liniowej zależności liczby gatunków (Species) od pozostałych zmiennych (z zadania - regresja liniowa, podpunkt a).

2.4 Testowanie założeń modelu liniowego

Przypomnijmy, że model liniowy opiera się na kilku kluczowych założeniach, które muszą być spełnione, aby wyniki estymacji były wiarygodne i pozwalały na poprawną interpretację.

Oprócz wykonania wykresów diagnostycznych bardzo często decydujemy się na wykonanie również testów statystycznych, które są mniej podatne na subiektywne interpretacje niż wykresy.

Przypomnijmy, że określając model regresji liniowej musimy przyjąć poniższe założenia:

1. Postać modelu liniowego - zależność między zmienną zależną, a zmiennymi niezależnymi jest liniowa.

2. Obserwacje \(X\) są niezależne od siebie. Oznacza to, że macierz \(X\) jest rzędu p. Założenie to jest konieczne, aby istniało jednoznaczne rozwiązanie dla \(\hat{\beta}\).

3. Homoskedastyczność - wariancja błędów (zmiennych \(\epsilon\)) jest stała.

4. Zmienne \(\epsilon\) są niezależne i mają rozkład normalny o średniej \(0\).

Testowanie postaci modelu liniowego

Przypomnijmy, że ten test wykonuje się przy budowaniu modelu liniowego za pomocą funkcji lm().

model_lin <- lm(mieszkania$cena~mieszkania$powierzchnia+mieszkania$pokoi)
summary(model_lin)

## 
## Call:
## lm(formula = mieszkania$cena ~ mieszkania$powierzchnia + mieszkania$pokoi)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39705  -9386   -863   9454  35097 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              82407.1     2569.9  32.066   <2e-16 ***
## mieszkania$powierzchnia   2070.9      149.2  13.883   <2e-16 ***
## mieszkania$pokoi          -840.1     2765.1  -0.304    0.762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14110 on 197 degrees of freedom
## Multiple R-squared:  0.8937, Adjusted R-squared:  0.8926 
## F-statistic: 828.3 on 2 and 197 DF,  p-value: < 2.2e-16

Tutaj p-wartość wynosi \(< 2.2e-16\), wiec odrzucamy hipotezę zerową (brak zależności) i stwierdzamy, że istnieje statystycznie istotna zależność liniowa między tymi zmiennymi.

Testowanie niezależności obserwacji X

Testowanie niezależności obserwacji sprowadza się do analizy reszt modelu. Niezależność reszt sugeruje, że błędy modelu nie są skorelowane ze sobą.

Testowanie homoskedastyczności

Test Breusch-Pagana służy sprawdzeniu, czy wariancja błędów jest stała. Wystarczy użyć funkcji bptest() z biblioteki lmtest.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(model_lin)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_lin
## BP = 0.5937, df = 2, p-value = 0.7432

Jeżeli p-wartość jest większa od przyjętego poziomu istotności, to nie odrzucamy hipotezy zerowej, że wariancja jest jednorodna.

Niezależność reszt

Przedstawimy, bez szczegółów matematycznych, test Durbin-Watsona, który bada czy występuje autokorelacja reszt. W tym celu użyjemy funkcji durbinWatsonTest() z biblioteki car.

library("car")

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

durbinWatsonTest(model_lin)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.1150583      1.764002   0.088
##  Alternative hypothesis: rho != 0

Test Durbin-Watsona testuje hipotezę zerową, że nie ma autokorelacji reszt (rho = 0), przeciwko alternatywie, że istnieje autokorelacja reszt (rho != 0). Zatem, wysoka p-wartość (większa od poziomu istotności 0.05) świadczy, że nie ma podstaw do odrzucenia hipotezy zerowej, wobec tego, możemy przyjąć, że obserwacje są niezależne.

Reszty mają rozkład normalny

Test Shapiro-Wilka jest jednym z najczęściej używanych testów normalności.

shapiro.test(model_lin$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_lin$residuals
## W = 0.99423, p-value = 0.6342

Jeżeli p-wartość jest większa od przyjętego poziomu istotności, to przyjmujemy hipotezę zerową, że reszty mają rozkład normalny. Natomiast, jeżeli p-wartość jest mała (mniejsza od przyjętego poziomu istotności) to odrzucamy hipotezę zerową na rzecz alternatywy i stwierdzamy, że reszty nie mają rozkładu normalnego.

2.5 Zadanie - testowanie załóżeń

Dla najlepszego modelu z zadania - regresja liniowa, przeprowadź odpowiednie testy statystyczne sprawdzając, czy są spełnione założenia modelu liniowego. Narysuj histogram dla reszt.

2.6 Automatyczny wybór zmiennych do modelu liniowego

Po co w ogóle się przejmować wyborem zmiennych?

Wymienimy tutaj kilka powodów dla których wybór “najlepszego” podzbioru predykatorów jest ważny:

1. Według zasady Brzytwy Ockhama najlepszym wyjaśnieniem danego zjawiska spośród kilku prawdopodobnych jest to najprostsze. W przypadku regresji oznacza to, że najlepszy jest najmniejszy model pasujący do danych.

2. Niepotrzebne zmienne zwiększają szum przy estymowaniu innych wielkości, np. predykcji.

3. Współliniowość (czyli sytuacja, gdy zmienne niezależne w modelu regresji są silnie skorelowane, co prowadzi do niestabilnych oszacowań współczynników) jest spowodowana zbyt dużą liczbą zmiennych.

4. Koszt: jeśli model ma być używany do przewidywania, możemy zaoszczędzić czas lub pieniądze.

Usuwanie i dodawanie zmiennych do modelu możemy wykonywać ręcznie lub automatycznie. Do automatycznego wybory istotnych zmiennych objaśniających służy funkcja step(). Stosuje ona metodę krokową zależną od wartości argumentu direction, która może być równa: forward, backward, both.

Przedstawimy działanie poszczególnych metod:

Eliminacja Forward

1. Rozpoczynamy od modelu z jednym składnikiem - wyrazem wolnym.

2. W każdym kroku dodajemy do modelu jedną zmienną uznaną za najbardziej istotną (z pozostałych zmiennych i ze względu na ustalone kryterium), tak długo, aż nie będzie już żadnych istotnych zmiennych do dodania.

Eliminacja Backward

1. Rozpoczynamy od pełnego modelu ze wszystkimi zmiennymi.

2. W każdym kroku usuwamy z modelu najmniej istotną zmienną, tak długo, aż wszystkie zmienne pozostałe w modelu będą istotne.

Eliminacja Forwadr-Backward

1. Rozpoczynamy od modelu z jednym składnikiem - wyrazem wolnym.

2. W każdym kroku dodajemy do modelu jedną zmienna z najmniejszą p-wartością, o ile jest ona istotna (w sensie poprawy danego kryterium).

3. Usuwamy z modelu zmienną z największą p-wartością, o ile jest ona nieistotna (w sensie danego kryterium).

4. Powtarzamy kroki 2 i 3 tak długo, aż model przestanie się zmieniać.

Funkcja step() wybiera zmienne ze względu na kryterium AIC, które chcemy minimalizować.

Przykładowe użycie funkcji step().

modelKoncowy <- step(lm(cena~., data=mieszkania), direction="forward")

## Start:  AIC=3679.97
## cena ~ pokoi + powierzchnia + dzielnica + typ.budynku

AIC(modelKoncowy) # Wyznacza wartość AIC dla ostatecznego modleu 

## [1] 4249.545

2.6.1 Zadanie - automatyczny wybór zmiennych

Korzystając z funkcji step(), wyznacz istotne zmienne dla zbiory danych gala korzystając w wszystkich przedstawionych metod eliminacji. Czy uzyskane w ten sposób modele zawierają ten sam zestaw zmiennych? Jeśli nie, to wybierz model, na podstawie najmniejszej wartości AIC oraz \(R^2_{adj}\).

3 Predykcja dla modelu liniowego

Mając dopasowany model liniowy możemy użyć obliczonych wartości współczynników \(\beta\) w celu wyznaczenia predykcji, czyli wyznaczenia oczekiwanej wartości zmiennej \(y\) dla zadanych wartości \(X\). Możemy to zrobić za pomocą funkcji predict().

Dla przykładu, zbudujmy model liniowy dla zbioru danych mieszkania z biblioteki Przewodnik, w którym cena zależy od powierzchni i dzielnicy. Następnie, za pomocą zbudowanego modelu określmy średnią cenę mieszkania o powierzchni 48 \(m^2\) w dzielnicy Biskupin i o powierzchni 68 \(m^2\) na Krzykach.

modelPD <- lm(cena~powierzchnia+dzielnica, data=mieszkania)
nowe_dane <- data.frame(powierzchnia=c(48,68), 
                        dzielnica=c("Biskupin", "Krzyki"))
predict(modelPD, newdata = nowe_dane)

##        1        2 
## 191415.9 210976.9

3.1 Zadanie - predykcja

Dla modelu z poprzedniego zadnia (zawierającego najlepszy zestaw zmiennych) przeprowadź predykcję dla dwóch nowych obserwacji (samodzielnie wymyślone).