Analiza regresji

Analiza regresji

Konspekt do zaj¦¢: Statystyczne metody analizy danych

Agnieszka Nowak-Brzezi«ska 28 pa¹dziernika 2009

1 Opis zaj¦¢

Celem zaj¦¢ jest realizacja praktyczna zagadnie« zwi¡zanych z analiz¡ regresji, wykresami rezyduów oraz obliczeniem warto±ci korelacji, omówionych na wykªadzie Pana Profesora Jacka Koronackiego ( http://www.ipipan.waw.pl/~korona/).

2 Wprowadzenie

Analiza regresji jest bardzo popularn¡ i ch¦tnie stosowan¡ technik¡ statystyczn¡

pozwalajac¡ opisywa¢ zwi¡zki zachodz¡ce pomi¦dzy zmiennymi wej±ciowymi (obja±niaj¡cymi) a wyj±ciowymi (obja±nianymi). Innymi sªowy dokonujemy es- tymacji jednych danych, korzystaj¡c z innych. Istnieje wiele ró»nych technik regresji. Niew¡tpliwie najpopularniejsza jest regresja liniowa. Zakªada ona, »e pomi¦dzy zmiennymi obja±niaj¡cymi i obja±nianymi istnieje mniej lub bardziej wyrazista zale»no±¢ liniowa.

Procedura liniowej regresji w R

Szukamy równania y = bˆ 0+ b1x.Wiadomym jest, »e mi¦dzy danymi, które s¡

rzeczywiste yi (które analizujemy) a tych, które by±my chcieli uzyska¢ (oczeki- wanymi) - ˆyi - istnieje zwykle pewna rozbie»no±¢. Wyra»amy j¡ wzorem ei = yi− ˆyi(rezydua, reszty).

Precyzyjniej powiemy, »e gdy po wykonaniu wykresu rozrzutu obserwujemy,

»e chmura punktów (xi, yi) ukªada si¦ wzdªu» prostej, mo»emy spróbowa¢

wyznaczy¢ jej równanie. Szukamy wtedy modelu regresji dla próbki i staramy si¦ tak wyznaczy¢ wspóªczynniki b1 i b0 w ukªadzie równo±ci

yi= b1xi+ b0+ ei, dla i = 1, ...N,

by suma warto±ci bezwzgl¦dnych bª¦dów ei byªa jak najmniejsza. Mówimy wt- edy to tzw. metodzie najmniejszych kwadratów (MNK). Metoda MNK polega na tym, »e szukamy takich warto±ci b0oraz b1które minimalizuj¡ sum¦ kwadratów rezyduów:

XN i=1

(yi− ˆyi)² Dla takich danych:

b1= sxy

s²_x = P_N

i=1(xi− xi)(yi− yi) P_N

i=1(xi− xi)²

za±

b0= ˆy − b1xˆ

W ±rodowisku R ta procedura wykonuje si¦ w 3 krokach:

1. wyrysowanie punktów na wykresie rozrzutu: > plot(x,y) 2. znalezienie warto±ci b1 i b0

3. dodanie linii do grafu: abline(lm(y ~ x))

Funkcja abline jest do±¢ trudna w zrozumieniu. Drukuje ona linie na bie»¡- cym oknie gracznym. Wykorzystuje do tego znan¡ z modelu liniowego funkcj¦

lm. Wyra»enie y ∼ x mówi, »e chcemy znale¹¢ model dla zmiennej y w za- le»no±ci od zmiennej x. Funkcja cor pomo»e obliczy¢ wspóªczynnik korelacji Pearsona:

R =

P(xi− ¯x)(yi− ¯y) pP(xi− ¯x)²(yi− ¯y)²,

,który pozwala stwierdzie¢ jak jedna zmienna zale»y od zmian na drugiej zmi- ennej. Warto±ci R² bliskie 1 ±wadcz¡ o silnej korelacji, za± bliskie 0 o braku korelacji. W ±rodowisku R u»yjemy do tego funkcji cor

> cor(x,y) # to find R [1] 0.881

> cor(x,y)^2 # to find R^2 [1] 0.776

Musimy jawnie wywoªa¢ summary(lm(y ∼ x)).

Gdy chcemy wyznacza¢ warto±ci zmiennej x w zale»no±ci od y to zamieni- amy we wzorach x z y. Otrzymane proste pokrywaj¡ si¦, gdy badana zale»no±ci jest zale»no±ci¡ funkcyjn¡ i nie ma losowo±ci.

Przypadki odstaj¡ce i znacz¡ce: W celu wykluczenia z analizy przypadków odstaj¡cych, które mog¡ na ni¡ niekorzystnie wpªyn¡¢, nale»y zrobi¢ wykresy skrzynkowe analizowanych zmiennych. Na wykresach tych kóªkiem i gwiazdk¡

zaznaczone s¡ przypadki odstaj¡ce, odpowiednio nietypowe i skrajne. Przypadki te sugeruje si¦ usuwa¢, a w przypadku du»ej ich liczby analizowa¢ osobno. Tu- taj nale»y uwzgl¦dni¢ uwagi Pana Profesora na temat zmiennych znacz¡cych i odstaj¡cych, pami¦taj¡c, »e nie ka»da zmienna odstaj¡ca jest nie znacz¡ca, i »e czasami zmienna znacz¡ca mo»e by¢ równie» odstaj¡ca i wtedy sugeruje sie jej nie wyrzuca¢ z analizy, a wr¦cz j¡ uwzgl¦dnia¢ po to by minimalizowa¢ zadanie MNK

Wspóªczynnik korelacji liniowej Pearsona (Karl Pearson (1896)) rxy=

P(xi− ¯x)(yi− ¯y) (n − 1)sxsy

,

.Przyjmuje warto±ci z przedziaªu [−1, 1]. Dodatnia warto±¢ tego wspóªczynnika oznacza, »e wzrost warto±ci jednej zmiennej generalnie poci¡ga za sob¡ wzrost warto±ci drugiej zmiennej; ujemna  spadek. r = 0, gdy nie ma zwi¡zku mi¦dzy zmiennymi, krk ≈ 1, gdy zwi¡zek jest bardzo silny.

Zwi¡zek regresji i wspóªczynnika Pearsona

Im bli»szy 1, tym dopasowanie lepsze.

Interpretacja rxy² (tzw. wspóªczynnik determinacji):

r_xy² = 1 −s²_y|x s²_y ,

gdzie s²_y|x jest bª¦dem kwadratowym regresji liniowej dla xi na yi danej rów- naniem y = a + bx. Wi¦c jest to:

s²_y|x= 1 n − 1

Xn

i=1

(yi− a − bxi)²,

i s²y jest wariancj¡ dla y:

s²_y = 1 n − 1

Xn i=1

(yi− ¯y)².

Je±li korelacja b¦dzie symetryczna dla xii yiotrzymamy te same warto±ci dopa- sowania dla zmiennych yi i xi:

r_xy² = 1 −s²_x|y s²_x .

3 Wykonajmy analiz¦...

1. krok 1: wczytujemy dane

> wiek<-c(18,19,20,21,22,23,24,25,26,27,28,29)

> wzrost<-c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8,82.8,83.5) 2. Chc¡c znale¹¢ równanie regresji dla zmiennych wzrost oraz wiek, gdzie

zakªadamy, »e zmienn¡ obja±nian¡ (y) jest wzrost a obja±niaj¡c¡ (x) jest wiek, powiemy, »e interesuje nas równanie:

wzrost = b0+ b1∗ wiek

Oczywi±cie musimy wliczy¢ w to jakie± zakªócenia, wi¦c powiemy, »e nasz model opiszemy równaniem:

wzrost = b0+ b1∗ wiek + e gdzie e b¦dzie bª¦dem resztowym.

3. Je±li u»yjemy polecenia: >lm.r = lm(formula = wzrost ~ wiek) to ut- worzony w ten sposób obiekt mo»e by¢ u»yty do wywoªania innych polece«

R. Je±li np u»yjemy teraz polecenia summary(lm.r) otrzymamy zestaw- ienie:

> summary(lm.r) Call:

lm(formula = wzrost ~ wiek) Residuals:

Min 1Q Median 3Q Max

-0.27238 -0.24248 -0.02762 0.16014 0.47238 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 64.9283 0.5084 127.71 < 2e-16 ***

wiek 0.6350 0.0214 29.66 4.43e-11 ***

---Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 0.256 on 10 degrees of freedom Multiple R-squared: 0.9888, Adjusted R-squared: 0.9876 F-statistic: 880 on 1 and 10 DF, p-value: 4.428e-11 4. Gdy u»yjemy polecenia coef(lm.r) otrzymamy zestawienie:

(Intercept) 64.9283

wiek 0.6350

5. Powiemy wtedy, »e nasz model regresji przedstawia si¦ nast¦puj¡co:

wzrost = 64.9283 + 0.6350 ∗ wiek

6. generujemy wykres

> plot(lm.r)

Waiting to confirm page change...

Klikaj¡c enter dostajemy kolejne wykresy...

7. obliczenie korelacji

> cor(wzrost,wiek) [1] 0.994366 8. obliczenie korelacji R²

> cor(wzrost,wiek)^2 [1] 0.988764

9. wyrysowanie wykresów rezyduów dla konkretnych zmiennych

> plot(wzrost,lm.r$residuals)

76 78 80 82

−0.20.00.20.4

wzrost

lm.r$residuals

Rysunek 1: wykres reszt dla zmiennej wzrost

18 20 22 24 26 28

−0.20.00.20.4

wiek

lm.r$residuals

Rysunek 2: wykres reszt dla zmiennej wiek

18 20 22 24 26 28

−0.20.00.20.4

wiek

lm.r$residuals

Rysunek 3: dodanie linii poziomej

Wówczas wynikiem b¦d¡ odpowiednio rysunki: 1 i 2.

10. dodaje poziom¡ linie

> abline(h = 0)

Wówczas wynikiem b¦dzie rysunek 3.

11. Zastosujemy zwykª¡ funkcje plot z argumentem b¦d¡cym rezultatem funkcji lm. Generuje ona 4 wykresy dla regresji. Zmieniaj¡c parametr mfrow uzyskujemy wszystkie wykresy w jednym oknie

> par(mfrow = c(2, 2));plot(lm.r);par(mfrow = c(1, 1)) Wówczas wynikiem b¦dzie rysunek 4.

12. Gdy chcemy u»y¢ R i regresji do predykcji:

U»ywamy do tego funkcji predict, której ogólna formuªa jest nast¦puj¡ca:

predict(model, data.frame(pred = new pred), level = 0.95, interval = ''confidence'') gdzie pred to zmienna dla której chcemy sprawdzi¢ now¡ warto±¢ by poda¢ jej

warto±¢ dla zmiennej obja±nianej. Podajemy tu tak»e odpowiedni przedziaª za- ufania.

Np chcemy przewidzie¢ jaka b¦dzie warto±¢ wzrostu dla wieku 28.5, to u»ywamy polecenia:

predict(lm.r,data.frame(wiek = 28.5), level = 0.9, interval = "confidence")

77 79 81 83

−0.20.2

Fitted values

Residuals

Residuals vs Fitted

3 8

10

−1.5 −0.5 0.5 1.5

−1.00.01.02.0

Theoretical Quantiles

Standardized residuals

Normal Q−Q

3 8

1

77 79 81 83

0.00.40.81.2

Fitted values

Standardized residuals

Scale−Location

3 8 1

0.00 0.10 0.20 0.30

−1012

Leverage

Standardized residuals

Cook’s distance

0.5 0.5 1

Residuals vs Leverage

3

10 1

Rysunek 4: rezultat funkcji lm - 4 wykresy na jednym

lm.r lwr upr

1 83.02483 82.78911 83.26054

gdzie jak widzimy mamy podane dolne i górne zakresy przedziaªów ufno±ci.

4 Zadania do wykonania

‚wiczenie 0: Dla zbioru:

> liczbaGodzin <-c(18,19,20,21,22,23,24,25,26,27,28,29,45,50,61,77)

> wynagrodzenie <-c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8, 82.8,83.5, 90,94,112,121)

wykonaj polecenia:

• podaj równanie regresji dla zmiennych liczbaGodzin oraz wynagrodzenie

?

• predykcja: jaka b¦dzie warto±¢ wynagrodzenia dla liczbyGodzin równej 39?

‚wiczenie 1. Wykonaj analiz¦ regresji liniowej dla par zmiennych x i y z pliku anscombedost¦pnego w pakiecie R. Porównaj wyniki otrzymane w tabelach oraz wykresy rozrzutu z zaznaczonymi prostymi regresji. Czy we wszystkich przy- padkach prosta regresji dobrze oddaje zale»no±¢ mi¦dzy zminnymi? ‚wicze- nie 2. Zaªó»my, ze sporzadzili±my krzyw¡ kalibracyjn¡ oznaczania pewnego zwi¡zku dla st¦»e« 1, 2, ..., 9, 10 mg/ml. Ka»de oznaczenie powtarzano trzykrot- nie. Wyniki zapisa¢ nale»y w dwóch wektorach x i y:

> x <-c(1, 1, 1, 2, 2, 2, 3 ,3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10)

> y <- c(2.00, 2.01, 2.00, 2.42, 2.42, 2.43, 2.72, 2.74, 2.74, 3.00, 3.00, 2.98, 3.23, 3.24, 3.23, 3.44,

3.44, 3.45, 3.66, 3.65, 3.62, 3.82, 3.85, 3.83, 4.00, 4.00, 4.00, 4.18, 4.16, 4.18)

. Wykonaj analiz¦ regresji dla danego zbioru.