1 Przykład
ilustrujący
Wygenerujemy dane do przykładu:
set.seed(123)
x <- rnorm(100, mean = 10, sd = 2) # zmienna objaśniająca
y <- 5 + 2 * x + rnorm(100, mean = 0, sd = 3) # zmienna zależna z szumem
dane <- data.frame(x = x, y = y)
Zbadajmy zależność pomiędzy zmiennymi x a y
wykorzystując wykres punktowy.
Na podstawie otrzymanego wykresu możemy przypuszczać, że między
zmienną x, a y występuje zależność liniowa.
Spróbujemy zatem dopasować do danych możliwie najlepiej prostą regresji,
czyli model Regresję liniową.
Zakładamy model postaci: \[\hat{y}=\beta_0+\beta_1x,\] gdzie \(\hat{y}\) jest dopasowaną (możliwie jak
najlepiej) krzywą prostą do naszych danych, \(\beta_0\) to wyraz wolny oraz \(\beta_1\) to współczynnik kierunkowy. Na
tym etapie nie będziemy jeszcze zagłębiać się w sposób wyznaczania tych
współczynników — użyjemy gotowej funkcji lm():
##
## Call:
## lm(formula = y ~ x, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7220 -2.0505 -0.2625 1.7419 9.8712
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4787 1.6579 3.304 0.00133 **
## x 1.9213 0.1603 11.984 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.912 on 98 degrees of freedom
## Multiple R-squared: 0.5944, Adjusted R-squared: 0.5903
## F-statistic: 143.6 on 1 and 98 DF, p-value: < 2.2e-16
W wyniku działania summary(model) interesują nas przede
wszystkim wartości znajdujące się w tabeli Coefficients, w
kolumnie Estimate:
pierwsza wartość to \(\beta_0\),
a druga to \(\beta_1\).
Można je również uzyskać za pomocą polecenia
coef(model). Oznacza to, że dopasowany model ma postać:
\[\hat{y}=\beta_0+\beta_1x = 2.8941+2.2320
x.\] Dodajmy teraz linię regresji do wykresu punktowego:
Pytania które należy sobie zadać to:
Jak dobrze dopasowana jest ta prosta do danych?
Na ile możemy ufać uzyskanym współczynnikom?
Czy różnice między obserwowanymi a przewidywanymi wartościami są
małe i rozłożone losowo?
W kolejnej sekcji przejdziemy do bardziej matematycznego
sformułowania problemu, uogólniając go na przypadek z więcej niż jedną
zmienną objaśniającą.
Obserwacja:
Otrzymaliśmy prostą postaci \[\hat{y}=\beta_0+\beta_1x = 2.8941+2.2320
x.\] podczas gdy dane wygenerowane zostały korzystając z
zależności \[
y=5+2x
\] z dodanym szumem. Różnica wydaje się istotna, stąd też
potrzeba dokładnej analizy statystycznej.
2 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ą.
3 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)^t\), \(\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.
3.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.\]
3.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).
3.1.1.1 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)\).
3.1.2 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.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.
##
## 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
(minimum, 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,
0.001<p<0.01 (**) — istotny współczynnik,
0.01<p<0.05 (*) — mniej istotny, ale akceptowalny,
0.05<p<0.1 (.) — marginalnie istotny,
p>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 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).
3.1.3 Współczynnik \(R^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, że 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_{adj}\) 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.
3.1.4 Zadanie - regresja
liniowa
Załaduj dane airquality które są dostępne w bazowym
pakiecie R. Zbiór ten zawiera dzienne pomiary jakości powietrza w Nowym
Jorku w 1973 roku. Dane obejmują 153 obserwacje (dni) i 6 zmiennych,
takich jak:
Ozone – poziom ozonu (ppb) – zmienna objaśniana,
Solar.R – promieniowanie słoneczne (lang),
Wind – prędkość wiatru (mph),
Temp – temperatura (°F),
Month – miesiąc (maj–wrzesień),
Day – dzień miesiąca.
Ponieważ ten zbiór danych zawiera brakujące dane, uzupełnij je
wywołując poniższy kod:
air_data <- airquality
library(Hmisc)
air_data$Ozone <- impute(air_data$Ozone)
air_data$Solar.R <- impute(air_data$Solar.R)
Używając funkcji lm() stwórz model liniowej zależności poziomu
ozonu (Ozone) od pozostałych zmiennych. Wyświetl szczegóły dotyczące
dopasowanego 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ą polecenia
nazwa_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ść zamiennej Ozone od zmiennych Wind, Temp i
również bez wyrazu wolnego. Porównaj go z modelami z podpunktów a) i d).
Który model jest najlepszy i dlaczego?
3.2 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 (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}\).
Wykres “Q-Q Residuals” (opisywany dokładniej w laboratorium 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ż poziomej czerwonej
linii.
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.
3.2.1 Zadanie - wykresy
diagnostyczne
Dla modelu utworzonego w punkcie e) z poprzedniego zadania
wykonaj wykresy diagnostyczne i opisz czy na ich podstawie można
powiedzieć, że model spełnia założenia.
Do zbioru danych air_data dodaj dwie wartości
odstające, tworząc nowy zbiór danych air_data_ext,
tj.:
outliers <- data.frame(
Ozone = c(300, 25),
Solar.R = c(50, 300),
Wind = c(50, 30.0),
Temp = c(85, 40),
Month = c(10, 10),
Day = c(1, 2)
)
air_data_ext <- rbind(air_data, outliers)
Zbuduj dla tego zbioru danych model liniowy, w którym zmienna Ozone
zależy w sposób liniowy od zmiennych Wind, Temp i bez wyrazu wolnego.
Sprawdź jak w tym przypadku wyglądają wykresy diagnostyczne.
- Oblicz siłę dźwigni, a następnie narysuj ją na wykresie. Następnie
znajdź obserwacje nietypowe korzystając z reguły kciuka (zaznacz próg na
utworzonym wykresie poziomą prostą). Ile jest obserwacji nietypowych?
Wypisz je. Jak wpłynęło dodanie wartości odstających na model liniowy
(porównaj z podpunktem a)?
Przykładowe rozwiązania
##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_data_ext)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.631 -19.586 -11.372* 2.855 247.327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind 0.15093 0.47915 0.315 0.753
## Temp 0.53090 0.07029 7.553 3.61e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.38 on 153 degrees of freedom
## Multiple R-squared: 0.6281, Adjusted R-squared: 0.6232
## F-statistic: 129.2 on 2 and 153 DF, p-value: < 2.2e-16
## [1] 2
## Ozone Solar.R Wind Temp Month Day
## 154 300 50 50 85 10 1
## 155 25 300 30 40 10 2
##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_bez)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.819 -15.320* -4.617* 11.443 102.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind -3.48007 0.42559 -8.177 1.1e-13 ***
## Temp 0.95718 0.05728 16.711 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.05 on 151 degrees of freedom
## Multiple R-squared: 0.8003, Adjusted R-squared: 0.7977
## F-statistic: 302.6 on 2 and 151 DF, p-value: < 2.2e-16
4 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 monetą 30 razy i uzyskaliśmy 22 razy orła i 8
razy reszkę. Zaczęliśmy się zastanawiać, czy taki wynik jest
statystycznie możliwy rzucając symetryczną monetą.
4.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\colon \
p=1/2\qquad \text{vs}\qquad H_1\colon \ p\neq 1/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 być
nieuczciwa.
4.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
\colon \ p\leq1/2\qquad \text{vs}\qquad H_1 \colon \ 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.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\).
4.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-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
różne 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\colon \ \hat{\beta_i}=0\qquad \text{vs}\qquad
H_1\colon \ \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.
4.3.1 Zadanie - istotność
współczynników
Sprawdź (wyświetlając podsumowanie), które współczynniki są istotne
statystycznie dla modelu liniowej zależności zmiennej Ozone od
pozostałych zmiennych (z zadania - regresja liniowa, podpunkt a).
4.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-16
Tutaj p-wartość wynosi \(\verb+<
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.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.
## lag Autocorrelation D-W Statistic p-value
## 1 0.1150583 1.764002 0.084
## 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-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.
4.4.1 Zadanie -
testowanie załóżeń
Dla modelu z podpunktu e), z zadania - regresja liniowa, przeprowadź
odpowiednie testy statystyczne sprawdzając, czy są spełnione założenia
modelu liniowego. Narysuj histogram dla reszt.
4.5 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.545
4.5.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}\).
5 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
5.1 Zadanie -
predykcja
Dla modelu z poprzedniego zadnia (zawierającego najlepszy zestaw
zmiennych) przeprowadź predykcję dla dwóch nowych obserwacji
(samodzielnie wybrane).
1 Przykład ilustrujący
Wygenerujemy dane do przykładu:
set.seed(123)
x <- rnorm(100, mean = 10, sd = 2) # zmienna objaśniająca
y <- 5 + 2 * x + rnorm(100, mean = 0, sd = 3) # zmienna zależna z szumem
dane <- data.frame(x = x, y = y)Zbadajmy zależność pomiędzy zmiennymi x a y
wykorzystując wykres punktowy.
Na podstawie otrzymanego wykresu możemy przypuszczać, że między
zmienną x, a y występuje zależność liniowa.
Spróbujemy zatem dopasować do danych możliwie najlepiej prostą regresji,
czyli model Regresję liniową.
Zakładamy model postaci: \[\hat{y}=\beta_0+\beta_1x,\] gdzie \(\hat{y}\) jest dopasowaną (możliwie jak
najlepiej) krzywą prostą do naszych danych, \(\beta_0\) to wyraz wolny oraz \(\beta_1\) to współczynnik kierunkowy. Na
tym etapie nie będziemy jeszcze zagłębiać się w sposób wyznaczania tych
współczynników — użyjemy gotowej funkcji lm():
##
## Call:
## lm(formula = y ~ x, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7220 -2.0505 -0.2625 1.7419 9.8712
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4787 1.6579 3.304 0.00133 **
## x 1.9213 0.1603 11.984 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.912 on 98 degrees of freedom
## Multiple R-squared: 0.5944, Adjusted R-squared: 0.5903
## F-statistic: 143.6 on 1 and 98 DF, p-value: < 2.2e-16W wyniku działania summary(model) interesują nas przede
wszystkim wartości znajdujące się w tabeli Coefficients, w
kolumnie Estimate:
pierwsza wartość to \(\beta_0\),
a druga to \(\beta_1\).
Można je również uzyskać za pomocą polecenia
coef(model). Oznacza to, że dopasowany model ma postać:
\[\hat{y}=\beta_0+\beta_1x = 2.8941+2.2320
x.\] Dodajmy teraz linię regresji do wykresu punktowego:
Pytania które należy sobie zadać to:
Jak dobrze dopasowana jest ta prosta do danych?
Na ile możemy ufać uzyskanym współczynnikom?
Czy różnice między obserwowanymi a przewidywanymi wartościami są małe i rozłożone losowo?
W kolejnej sekcji przejdziemy do bardziej matematycznego sformułowania problemu, uogólniając go na przypadek z więcej niż jedną zmienną objaśniającą.
Obserwacja:
Otrzymaliśmy prostą postaci \[\hat{y}=\beta_0+\beta_1x = 2.8941+2.2320 x.\] podczas gdy dane wygenerowane zostały korzystając z zależności \[ y=5+2x \] z dodanym szumem. Różnica wydaje się istotna, stąd też potrzeba dokładnej analizy statystycznej.2 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ą.
3 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)^t\), \(\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.
3.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.\]
3.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).
3.1.1.1 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)\).
3.1.2 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.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.
##
## 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
(minimum, 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,
0.001<p<0.01 (**) — istotny współczynnik,
0.01<p<0.05 (*) — mniej istotny, ale akceptowalny,
0.05<p<0.1 (.) — marginalnie istotny,
p>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 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).
3.1.3 Współczynnik \(R^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, że 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_{adj}\) 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.
3.1.4 Zadanie - regresja
liniowa
Załaduj dane airquality które są dostępne w bazowym
pakiecie R. Zbiór ten zawiera dzienne pomiary jakości powietrza w Nowym
Jorku w 1973 roku. Dane obejmują 153 obserwacje (dni) i 6 zmiennych,
takich jak:
Ozone – poziom ozonu (ppb) – zmienna objaśniana,
Solar.R – promieniowanie słoneczne (lang),
Wind – prędkość wiatru (mph),
Temp – temperatura (°F),
Month – miesiąc (maj–wrzesień),
Day – dzień miesiąca.
Ponieważ ten zbiór danych zawiera brakujące dane, uzupełnij je
wywołując poniższy kod:
air_data <- airquality
library(Hmisc)
air_data$Ozone <- impute(air_data$Ozone)
air_data$Solar.R <- impute(air_data$Solar.R)
Używając funkcji lm() stwórz model liniowej zależności poziomu
ozonu (Ozone) od pozostałych zmiennych. Wyświetl szczegóły dotyczące
dopasowanego 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ą polecenia
nazwa_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ść zamiennej Ozone od zmiennych Wind, Temp i
również bez wyrazu wolnego. Porównaj go z modelami z podpunktów a) i d).
Który model jest najlepszy i dlaczego?
3.2 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 (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}\).
Wykres “Q-Q Residuals” (opisywany dokładniej w laboratorium 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ż poziomej czerwonej
linii.
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.
3.2.1 Zadanie - wykresy
diagnostyczne
Dla modelu utworzonego w punkcie e) z poprzedniego zadania
wykonaj wykresy diagnostyczne i opisz czy na ich podstawie można
powiedzieć, że model spełnia założenia.
Do zbioru danych air_data dodaj dwie wartości
odstające, tworząc nowy zbiór danych air_data_ext,
tj.:
outliers <- data.frame(
Ozone = c(300, 25),
Solar.R = c(50, 300),
Wind = c(50, 30.0),
Temp = c(85, 40),
Month = c(10, 10),
Day = c(1, 2)
)
air_data_ext <- rbind(air_data, outliers)
Zbuduj dla tego zbioru danych model liniowy, w którym zmienna Ozone
zależy w sposób liniowy od zmiennych Wind, Temp i bez wyrazu wolnego.
Sprawdź jak w tym przypadku wyglądają wykresy diagnostyczne.
- Oblicz siłę dźwigni, a następnie narysuj ją na wykresie. Następnie
znajdź obserwacje nietypowe korzystając z reguły kciuka (zaznacz próg na
utworzonym wykresie poziomą prostą). Ile jest obserwacji nietypowych?
Wypisz je. Jak wpłynęło dodanie wartości odstających na model liniowy
(porównaj z podpunktem a)?
Przykładowe rozwiązania
##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_data_ext)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.631 -19.586 -11.372* 2.855 247.327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind 0.15093 0.47915 0.315 0.753
## Temp 0.53090 0.07029 7.553 3.61e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.38 on 153 degrees of freedom
## Multiple R-squared: 0.6281, Adjusted R-squared: 0.6232
## F-statistic: 129.2 on 2 and 153 DF, p-value: < 2.2e-16
## [1] 2
## Ozone Solar.R Wind Temp Month Day
## 154 300 50 50 85 10 1
## 155 25 300 30 40 10 2
##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_bez)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.819 -15.320* -4.617* 11.443 102.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind -3.48007 0.42559 -8.177 1.1e-13 ***
## Temp 0.95718 0.05728 16.711 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.05 on 151 degrees of freedom
## Multiple R-squared: 0.8003, Adjusted R-squared: 0.7977
## F-statistic: 302.6 on 2 and 151 DF, p-value: < 2.2e-16
3.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.\]
3.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).
3.1.1.1 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)\).
3.1.2 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.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.
##
## 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
(minimum, 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,
0.001<p<0.01 (**) — istotny współczynnik,
0.01<p<0.05 (*) — mniej istotny, ale akceptowalny,
0.05<p<0.1 (.) — marginalnie istotny,
p>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 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).
3.1.3 Współczynnik \(R^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, że 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_{adj}\) 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.
3.1.4 Zadanie - regresja
liniowa
Załaduj dane airquality które są dostępne w bazowym
pakiecie R. Zbiór ten zawiera dzienne pomiary jakości powietrza w Nowym
Jorku w 1973 roku. Dane obejmują 153 obserwacje (dni) i 6 zmiennych,
takich jak:
Ozone – poziom ozonu (ppb) – zmienna objaśniana,
Solar.R – promieniowanie słoneczne (lang),
Wind – prędkość wiatru (mph),
Temp – temperatura (°F),
Month – miesiąc (maj–wrzesień),
Day – dzień miesiąca.
Ponieważ ten zbiór danych zawiera brakujące dane, uzupełnij je
wywołując poniższy kod:
air_data <- airquality
library(Hmisc)
air_data$Ozone <- impute(air_data$Ozone)
air_data$Solar.R <- impute(air_data$Solar.R)
Używając funkcji lm() stwórz model liniowej zależności poziomu
ozonu (Ozone) od pozostałych zmiennych. Wyświetl szczegóły dotyczące
dopasowanego 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ą polecenia
nazwa_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ść zamiennej Ozone od zmiennych Wind, Temp i
również bez wyrazu wolnego. Porównaj go z modelami z podpunktów a) i d).
Który model jest najlepszy i dlaczego?
3.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).
3.1.1.1 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)\).
3.1.2 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 (minimum, 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,
0.001<p<0.01 (**) — istotny współczynnik,
0.01<p<0.05 (*) — mniej istotny, ale akceptowalny,
0.05<p<0.1 (.) — marginalnie istotny,
p>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 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).
3.1.3 Współczynnik \(R^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, że 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_{adj}\) 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.
3.1.4 Zadanie - regresja liniowa
Załaduj dane airquality które są dostępne w bazowym
pakiecie R. Zbiór ten zawiera dzienne pomiary jakości powietrza w Nowym
Jorku w 1973 roku. Dane obejmują 153 obserwacje (dni) i 6 zmiennych,
takich jak:
Ozone – poziom ozonu (ppb) – zmienna objaśniana,
Solar.R – promieniowanie słoneczne (lang),
Wind – prędkość wiatru (mph),
Temp – temperatura (°F),
Month – miesiąc (maj–wrzesień),
Day – dzień miesiąca.
Ponieważ ten zbiór danych zawiera brakujące dane, uzupełnij je wywołując poniższy kod:
air_data <- airquality
library(Hmisc)
air_data$Ozone <- impute(air_data$Ozone)
air_data$Solar.R <- impute(air_data$Solar.R)Używając funkcji lm() stwórz model liniowej zależności poziomu ozonu (Ozone) od pozostałych zmiennych. Wyświetl szczegóły dotyczące dopasowanego 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ść zamiennej Ozone od zmiennych Wind, Temp i również bez wyrazu wolnego. Porównaj go z modelami z podpunktów a) i d). Który model jest najlepszy i dlaczego?
3.2 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 (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}\).
Wykres “Q-Q Residuals” (opisywany dokładniej w laboratorium 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ż poziomej czerwonej linii.
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.
3.2.1 Zadanie - wykresy diagnostyczne
Dla modelu utworzonego w punkcie e) z poprzedniego zadania wykonaj wykresy diagnostyczne i opisz czy na ich podstawie można powiedzieć, że model spełnia założenia.
Do zbioru danych
air_datadodaj dwie wartości odstające, tworząc nowy zbiór danychair_data_ext, tj.:
outliers <- data.frame(
Ozone = c(300, 25),
Solar.R = c(50, 300),
Wind = c(50, 30.0),
Temp = c(85, 40),
Month = c(10, 10),
Day = c(1, 2)
)
air_data_ext <- rbind(air_data, outliers)Zbuduj dla tego zbioru danych model liniowy, w którym zmienna Ozone zależy w sposób liniowy od zmiennych Wind, Temp i bez wyrazu wolnego. Sprawdź jak w tym przypadku wyglądają wykresy diagnostyczne.
- Oblicz siłę dźwigni, a następnie narysuj ją na wykresie. Następnie znajdź obserwacje nietypowe korzystając z reguły kciuka (zaznacz próg na utworzonym wykresie poziomą prostą). Ile jest obserwacji nietypowych? Wypisz je. Jak wpłynęło dodanie wartości odstających na model liniowy (porównaj z podpunktem a)?
Przykładowe rozwiązania
##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_data_ext)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.631 -19.586 -11.372* 2.855 247.327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind 0.15093 0.47915 0.315 0.753
## Temp 0.53090 0.07029 7.553 3.61e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.38 on 153 degrees of freedom
## Multiple R-squared: 0.6281, Adjusted R-squared: 0.6232
## F-statistic: 129.2 on 2 and 153 DF, p-value: < 2.2e-16
## [1] 2## Ozone Solar.R Wind Temp Month Day
## 154 300 50 50 85 10 1
## 155 25 300 30 40 10 2##
## Call:
## lm(formula = Ozone ~ Wind + Temp - 1, data = air_bez)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.819 -15.320* -4.617* 11.443 102.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Wind -3.48007 0.42559 -8.177 1.1e-13 ***
## Temp 0.95718 0.05728 16.711 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.05 on 151 degrees of freedom
## Multiple R-squared: 0.8003, Adjusted R-squared: 0.7977
## F-statistic: 302.6 on 2 and 151 DF, p-value: < 2.2e-164 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 monetą 30 razy i uzyskaliśmy 22 razy orła i 8 razy reszkę. Zaczęliśmy się zastanawiać, czy taki wynik jest statystycznie możliwy rzucając symetryczną monetą.
4.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\colon \
p=1/2\qquad \text{vs}\qquad H_1\colon \ p\neq 1/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 być
nieuczciwa.
4.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
\colon \ p\leq1/2\qquad \text{vs}\qquad H_1 \colon \ 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.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\).
4.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-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
różne 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\colon \ \hat{\beta_i}=0\qquad \text{vs}\qquad
H_1\colon \ \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.
4.3.1 Zadanie - istotność
współczynników
Sprawdź (wyświetlając podsumowanie), które współczynniki są istotne
statystycznie dla modelu liniowej zależności zmiennej Ozone od
pozostałych zmiennych (z zadania - regresja liniowa, podpunkt a).
4.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-16
Tutaj p-wartość wynosi \(\verb+<
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.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.
## lag Autocorrelation D-W Statistic p-value
## 1 0.1150583 1.764002 0.084
## 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-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.
4.4.1 Zadanie -
testowanie załóżeń
Dla modelu z podpunktu e), z zadania - regresja liniowa, przeprowadź
odpowiednie testy statystyczne sprawdzając, czy są spełnione założenia
modelu liniowego. Narysuj histogram dla reszt.
4.5 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.545
4.5.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}\).
4.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\colon \ p=1/2\qquad \text{vs}\qquad H_1\colon \ p\neq 1/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 być nieuczciwa.
4.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 \colon \ p\leq1/2\qquad \text{vs}\qquad H_1 \colon \ 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\).
4.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
różne 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\colon \ \hat{\beta_i}=0\qquad \text{vs}\qquad H_1\colon \ \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.
4.3.1 Zadanie - istotność współczynników
Sprawdź (wyświetlając podsumowanie), które współczynniki są istotne statystycznie dla modelu liniowej zależności zmiennej Ozone od pozostałych zmiennych (z zadania - regresja liniowa, podpunkt a).
4.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 \(\verb+< 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.084
## 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.
4.4.1 Zadanie - testowanie załóżeń
Dla modelu z podpunktu e), z zadania - regresja liniowa, przeprowadź odpowiednie testy statystyczne sprawdzając, czy są spełnione założenia modelu liniowego. Narysuj histogram dla reszt.
4.5 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.5454.5.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}\).
5 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.95.1 Zadanie - predykcja
Dla modelu z poprzedniego zadnia (zawierającego najlepszy zestaw zmiennych) przeprowadź predykcję dla dwóch nowych obserwacji (samodzielnie wybrane).