1.1 Obserwacje odstaj¡ce

Zagadnienia regresji. Cz¦±¢ II

Konspekt do zaj¦¢: Statystyczne metody analizy danych

Agnieszka Nowak-Brzezi«ska 16 listopada 2009

Konsultacje z Panem Profesorem Jackiem Koronackim przy- czyniªy si¦ do kilku zmian w sporz¡dzonym konspekcie (z dnia 4 listopada). Niniejsza wersja niech b¦dzie zatem obo- wi¡zuj¡ca.

Z analiz¡ regresji wi¡»¡ si¦ nast¦puj¡ce zagadnienia:

1. Maj¡c równanie prostej y = β1x + β0 powiemy, »e β1 jest nachyleniem linii regresji,za± β0 punktem przeci¦cia linii regresji z osi¡ x (wyrazem wolnym). Wyraz wolny β0 jest miejscem na osi y, gdzie linia regresji przecina t¦ o±, czyli jest to przewidywana warto±¢ zmiennej odpowiedzi, gdy zmienna obja±niaj¡ca jest równa 0.

2. wspóªczynnik determinacji - okre±la stopie«, w jakim linia regresji naj- mniejszych kwadratów wyja±nia zmienno±¢ obserwowanych danych. Do- kªadniej mówi¡c mierzy on stopie« dopasowania regresji jako przybli»enia liniowej zale»no±ci pomi¦dzy zmienn¡ celu a zmienn¡ obja±niaj¡c¡:

R²= SSR SST

, gdzie SSR to regresyjna suma kwadratów (SSR =P_n

i=1(ˆy − ¯y)²) za±

SST to caªkowita suma kwadratów (SST =P_n

i=1(y − ¯y)²). Maksymalna warto±¢ R²(R²= 1) jest osi¡gana, gdy regresja idealnie pasuje do danych, co ma miejsce wówczas, gdy ka»dy z punktów danych le»y dokªadnie na oszacowanej linii regresji. Za±, warto±¢ minimalna (0) b¦dzie oznacza¢

sªabe dopasowanie regresji do zbioru danych.

3. obsewacje odstaj¡ce oraz obserwacje wpªywowe - czyli takie obserwacje, które mog¡, ale nie musz¡, wywiera¢ nadmierny nacisk na wyniki regresji.

Zostan¡ przedmiotem niniejszych zaj¦¢.

1 Obserwacje odstaj¡ce a obserwacje wpªywowe

1.1 Obserwacje odstaj¡ce

Obserwacja odstaj¡ca czyli nietypowa (ang. outlier) jest obserwacj¡, która nie speªnia równo±ci

Yi= β0+ β1xi+ ²i

dla i = 1, 2, . . . , n, gdzie ²i s¡ niezale»nymi zmiennymi losowymi o tym samym rozkªadzie ze ±redni¡ 0 i wariancj¡ σ². Mówimy cz¦sto, »e obserwacja taka nie nale»y do modelu prostej regresji.

Chcemy wykrywa¢ zmienne odstaj¡ce gdy» mog¡ one znacz¡co wpªwa¢ na posta¢ prostej regresji (prostej MNK):

b0+ b1x dla której warto±¢ sumy: P_n

i=1(yi− ˆyi)² a wi¦c i sumyP_n

i=1(yi− (b0+ b1xi))² jest (ma by¢) najmniejsza.

Je±li obserwacja wpªywa na zmian¦ wspóªczynnika kierunkowego prostej, na- zwiemy j¡ obserwacj¡ wpªywow¡. Przy tym powiemy, »e je±li warto±¢ zmiennej obja±niaj¡cej dla danej obserwacji znacz¡co odbiega od typowych warto±ci tej zmiennej to uznamy j¡ za potencjalnie wpªywow¡.

1.1.1 Obserwacje odstaj¡ce dla zbioru wielu zmiennych obja±niaj¡- cych

Je±li analizujemy tylko pojedyncze zmienne obja±niaj¡ce, to identykacja ob- serwacji odstaj¡cych jest do±¢ prosta. Wystarczy generowa¢ wykresy rozrzutu b¡d¹ histogramy.

Je±li chcemy szuka¢ obserwacji odstaj¡cych globalnie (nie dla pojedynczej zmiennej obja±niaj¡cej ale dla wielu) wówczas mo»emy analizowa¢ rezydua lub rezydua studentyzowane, i w±ród nich szuka¢ warto±ci odstaj¡cych.

1. Maj¡c wektor warto±ci resztowych(rezyduów) e = (e1, e2, . . . , en), gdzie warto±¢ resztowa ei = Yi− ˆYi (ei = Yi− (b0xi+ b1)) powiemy, »e bª¡d standardowy rezyduum ei jest równy:

SEei= S ∗ s

1 − (1

n+ (xi− ¯x)² P_n

i=1(xi− ¯x)²).

Wtedy studentyzowana warto±¢ resztowa b¦dzie odpowiada¢ warto±ci ri= e_i

SEei

2. Sporz¡dzaj¡c wykres warto±ci studentyzowanych rezyduów ri wzgl¦dem ich indeksu b¦dziemy potrali rozpoznawa¢ te du»e warto±ci, które przy- puszczalnie b¦d¡ odstaj¡cymi.

Podsumowuj¡c powiemy, »e nowa obserwacja b¦dzie punktem odstaj¡cym je-

±li b¦dzie si¦ cechowa¢ du»¡ warto±ci¡ studentyzowanej (standaryzowanej) reszty.

W praktyce, obserwacje odstaj¡ce to takie, których warto±¢ bezwzgl¦dnych stu- dentyzowanych reszt przekracza 2.

1.2 Obserwacje wpªywowe

Obserwacja jest wpªywowa (ang. inuential), je±li jej obecno±¢ wpªywa na pro- st¡ regresji, w taki sposób, »e zmienia si¦ wspóªczynnik kierunkowy tej prostej.

Inaczej powiemy, »e je±li obserwacja jest wpªywowa to inaczej wygl¡da prosta regresji w zale»no±ci od tego czy ta obserwacja zostaªa uj¦ta w zbiorze, czy te»

nie (zostaªa usuni¦ta).

1.2.1 Identykacja obserwacji wpªywowych

W praktyce, je±li obserwowana warto±¢ le»y w I-ym kwartylu rozkªadu (czyli ma warto±¢ mniejsz¡ ni» 25 centyl), to mówimy, »e ma ona maªy wpªyw na regresj¦.

Obserwacje le»¡ce mi¦dzy I a III kwartylem nazywamy wpªywowymi.

Mówimy tak»e, »e czynnik¹_n+P_n^(xi−¯x)2

i=1(xi−¯x)² we wzorze na SEeito tzw. wpªyw tej obserwacji (czasami nazywany w literaturze d¹wigni¡). Zwykle obserwacje cechuj¡ce si¦ wysok¡ warto±ci¡ d¹wigni b¦d¡ uznawane za wpªywowe. Dodat- kowo powiemy, »e nawet je±li obserwacja jest odstaj¡ca, ale ma maª¡ warto±¢

wpªywu to uznamy, »e nie jest ona wpªywowa.

Wpªyw i-tej obserwacji hi = _n¹ + P_n^(xi−¯x)2

i=1(xi−¯x)² b¦dziemy okre±la¢ jako od- st¦pstwo obserwacji xi od ¯x. Wiadomo, »e dla modelu, który ma p parametrów powiemy, »eP_n

i=1hi= poraz dla ka»dego i 1 ≥ hi≥ ¹_n. To oznacza, »e typowa warto±¢ wpªywu hi nie powinna przekracza¢ warto±ci ^p_n. Je±li za± warto±¢ ta dla analizowanej i-tej obsewacji przekracza warto±¢ ^2p_n (a dla maªych prób ^3p_n) wówczas tak¡ zmienn¡ uznamy za potencjalnie wpªywow¡.

Tak naprawd¦ warto±¢ wpªywu dla i-tej obserwacji zale»y jedynie od (xi−

¯

x)². Czyli, im wi¦ksza ta ró»nica (podnoszona do kwadratu) tym wi¦ksza war- to±¢ wpªywu.

Inny sposób na wykrycie obserwacji wpªywowych to pomiar odlegªo±ci Co- oka, w której wykorzystujemy tzw. modykowane rezydua. Usuwaj¡c obserwa- cj¦, któr¡ chcemy uzna¢ za wpªywow¡ ze zbioru obserwacji, i obliczaj¡c ró»nic¦

Yˆj i ˆYj(i), je±li b¦dzie ona wysoka to wtedy powiemy, »e obserwacja (usuni¦ta) jest wpªywowa.

Odlegªo±¢ Cooka mierzy poziom wpªywu obserwacji, uwzwgl¦dniaj¡c za- równo wielko±¢ reszty, jak i wysoko±¢ wpªywu dla tej obserwacji. Dla i-tej obserwacji odlegªo±¢ Cooka jest obliczana jako:

Di= P_n

j=1( ˆYj− ˆY_j(i))² pS² = e²_i

pS² hi

(1 − hi)²

, gdzie ˆY_j(i) jest obserwacj¡ przewidywan¡ dla j-tej obserwacji obliczon¡ na podstawie danych z usuni¦t¡ obserwacj¡ i-t¡, za± ˆYj b¦dzie warto±ci¡ przewi- dywan¡ dla j-tej obserwacji gdy i-tej obserwacji nie usuni¦to. Du»a warto±¢ Di

mówi o du»ym wpªywie usuni¦cia i-tej obserwacji, a tym samym obserwacj¦ i-t¡

uznajemy za wpªywow¡.

2 Przykªad analizy krok po kroku

Zaªó»my,»e mamy do analizy dane zbiorów a i b, gdzie a oznacza liczb¦ godzin pracy, a b - wynagrodzenie za t¦ prac¦. Komendy do wywoªania modelu regresji dla takich danych wygl¡daj¡ nast¦puj¡co:

> a<-c(2,2,3,4,4,5,6,7,8,9)

> b<-c(10,11,12,13,14,15,20,18,22,25)

> model = lm(b ~ a)

> summary(model) Efektem b¦dzie:

Call:

lm(formula = b ~ a) Residuals:

Min 1Q Median 3Q Max

-2.000e+00 -7.500e-01 4.876e-16 7.500e-01 2.000e+00 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 6.0000 0.9189 6.529 0.000182 ***

a 2.0000 0.1667 12.000 2.14e-06 ***

---Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 1.225 on 8 degrees of freedom

Multiple R-squared: 0.9474, Adjusted R-squared: 0.9408 F-statistic: 144 on 1 and 8 DF, p-value: 2.144e-06

Z tego wynika jednoznacznie, »e równanie regresji wygl¡da¢ b¦dzie nast¦puj¡co:

ˆ

y = 2x + 6

Zreszt¡ wspóªczynniki β0i β1mo»na wyznaczy¢ wywoluj¡c komend¦ coef():czego efekt b¦dzie nast¦puj¡cy:

> coef(model) (Intercept) a

6 2

Teraz chc¡c pozna¢ warto±ci rezyduów dla ka»dej obserwacji wystarczy wy- woªa¢: resid(model) czego efekt b¦dzie nast¦puj¡cy:

1 2 3 4 5

.123139e-15 1.000000e+00 5.447839e-16 -1.000000e+00 4.753475e-16 6 7 8 9 10

-1.000000e+00 2.000000e+00 -2.000000e+00 4.998166e-16 1.000000e+00 Widzimu wi¦c warto±ci poszczególnych reszt dla kolejnych obserwacji w zbiorze a. Aby wyznaczy¢ warto±ci oczekiwane dla ka»dej obserwacji u»ywamy ko- mendy:

fitted(model), czego efektem b¦d¡ warto±ci:

> fitted(model)

1 2 3 4 5 6 7 8 9 10

10 10 12 14 14 16 18 20 22 24

Jak ªatwo zauwa»y¢ b¦d¡ to po prostu warto±ci ˆyi, a wi¦c dla ka»dej z 10- ciu obserwacji w zbiorze a (dla zmiennej obja±nianej a) za pomoc¡ rowniania regresji: ˆy = 2x + 6 jeste±my w stanie wskaza¢ warto±ci zmiennej obja±niaj¡cej - nie te rzeczywiste (podane w zbiorze b) lecz te oczekiwane obliczone z wzoru.

We¹my np piersz¡ obserwacj¦, która odpowiednio dla zbiorów a i b przyjmuje warto±ci 2 i 10. A wi¦c tutaj x = 2 a y = 10. Je±li teraz sprawdzimy jaka bylaby oczekiwana warto±ci y a wi¦c ˆy to podstawimy odpowiednio te warto±ci do równania: ˆy = 2x + 6 = 2 ∗ 2 + 6 = 10 i otrzymamy warto±¢ 10. W tym konkretnym przypadku warto±ci rzeczywiste i te oczekiwane s¡ identyczne. Nie zawsze tak jest, i je±li nie jest to mówimy o tzw. bª¦dzie. Czasami okre±lamy go poj¦ciem reszty (ang. residua). Tak b¦dzie w przypadku obserwacji 2, gdzie odpowednio a2 = 2 i b2 = 11. Dla a = 2 warto±ci¡ oczekiwan¡ jest przecie»

(zgodnie ze wzorem) warto±¢ 10 a nie 11, st¡d dla tej warto±ci wyst¡pi ró»nica mi¦dzy warto±cia rzeczywist¡ a t¡ oczekiwan¡ o warto±ci 1. Pami¦tajmy, ze potem ujmujemy wszystkie wyst¦puj¡ce ró»nice, i sumujemy ich kwadraty (tutaj nawi¡zanie do metody MNK z poprzednich zaj¦¢ a przede wszystkim wykªadu).

2.1 Identykacja obserwacji odstaj¡cych

Jedn¡ z fukcji pozwalaj¡cych odnajdywa¢ obserwacje odstaj¡ce jest rstudent.

Wywoªuj¡c komendy:

> jack<-rstudent(model)

> jack[which.max(abs(jack))]

uzyskamy informacje o tym, które obserwacje byªy odstaj¡ce.

-2.1764298

Mo»na to zobrazowa¢ wykresem. Nast¦puj¡cy kod R-a:

> d<-rstudent(model)

> plot(d,ylab="Jacknife Residuals",main="Jacknife Residuals") da w wyniku wykres 1 na którym widoczne s¡ punkty odstaj¡ce.

2 4 6 8 10

−2−1012

Jacknife Residuals

Index

Jacknife Residuals

Rysunek 1: Wykres obserwacji odstaj¡cych

Doskonale do wykrywania obserwacji odstaj¡cych nadaje si¦ biblioteka car, w ktorej jest m.in. funkcja outlier.test.

Jej u»ycie dla naszego modelu o nazwie model wygl¡da nast¦puj¡co:

> library(car)

> outlier.test(model)

a w efekcie spowoduje, »e otrzymamy nast¦puj¡ce warto±ci:

max|rstudent| = 2.176429, degrees of freedom = 7, unadjusted p = 0.06598772, Bonferroni p = 0.6598772 Observation: 8

Otrzymujemy zatem informacj¦, »e obserwacj¡ nietypow¡ jest tutaj obser- wacja numer 8 przy 7 stopniach swobody (bo T − m − 2).

Bardzo wa»ne dla wykrycia obserwacji odstaj¡cych s¡ tzw. studentyzowane reszty. Wykres dla nich mo»emy wykona¢ wywoªuj¡c komend¦:

> qq.plot(model, main="QQ Plot")

Efektem b¦dzie wykres, na którym zobaczymy rozkªad obserwacji mi¦dzy I i III kwartylem, st¡d nazwa wykresy "mi¦dzykwartylny"(rysunek 2).

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

−2−1012

QQ Plot

t Quantiles

Studentized Residuals(model)

6

7

Rysunek 2: Wykres mi¦dzykwartylny dla obserwacji nietypowych

Efektem b¦dzie wykres z tytuªem QQP lot, na którym klikaj¡c myszk¡ w dowolny punkt, znaleziony i zidentykowany zostanie najbli»szy mu punkt w zbiorze danych (po wyborze opcji ST OP zostanie zwrócony indeks tej obserwa- cji).

2.2 Identykacja obserwacji wpªywowych

Mówi si¦, obserwacje wpªywowe mo»na wyznaczy¢ u»ywaj¡c tzw. odleglo±ci Co- oka. Je±liby±my chcieli aby punktem granicznym byª poziom, gdy odlegªo±¢ jest wi¦ksza ni» warto±¢ 4/(n − k − 1). Najpierw zatem przypisujemy do zmiennej cutoff poziom

cutoff <- 4/((nrow(model)-length(model\$coefficients)-2)) Wówczas wywoªuj¡c komend¦

> plot(model, which=4, cook.levels=cutoff) otrzymamy wykres 3

2 4 6 8 10

0.00.10.20.30.4

Obs. number

Cook’s distance

lm(b ~ a) Cook’s distance

10 8

7

Rysunek 3: Wykres obserwacji wpªywowych

Teraz je±li chcemy pozna¢ wykres dla samych danych wpªywowych mozemy u»y¢ komendy:

> influencePlot(model, main="Influence Plot",sub="Circle size is proportial to Cook's Distance")

którego efektem b¦dzie wykres 4:

0.10 0.15 0.20 0.25 0.30 0.35 0.40

−2−1012

Influence Plot

Circle size is proportial to Cook’s Distance Hat−Values

Studentized Residuals

2

4 6

7

8

10

Rysunek 4: Wykres obserwacji wpªywowych z zaznaczeniem odlegªo±ci Cooka

Do wykrycia obserwacji wpªywowych mo»emy tak»e u»y¢ funkcji

> influence.measures(model) której efekty b¦dzie nast¦puj¡cy:

Influence measures of lm(formula = b ~ a) :

dfb.1_ dfb.a dffit cov.r cook.d hat inf

1 0.0000 0.00e+00 -1.14e-15 1.781 7.45e-31 0.267 * 2 0.5570 -4.52e-01 5.71e-01 1.399 1.65e-01 0.267 3 0.0000 0.00e+00 2.10e-16 1.581 2.52e-32 0.174 4 -0.2337 1.24e-01 -3.13e-01 1.215 5.08e-02 0.119 5 0.0000 0.00e+00 1.42e-16 1.482 1.15e-32 0.119 6 -0.1187 2.05e-18 -2.82e-01 1.195 4.12e-02 0.100 7 0.0217 2.99e-01 7.57e-01 0.573 2.03e-01 0.119 8 0.2719 -6.52e-01 -9.99e-01 0.563 3.40e-01 0.174 9 0.0000 0.00e+00 2.69e-16 1.781 4.13e-32 0.267 * 10 -0.4910 7.42e-01 8.58e-01 1.607 3.62e-01 0.396

>

Jak wida¢, ostatnia kolumna wskazuje na obserwacje wpªywowe zaznaczaj¡c przy nich symbol ∗. Z naszych danych wynika, »e w zbiorze 10 obserwacji mamy 2 wpªywowe. S¡ to obserwacie 1 i 9. Funkcja lm.influence dostarcza informacji o 4 parametrach: hat, coefficients, sigma i wt.res (weighted residuals).

> lm.influence(lm(model))

$hat 1 2 3 4 5 6 7

8 9 10

0.2666667 0.2666667 0.1740741 0.1185185 0.1185185 0.1000000 0.1185185 0.1740741 0.2666667 0.3962963

$coefficients

(Intercept) a

1 0.00000000 0.000000e+00 2 0.51515152 -7.575758e-02 3 0.00000000 0.000000e+00 4 -0.21848739 2.100840e-02 5 0.00000000 0.000000e+00 6 -0.11111111 3.485662e-19 7 0.01680672 4.201681e-02 8 0.20627803 -8.968610e-02 9 0.00000000 0.000000e+00 10 -0.44785276 1.226994e-01

$sigma

1 2 3 4 5 6 7 8

1.309307 1.232672 1.309307 1.245882 1.309307 1.247219 1.032486 1.011149

$wt.res

1 2 3 4 5 6 7 8 9 10 0 1 0 -1 0 -1 2 -2 0 1

3 Zadanie do wykonania

Dla zbioru: http://lib.stat.cmu.edu/DASL/Datafiles/Cereals.html

• wyznacz równanie regresji sugeruj¡c si¦ tym, »e chcemy wyznaczy¢ war- to±ci od»ywcze (rating) platków maj¡c dane zawarto±ci cukrów (sugar).

• jaka b¦dzie przewidywana warto±¢ od»ywcza pªatków z zerow¡ zawarto±ci¡

cukrów?

• oblicz i zinterpretuj wspóªczynnik korelacji

• o ile wzro±nie lub zmaleje wartos¢ od»ywcza pªatków je±li zawarto±¢ cu- krów wzro±nie o 1 gram ?

• wyznacz obserwacje nietypowe, wpªywowe i o wysokiej d¹wigni - je±li ist- niej¡. Podaj ich nazwy (identydikatory). Ile jest takich warto±ci ?

4 Bibliograa

Opracowanie przygotowano w oparciu o prace:

1. J. Koronacki i J. ‚wik, Statystyczne systemy ucz¡ce si¦, wyd. II, Exit 2008

2. J. Koronacki i J. Mielniczuk, Statystyka dla studentów kierunków tech- nicznych i przyrodniczych, WNT 2006

3. Daniel T. Larose, Metody i modele eksploracji danych, Tytuª oryginalny:

Data Mining Methods and Models, Wydawnictwo Naukowe PWN 2008 4. Redakcja naukowa: Marek Walesiak, Eugeniusz Gatnar, Statystyczna ana-

liza danych z wykorzystaniem programu R, Wydawnictwo Naukowe PWN 2009