TEMAT 2
Wyrównanie odporne na błędy grube
Wpływ danych na modele prostej regresji

Król Mariusz
GiK, GIP 1
Rok akad. 2015/16

SPRAWOZDANIE TECHNICZNE

Dane formalno - prawne:

 Zleceniodawca: Akademia Górniczo-Hutnicza im. S. Staszica w Krakowie,
Wydział Geodezji Górniczej i Inżynierii Środowiska;
 Wykonawca: Mariusz Król, GiK - GIP,
 Okres wykonywania zlecenia:
- termin rozpoczęcia prac: 14.10.2015r.
- termin zakończenia prac: 21.10.2015r.
 Przedmiot zlecenia: Wyrównanie odporne na błędy grube. Wpływ danych na modele prostej regresji.
 Wykorzystane oprogramowanie: The R Project for Statistical Computing.

Opracowanie tematu:

1. Przygotowanie danych: W programie „R” przygotowano dane dla 30 punktów w postaci równo rozłożonych współrzędnych x,y będących obrazem funkcji y=x. Do wyznaczonych współrzędnych wprowadzono losowy szum o wartości 0.5.
 Dodatkowo wprowadzono ręcznie 6 punktów, których współrzędne będą obrazować błędy grube niezbędne do wykonania tego doświadczenia.

len=(30)
x=30*runif(len)
y=x+rnorm(len,0,0.5)

OTRZYMANE DANE (wraz z ręcznie wpisanymi współrzędnymi-kolor pomarańczowy):

nr	x	y
1	3.923731	3.196332
2	20	0
3	13.01354	13.70935
4	12.72527	12.54677
5	21.77191	22.22374
6	7.237774	7.547028
7	25	7
8	21.04206	20.58468
9	12.8966	12.87582
10	20.4024	19.80684
11	19.3239	18.23782
12	24.3895	24.2809
13	9.131399	9.235892
14	18.81731	19.56509
15	20.42104	20.80813
16	7	21
17	25.80125	25.88693
18	12.56871	13.09454
19	1.500314	2.003454
20	29.40618	29.0336
21	23.73113	23.31034
22	19.84401	20.2216
23	30	5
24	22.10374	20.88234
25	25.03414	25.74214
26	14.77456	14.35524
27	23.17626	22.87097
28	10.10339	10.89825
29	29.0214	29.2469
30	11.0347	10.62019
31	22.13243	22.64489
32	5	31
33	28.19575	27.8548
34	7.488908	7.246011
35	0	20
36	10.33543	9.984817

2. Wczytanie danych: Po przygotowaniu danych w postaci plików tekstowych ze współrzędnymi x,y i zmianie domyślnego katalogu w programie „R” na ten, w którym znajdują się wymienione wyżej pliki tekstowe dokonano wczytania danych.

x =scan(file="X.txt")
y=scan(file="Y.txt")

 Zobrazowanie danych wejściowych na wykresie

plot(y~x)

3. Model regresji liniowej: stworzono przy tym ramkę danych „ds3”

ds3=data.frame(x=x,y=y)
m.lm=lm(y~x,data=ds3)
summary(m.lm)

 Otrzymano następujące wyniki regresji liniowej:

Call:
lm(formula = y ~ x, data = ds3)

Residuals:
Min 1Q Median 3Q Max
-18.220 -3.551 1.548 3.593 19.685

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.0130 2.8365 3.177 0.00316 **
x 0.4603 0.1507 3.054 0.00436 **
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.493 on 34 degrees of freedom
Multiple R-squared: 0.2153, Adjusted R-squared: 0.1922
F-statistic: 9.33 on 1 and 34 DF, p-value: 0.004363

Z powyższych zestawień wyników można odczytać wartości parametrów:
a = 0,4603, b = 9,0130 będących rozwiązaniem funkcji regresji liniowej y=ax+b.

 Wydruk obliczonej funkcji regresji:

l ayout(1)
plot(x~y)
lines(lowess(y,x))

 Badanie reszt:

p loop=par(mfrow=c(2,2),pty="s")
plot(m.lm)

 Komentarz:
Wykres „Residuals vs Fitted” obrazuje jak ma się wpasowanie do reszt. Widać dobrze odstające punkty.
Wykres „Scale-Location” to wyciągnięty pierwiastek z wartości pierwszego wykresu.
Wykres „Normal Q-Q” to reszty zestandaryzowane z kwantylami. Gdy model jest bezbłędny punkty są wzdłuż jednej prostej. W tym przypadku widać błędy w danych, są to odstające punkty.
Wykres „Residuals vs Leverage” porównuje zestandaryzowane reszty z obliczonymi resztami. Linia przerywana to odległość Cook’a. Jest ona łączną miarą wpływu poszczególnych obserwacji na linię regresji.

 Obliczenie odległości Cook’a i wartości zestandaryzowanych reszt:

library(MASS)
dl=cooks.distance(m.lm)
r=stdres(m.lm)
a=cbind(ds3,dl,r)
a[dl>4/36,]

x y dl r
23 30 5 0.3373892 -2.503311
32 5 31 0.3505590 2.746406
35 0 20 0.2098580 1.584105

Komentarz:
Punkty 23,32,35 to obserwacje, które posiadają błędy grube..Są to jedne z tych punktów które zostały ręcznie dopisane. Wyświetlone zostały tylko trzy z sześciu dlatego, że prawdopodobnie ich wpływ na regresję jest największy.

 Wyznaczenie wartości bezwzględnych reszt zestandaryzowanych:

rabs=abs(r)
a=cbind(ds3,dl,r,rabs)
asorted=a[order(-rabs),]
asorted[1:10,]

x y dl r rabs
32 5.000000 31.000000 0.35055898 2.7464056 2.7464056
23 30.000000 5.000000 0.33738917 -2.5033107 2.5033107
2 20.000000 0.000000 0.09983541 -2.4708931 2.4708931
7 25.000000 7.000000 0.09890090 -1.8555669 1.8555669
35 0.000000 20.000000 0.20985799 1.5841054 1.5841054
16 7.000000 21.000000 0.05302037 1.2111836 1.2111836
19 1.500314 2.003454 0.08504384 -1.0977331 1.0977331
1 3.923731 3.196332 0.06068967 -1.0698691 1.0698691
29 29.021398 29.246899 0.04405344 0.9602024 0.9602024
20 29.406176 29.033595 0.04125971 0.9075739 0.9075739

Zestawiono 10 punktów, których wartości bezwzględne reszt zestandaryzowanych są największe. Pierwsze 6 z nich to wpisane ręcznie wartości z błędami grubymi.

 Wykres odległości Cook’a

library(car)
influencePlot(lm(y~x), main="InfluencePlot",sub="Promień kółka jest proporcjonalny do odległości Cooka")

StudRes Hat CookD
32 3.067251 0.08504727 0.5920802
35 1.621625 0.14329151 0.4581026

influence.measures(lm(y~x))

Komentarz:
Gwiazdki na końcu wiersza wskazują, że dany punkt zawiera błąd gruby. Świadczą o tym największe odległości Cook’a.

lm(formula=y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x
9.0130 0.4603

4. Regresja liniowa z „ręcznym” usunięciem danych: zdecydowano się usunąć następujące numery punktów: 2,7,16,23,32,35 wobec których zachodzi podejrzenie występowania błędów grubych w obserwacjach.
Po usunięciu tych punktów ponownie policzono odległości Cook’a, wartości zlinearyzowanych reszt oraz wartości bezwzględne reszt zestandaryzowanych.

m.lm2=lm(y~x,subset=-c(2,7,16,23,32,35))
summary(m.lm2)

Call:
lm(formula = y ~ x, subset = -c(2, 7, 16, 23, 32, 35))

Residuals:
Min 1Q Median 3Q Max
-1.15903 -0.36742 -0.02765 0.43371 0.79182

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09943 0.25170 0.395 0.696
x 0.99268 0.01329 74.667 <2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5471 on 28 degrees of freedom
Multiple R-squared: 0.995, Adjusted R-squared: 0.9948
F-statistic: 5575 on 1 and 28 DF, p-value: < 2.2e-16

Z powyższych zestawień wyników można odczytać wartości parametrów:
a = 0,99268, b = 0,09943 będących rozwiązaniem funkcji regresji liniowej y=ax+b.

ploop=par(mfrow=c(2,2),pty="s")
plot(m.lm2)

dl.2=cooks.distance(m.lm2)
r2=stdres(m.lm2)
s=seq(1,36,b=1)
s=s[-c(2,7,16,23,32,35)]
a2=cbind(s,dl.2,r2)
a2[dl.2>4/33,]

s dl.2 r2
1.0000000 0.2019205 -1.5734183

rabs2=abs(r2)
a2=cbind(s,dl.2,r2,rabs2)
asorted2=a2[order(-rabs2),]
asorted2[1:10,]

s dl.2 r2 rabs2
24 24 0.11484968 -2.1697452 2.1697452
11 11 0.06964579 -1.9433998 1.9433998
1 1 0.20192054 -1.5734183 1.5734183
25 25 0.08193567 1.4992595 1.4992595
14 14 0.03827684 1.4624554 1.4624554
28 28 0.07300768 1.4541711 1.4541711
3 3 0.03903677 1.2934620 1.2934620
31 31 0.02837793 1.0765868 1.0765868
10 10 0.02085209 -1.0173308 1.0173308
18 18 0.02323247 0.9706964 0.9706964

Komentarz:
W tym przypadku brak jest gwiazdek, tak więc odległości Cook’a nie są duże. Odrzucenie obserwacji z błędami grubymi sprawiło, że na wykresach próbki rozłożone są równomiernie. W odróżnieniu od poprzedniego przypadku teraz brak jest większych grup punktów i oddalonych od nich pojedynczych punktów, które odstają od reszty danych.

5. Regresja liniowa z modelem Hubera

library(MASS)
m.hub=rlm(y~x)
summary(m.hub)

Call: rlm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-24.34684 -0.57019 -0.03898 0.53549 25.39636

Coefficients:
Value Std. Error t value
(Intercept) 0.8550 0.3580 2.3880
x 0.9497 0.0190 49.9247

Residual standard error: 0.8086 on 34 degrees of freedom

Z powyższych zestawień wyników można odczytać wartości parametrów:
a = 0,9497, b = 0,8550 będących rozwiązaniem funkcji regresji liniowej y=ax+b.

ploop=par(mfrow=c(2,2),pty="s")
plot(m.hub)

plot(m.hub$w, ylab="Wagi Hubera")

a=cbind(ds3,dl,r,m.hub$w)
asorted3=a[order(m.hub$w),]
asorted3[1:10,]

x y dl r m.hub$w
32 5.000000 31.000000 0.3505589836 2.7464056 0.04282529
23 30.000000 5.000000 0.3373891661 -2.5033107 0.04467126
2 20.000000 0.000000 0.0998354119 -2.4708931 0.05479039
35 0.000000 20.000000 0.2098579937 1.5841054 0.05681113
7 25.000000 7.000000 0.0989008954 -1.8555669 0.06180146
16 7.000000 21.000000 0.0530203650 1.2111836 0.08058417
1 3.923731 3.196332 0.0606896698 -1.0698691 0.78442610
25 25.034141 25.742143 0.0147223018 0.7143621 0.97788374
3 13.013539 13.709349 0.0005415254 -0.1757302 1.00000000
4 12.725274 12.546774 0.0017980781 -0.3157111 1.00000000

Komentarz:
Ten model regresji, w odróżnieniu od metody Cook’a, bardzo wyraźnie pokazuje błędy grube w obserwacjach. Wykresy bardzo obrazowo ukazują te punkty, ale co jest najistotniejsze, widać jak małe wartości wag przyjmowane są w tym modelu dla wartości odstających (wartości m.hub$w)

6. Regresja z modelem Bisquare

m.bi=rlm(y~x,method='MM')
m.bi

Call:
rlm(formula = y ~ x, method = "MM")
Converged in 3 iterations

Coefficients:
(Intercept) x
0.1018005 0.9929018

Degrees of freedom: 36 total; 34 residual
Scale estimate: 0.792

Z powyższych zestawień wyników można odczytać wartości parametrów:
a = 0,9929018, b = 0,1018005 będących rozwiązaniem funkcji regresji liniowej y=ax+b.

plot(m.bi$w, ylab="Wagi Bis",type="p",lwd=10)

 Podsumowanie
Badane modele regresji spełniają swoją funkcję właściwie wykrywają obserwacje obarczone błędami grubymi.
Najlepszą metodą (z punktu widzenia czytelności i jasności otrzymanych wyników) jest regresja liniowa z modelem Hubera. W tym modelu, na podstawie wag Hubera, doskonale można odczytać, które obserwacje są znacząco odstające w zbiorze danych a które nie. Punkty z błędami grubymi otrzymały bardzo małe wagi w przedziale 0,4 - 0,8, natomiast obserwacje poprawne otrzymały wagi równe 1 lub zbliżone do jedności. Dzięki tym wartością szybko można wychwycić nieprawidłowe obserwacje w analizowanym zbiorze danych.

Prosta regresji w opracowanym temacie zbliżona była swym kształtem do linii prostej. Wynika to z faktu iż wprowadzone błędy grube do obserwacji, ich wpływ na cały zbiór miar charakter symetryczny, był to szczególny przypadek, kiedy to błędy grube, mimo że znajdowały się w zbiorze nie wpływały znacząco na proces regresji liniowej.

Kraków 20.10.2015 Król Mariusz

13