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\). 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ą.
1.1 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:
Postać modelu liniowego - zależność między zmienną zależną, a zmiennymi niezależnymi jest liniowa.
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}\).
Homoskedastyczność - wariancja błędów (zmiennych \(\epsilon\)) jest stała.
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 &\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.2 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.2.1 Wyprowadzenie (dla odważnych)
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.2.2 Uzasadnienie geometryczne (dla odważnych)
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 \(\mathrm{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.2.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.
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:
## (Intercept) mieszkania$powierzchnia mieszkania$pokoi
## 82407.0883 2070.8966 -840.1008Zilustrujmy 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.
##
## 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-16W 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 kwartyl, medianę, trzeci kwartyl 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,
0.001 < p <0.01 (**) — istotny współczynnik,
0.01 < p <0.05 (*) — mniej istotny, ale akceptowalny,
0.05 < p <0.1 (.) — marginalnie istotny,
0.1 < p < 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ę (i spokojnie można było pominąć na wykresie).
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 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.4 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.
1.2.5 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.6 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)
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.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()).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ą polecenianazwa_modelu_liniowego$residuals.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?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.2.7 Wskazówki/przykładowe wyniki
## [1] "a"##
## Call:
## lm(formula = Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
## data = gala)
##
## Residuals:
## Min 1Q Median 3Q Max
## -111.679 -34.898 -7.862 33.460 182.584
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.068221 19.154198 0.369 0.715351
## Area -0.023938 0.022422 -1.068 0.296318
## Elevation 0.319465 0.053663 5.953 3.82e-06 ***
## Nearest 0.009144 1.054136 0.009 0.993151
## Scruz -0.240524 0.215402 -1.117 0.275208
## Adjacent -0.074805 0.017700 -4.226 0.000297 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.98 on 24 degrees of freedom
## Multiple R-squared: 0.7658, Adjusted R-squared: 0.7171
## F-statistic: 15.7 on 5 and 24 DF, p-value: 6.838e-07## [1] "b"## [,1]
## 1 7.068220709
## Area -0.023938338
## Elevation 0.319464761
## Nearest 0.009143961
## Scruz -0.240524230
## Adjacent -0.074804832## [1] "c"## [1] 0.7047635## [1] "d"##
## Call:
## lm(formula = Species ~ Area + Elevation + Nearest + Scruz + Adjacent -
## 1, data = gala)
##
## Residuals:
## Min 1Q Median 3Q Max
## -116.638 -31.142 -7.858 37.744 182.422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Area -0.02664 0.02082 -1.280 0.212373
## Elevation 0.33065 0.04351 7.600 5.9e-08 ***
## Nearest 0.02590 1.03480 0.025 0.980232
## Scruz -0.21359 0.19913 -1.073 0.293682
## Adjacent -0.07646 0.01682 -4.545 0.000121 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 59.91 on 25 degrees of freedom
## Multiple R-squared: 0.8502, Adjusted R-squared: 0.8202
## F-statistic: 28.38 on 5 and 25 DF, p-value: 1.515e-09## [1] "e"##
## Call:
## lm(formula = Species ~ Elevation + Adjacent - 1, data = gala)
##
## Residuals:
## Min 1Q Median 3Q Max
## -103.44 -33.22 -10.31 23.89 203.42
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Elevation 0.27849 0.02413 11.539 3.73e-12 ***
## Adjacent -0.06910 0.01505 -4.593 8.43e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 59.77 on 28 degrees of freedom
## Multiple R-squared: 0.833, Adjusted R-squared: 0.8211
## F-statistic: 69.83 on 2 and 28 DF, p-value: 1.312e-111.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.
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 (to 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}\).
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\).
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.
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 funkcjimodel.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
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.
Znajdź nazwy wysp, które odpowiadają najmniejszej i największej wartości reszt.
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 ?
1.3.2 Przykładowe wykresy - odpowiedzi
## [1] "a"
## [1] "b"## Fernandina SantaCruz
## 0.3898219 203.4226235## [1] "c"
## [1] 3## Species Endemics Area Elevation Nearest Scruz Adjacent
## Fernandina 93 35 634.49 1494 4.3 95.3 4669.32
## Isabela 347 89 4669.32 1707 0.7 28.1 634.49
## SanSalvador 237 81 572.33 906 0.2 19.8 4.89##
## Call:
## lm(formula = Species ~ Elevation + Adjacent - 1, data = gala_bez)
##
## Residuals:
## Min 1Q Median 3Q Max
## -141.57 -41.87 -17.27 20.57 164.17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Elevation 0.32391 0.03301 9.813 4.69e-10 ***
## Adjacent -0.04713 0.04837 -0.974 0.339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 58.05 on 25 degrees of freedom
## Multiple R-squared: 0.7964, Adjusted R-squared: 0.7801
## F-statistic: 48.89 on 2 and 25 DF, p-value: 2.29e-09
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:
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.7333333W 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.
##
## 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.7333333Zauważ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.2.1 Zadanie - test jednostronny
W powyższym teście sprawdzić wynik, gdy jako alternatywę przyjmiemy
less. Jaka jest p-wartość? Jak zinterpretować wynik?
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:
##
## 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-16Wyś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:
Postać modelu liniowego - zależność między zmienną zależną, a zmiennymi niezależnymi jest liniowa.
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}\).
Homoskedastyczność - wariancja błędów (zmiennych \(\epsilon\)) jest stała.
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().
##
## 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-16Tutaj 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.
##
## studentized Breusch-Pagan test
##
## data: model_lin
## BP = 0.5937, df = 2, p-value = 0.7432Jeż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.
## lag Autocorrelation D-W Statistic p-value
## 1 0.1150583 1.764002 0.076
## Alternative hypothesis: rho != 0Test 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-Wilk normality test
##
## data: model_lin$residuals
## W = 0.99423, p-value = 0.6342Jeż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:
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.
Niepotrzebne zmienne zwiększają szum przy estymowaniu innych wielkości, np. predykcji.
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.
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
Rozpoczynamy od modelu z jednym składnikiem - wyrazem wolnym.
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
Rozpoczynamy od pełnego modelu ze wszystkimi zmiennymi.
W każdym kroku usuwamy z modelu najmniej istotną zmienną, tak długo, aż wszystkie zmienne pozostałe w modelu będą istotne.
Eliminacja Forward-Backward
Rozpoczynamy od modelu z jednym składnikiem - wyrazem wolnym.
W każdym kroku dodajemy do modelu jedną zmienna z najmniejszą p-wartością, o ile jest ona istotna (w sensie poprawy danego kryterium).
Usuwamy z modelu zmienną z największą p-wartością, o ile jest ona nieistotna (w sensie danego kryterium).
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().
## Start: AIC=3679.97
## cena ~ pokoi + powierzchnia + dzielnica + typ.budynku## [1] 4249.5452.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.93.1 Zadanie - predykcja
Dla modelu z poprzedniego zadnia (zawierającego najlepszy zestaw zmiennych) przeprowadź predykcję dla dwóch nowych obserwacji (samodzielnie wymyślone).