Elementy analizy
korelacji i regresji
ANALIZA KORELACJI
" Czy między cechami zachodzi współzale\
ność w tym
sensie, \
e zmianom jednej cechy towarzyszÄ… zmiany
drugiej cechy?
" Jaki jest kierunek tej współzale\
ności?
" Jaka jest jej siła?
Problem badawczy
Czy istnieje związek korelacyjny między temperaturą wody
w pewnym jeziorze i ilością tlenu rozpuszczonego w wodzie?
Tlen
Temperatura rozpuszczony
Lp. (Ä…C) (mg/dm3)
1. 11,1 9,3
2. 11,3 9,6
3. 11,6 9,2
4. 11,7 7,3
5. 12,7 6,7
6. 13,8 6,0
7. 14,8 6,2
8. 15,2 6,2
9. 15,7 5,9
10. 16,3 4,3
> temperatura<-c(11.1,11.3,11.6,11.7,12.7,13.8,14.8,15.2,15.7,16.3)
> temperatura
[1] 11.1 11.3 11.6 11.7 12.7 13.8 14.8 15.2 15.7 16.3
> tlen<-c(9.3,9.6,9.2,7.3,6.7,6.0,6.2,6.2,5.9,4.3)
> tlen
[1] 9.3 9.6 9.2 7.3 6.7 6.0 6.2 6.2 5.9 4.3
Wykres rozrzutu
> plot(x=temperatura,y=tlen,
xlab="Temperatura wody",
ylab="Tlen rozpuszczony",
main="Wykres rozrzutu",
pch=19)
> help(points)#opis dla pch
11 12 13 14 15 16
Temperatura wody
9
8
7
Tlen rozpuszczony
6
5
Wykres rozrzutu
Wy\
szej temperaturze odpowiada
na ogół ni\
sza zawartość tlenu
Zale\
ność między temperaturą
i zawartością tlenu ma charakter
ujemny.
Między temperaturą i zawartością
tlenu zachodzi zale\
ność
korelacyjna o charakterze
ujemnym.
11 12 13 14 15 16
Temperatura wody
Zale\
ność korelacyjna między cechami występuje wtedy, gdy ustalonym
wartościom jednej cechy towarzyszą niejednakowe średnie wartości drugiej cechy.
9
8
7
Tlen rozpuszczony
6
5
Zale\
ność korelacyjna o charakterze dodatnim:
 małym wartościom x-ów odpowiadają średnio  małe wartości y-ków
 du\
ym wartościom x-ów odpowiadają średnio  du\
e wartości y-ków
y
y
x x
y
x
Zale\
ność korelacyjna o charakterze ujemnym:
 małym wartościom x-ów odpowiadają średnio  du\
e wartości y-ków
 du\
ym wartościom x-ów odpowiadają średnio  małe wartości y-ków
y
y
x x
y
x
Brak zale\
ności korelacyjnej
y y
x x
y
x
Zale\
ność korelacyjna liniowa o charakterze dodatnim
r = 0.18 r = 0.47
y y
x x
r = 0.83 r = 1
y y
x x
Zale\
ność korelacyjna liniowa o charakterze ujemnym
r = -0.18 r = -0.47
y y
x x
r = -0.83 r = -1
y y
x x
Współczynnik korelacji liniowej Pearsona
Szereg dwuwymiarowy
n
Lp. xi yi
"(x - x)(yi - y)
i
i=1
r =
1. x1 y1
n n
2. x2 y2
"(x - x)2 Å""(y - y)2
i i
M
i=1 i=1
M
M
n. xn yn
> cor(x=temperatura, y =tlen, method="pearson")
[1] -0.9038088
Własności współczynnika r
0 d" | r | d" 0,2
Bardzo słaba zale\
ność liniowa
" Współczynnik r ocenia siłę
0,2 < | r | d" 0,4
liniowej zale\
ności między
cechami. SÅ‚aba zale\
ność liniowa
" Przyjmuje wartości pomiędzy
0,4 < | r | d" 0,6
 1 i 1
Umiarkowana zale\
ność liniowa
- 1 d" r d" 1
0,6 < | r | d" 0,8
" Znak współczynnika wskazuje
Silna zale\
ność liniowa
na kierunek zale\
ności.
| r | > 0,8
Bardzo silna zale\
ność liniowa
Liniowa funkcyjna zale\
ność między cechami
y
y
r =  1
r = 1
x
x
Własności współczynnika r
r = 0 oznacza brak zale\
ności r = 0 nie wyklucza występowania
liniowej. zale\
ności nieliniowej
Mo\
e to oznaczać brak korelacji
między cechami.
y
x
-6 -4 -2 0 2 4 6
Zale\
ność funkcyjna: y = x2
20
15
10
5
0
r = -0,9
" Współczynnik korelacji liniowej
Pearsona wynosi w przybli\
eniu -0,90
" Między temperaturą wody,
a ilością tlenu rozpuszczonego
zachodzi silna liniowa zale\
ność
o charakterze ujemnym
11 12 13 14 15 16
Temperatura wody
> plot(x=temperatura,y=tlen,xlab="Temperatura wody",
ylab="Tlen rozpuszczony", main="r = -0,9", pch=19)
> abline(lm(tlen~temperatura),col="red", lty=2)
9
8
7
Tlen rozpuszczony
6
5
Czy korelacja między temperaturą i ilością tlenu jest statystycznie istotna?
Hipoteza zerowa:
r = 0 (brak zale\
ności liniowej między temperaturą i ilością tlenu )
Hipoteza alternatywna:
r Ä„ 0 (istnieje zale\
ność liniowa między temperaturą i ilością tlenu,
przy czym nie przesÄ…dzamy o kierunku tej zale\
ności)
> cor.test(x=temperatura,y=tlen
method="pearson",
alternative="two.sided")
Pearson's product-moment correlation
data: temperatura and tlen
t = -5.9737, df = 8, p-value = 0.000333
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
Na poziomie istotności 0,05 nale\
y
-0.9772943 -0.6362455
odrzucić hipotezę zerową i przyjąć
sample estimates:
hipotezÄ™ alternatywnÄ….
cor
Zale\
ność jest statystycznie istotna.
-0.9038088
Czy ujemna korelacja między temperaturą i ilością tlenu jest statystycznie
istotna?
Hipoteza zerowa:
r e" 0 (brak zale\
ności liniowej lub zale\
ność dodatnia )
Hipoteza alternatywna:
r < 0 (istnieje ujemna zale\
ność liniowa między temperaturą i ilością tlenu)
> cor.test(x=temperatura,y=tlen,
method="pearson",
alternative="less")
Pearson's product-moment correlation
data: temperatura and tlen
t = -5.9737, df = 8, p-value = 0.0001665
alternative hypothesis: true correlation is less than 0
95 percent confidence interval:
-1.0000000 -0.7018518
Na poziomie istotności 0,05 nale\
y
sample estimates:
odrzucić hipotezę zerową i przyjąć
cor
hipotezÄ™ alternatywnÄ….
-0.9038088
Zale\
ność ujemna jest statystycznie
istotna.
> cor.test(x=temperatura,y=tlen,
method="pearson",
alternative="two.sided")
Pearson's product-moment correlation
data: temperatura and tlen
t = -5.9737, df = 8, p-value = 0.000333
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.9772943 -0.6362455
sample estimates:
cor
95% przedział ufności
-0.9038088
dla współczynnika korelacji ª [-0,98; -0,64]
Dla zupełnie ró\
nych układów
punktów na wykresie rozrzutu
W ka\
dym przypadku r = 0,82
współczynnik korelacji mo\
e
być taki sam, dlatego zawsze
konieczna jest analiza
graficzna zale\
ności między
cechami.
5 10 15 5 10 15
x1 x2
5 10 15 5 10 15
x3 x4
y1
y2
4
6
8
10
12
4
6
8
10
12
y3
y4
4
6
8
10
12
4
6
8
10
12
Nawet bardzo silna zale\
ność liniowa między cechami nie przesądza
o zale\
ności przyczynowo-skutkowej między tymi cechami.
Z
X
Y
Pozorny związek między X i Y
ANALIZA REGRESJI
" Czy istnieje metoda, która pozwala przewidzieć wartości
jednej cechy na podstawie znajomości wartości drugiej
cechy?
" Jeśli tak, to jaka jest dokładność tego przewidywania?
" Czy zale\
ność między cechami da się opisać analitycznie (wzorem)?
" Najprostszy rodzaj zale\
ności to zale\
ność liniowa.
Jak wyznaczyć równanie prostej, która najlepiej będzie
reprezentować analizowaną chmurę punktów na wykresie
rozproszenia?
Analiza regresji  terminologia
" Cecha (zmienna), której wartości
chcemy przewidywać
= zmienna objaśniana
(lub zmienna zale\
na).
Oznaczamy jÄ… Y.
" Cecha (zmienna), za pomocÄ…
której chcemy dokonać
przewidywania
= zmienna objaśniająca
zmienna objaśniająca X
(lub zmienna niezale\
na).
Oznaczamy jÄ… X.
zmienna objaśniana Y
" Załó\
my, \
e do analizowanej chmury punktów najlepiej pasuje prosta
o równaniu y = a0 + a1x.
X 1 2 3 4 5
Y 20 40 30 90 80
y = a0 + a1x
0 1 2 3 4 5 6 0 1 2 3 4 5 6
x x
Cel: wyznaczenie współczynników a0 oraz a1
y
y
0
20
40
60
80
100
0
20
40
60
80
100
Wartość y zmiennej objaśnianej przewidywana na podstawie prostej
w
i = a0 + a1xi
y = a0 + a1x dla x = xi to
90
86
80
69
Przewidywana
52
wartość zmiennej Y
dla x = 3
40
35
Obserwowana
30
wartość zmiennej Y
20
18
dla x = 3
0 1 2 3 4 5 6
x
Wartości przewidywane obliczono przy zało\
eniu, \
e y = 1 + 17x
100
80
60
y
40
20
0
90
Ró\
nica między wartością obserwowaną
86
i przewidywaną to wartość resztowa
80
(reszta, rezyduum)
69
Lp. xi yi w
i ei = yi - w
i
52
1. 1 20 18 2
40
2. 2 40 35 5
35
3. 3 30 52 -22
30
4. 4 90 69 21
20
18
5. 5 80 86 -6
Suma 15 260 260 0
0 1 2 3 4 5 6
Reszty traktujemy jako
x
błędy przewidywania
100
80
60
y
40
20
0
Metoda Najmniejszych Kwadratów
Jako prostą najlepiej reprezentującą chmurę punktów wybieramy
prostą y = a0 + a1x wyznaczoną przez takie współczynniki a1 i a0 ,
dla których suma kwadratów reszt jest najmniejsza.
n n
S(a0,a1) =
"(y - w
i )2 = "( yi - (a0 + a1xi ))2
i
i=1 i=1
Lp. xi yi w
i ei = yi - w
i
ei 2
1. 1 20 18 2 4
2. 2 40 35 5 25
Obliczenia wykonano
3. 3 30 52 -22 484
przy zało\
eniu,
4. 4 90 69 21 441
\
e a0 = 1, a1 = 17,
5. 5 80 86 -6 36
czyli y = 1 + 17x
Suma 15 260 260 0 990
" Problem:
wyznaczyć minimum funkcji
n n
S(a0,a1) =
"(y - w
i )2 = "( yi - (a0 + a1xi ))2
i
i=1 i=1
" Narzędzie do rozwiązania problemu:
rachunek ró\
niczkowy funkcji dwóch zmiennych.
" RozwiÄ…zanie problemu:
n
"(x - x)(yi - y)
i
i=1
a1 =
n
"(x - x)2
i
i=1
a0 = y - a1x
Prostą regresji opartą na Metodzie Najmniejszych Kwadratów
nazywamy prostą y = a0 + a1x, dla której wartość sumy
n n
S(a0,a1) =
"(y - w
i )2 = "( yi - (a0 + a1xi ))2
i
i=1 i=1
jest najmniejsza.
Niektóre własności prostej regresji
(x, y)
" Prosta regresji przechodzi przez punkt
" Suma reszt jest równa zero
n
"e = 0
i
i=1
Prosta regresji dla temperatury i ilości tlenu rozpuszczonego w wodzie
Wykres rozrzutu
Tlen
Temperatura rozpuszczony
Lp. (Ä…C) (mg/dm3)
1. 11,1 9,3
2. 11,3 9,6
3. 11,6 9,2
4. 11,7 7,3
5. 12,7 6,7
6. 13,8 6,0
7. 14,8 6,2
8. 15,2 6,2
9. 15,7 5,9
10. 16,3 4,3
11 12 13 14 15 16
Temperatura wody
9
8
7
Tlen rozpuszczony
6
5
Wpółczynniki prostej regresji wyznacza komenda lm (skrót od linear model)
> lm(tlen~temperatura)
Zmienna Zmienna
objaśniana objaśniająca
Call:
lm(formula = tlen ~ temperatura)
Coefficients:
(Intercept) temperatura
17.8224 -0.8012
Współczynnik a0 Współczynnik a1
Poszukiwana prosta regresji ma równanie y = 17,8 + ( 0,8)x,
czyli y = 17,8  0,8x
Przewidywana ilość tlenu = 17,8  0,8·temperatura
y =17,8  0,8x Przy wzroście temperatury
o 1 stopień Celsjusza
mo\
na oczekiwać spadku
ilości tlenu o 0,8 mg/dm3
11 12 13 14 15 16
Temperatura wody
9
8
7
Tlen rozpuszczony
6
5
Rozkład całkowitej zmienności zmiennej objaśnianej
" Całkowita suma kwadratów
n
" Zachodzi, \
e
SST =
"(y - y)2
i
i=1
SST = SSE + SSR
" Suma kwadratów błędów
n
SSE =
"(y - w
i )2
i
i=1
" Regresyjna suma kwadratów
n
SSR =
"( w
i - y)2
i=1
SST = SSE + SSR
Całkowita zmienność cechy Y =
= zmienność wyjaśniona przez zale\
ność liniową (względem cechy X)
+ zmienność nie wyjaśniona przez zale\
ność liniową (względem cechy X)
Współczynnik determinacji
SSR SSE
R2 = = 1-
SST SST
Współczynnik determinacji
SSR SSE
R2 = = 1-
SST SST
określa stopień, w jakim zale\
ność liniowa pomiędzy cechami Y i X
tłumaczy zmienność obserwowaną na wykresie rozproszenia.
Im R2 większy, tym skonstruowany model zale\
ności pomiędzy Y i X jest
lepszy.
Współczynnik determinacji jest powiązany ze współczynnikiem korelacji
liniowej Pearsona:
R2 = r2
a1 = 1.05 , r = 0.97 a1 = 1.2 , r = 0.78
Współczynnik determinacji = 0.94 Współczynnik determinacji = 0.6
Współczynnik a1 odpowiada
SST = 181.97 SST = 372.17
za szybkość zmiany zmiennej
SSR = 171.89 SSR = 223.85
SSE = 10.08 SSE = 148.32
objaśnianej.
2 4 6 8 10 2 4 6 8 10 Współczynnik korelacji opisuje
stopień skupiania się punktów
x x
wokół prostej regresji.
a1 = 0.52 , r = 0.92 a1 = 0.49 , r = 0.56
Współczynnik determinacji = 0.85 Współczynnik determinacji = 0.31
SST = 49.12 SST = 121.85
SSR = 41.62 SSR = 37.55
SSE = 7.5 SSE = 84.31
Znajomość jedynie
współczynnika korelacji nie
pozwala określić współczynnika
a1 i odwrotnie: znajomość a1 nie
pozwala określić współczynnika
korelacji.
2 4 6 8 10 2 4 6 8 10
x x
y
y
0
2
4
6
8
10
12
0
2
4
6
8
10
12
y
y
0
2
4
6
8
10
12
0
2
4
6
8
10
12
Model zale\
ności liniowej między zmiennymi
Wartość zmiennej objaśnianej =
= funkcja liniowa zmiennej objaśniającej + błąd losowy
Yi = a0 + a1 ·xi + ei , i = 1, 2, ..., n
" Oszacowaniami (estymatorami) nieznanych współczynników a0 i a1 są
odpowiednio a0 i a1. Współczynniki a0 i a1 nazywa się parametrami modelu.
" O błędzie losowym ei zakłada się, \
e jest zmienną losową o rozkładzie
normalnym z wartością oczekiwaną równą zero i nieznaną wariancją s2.
" Oszacowaniem nieznanej wariancji s2 jest błąd średniokwadratowy S2
n
2
"e
i
SSE
2
i=1
S = =
n - 2 n - 2
" Średni błąd szacunku modelu (odchylenie standardowe reszt) to
2
S = S
" Średni błąd szacunku modelu określa, jakie jest przeciętne odchylenie
obserwowanych wartości zmiennej objaśnianej od wartości
przewidywanych przez model (zale\
ności liniowej).
" Im S jest mniejsze, tym  dobroć dopasowania prostej regresji do
danych empirycznych (czyli punktów na wykresie rozrzutu) jest
lepsza.
" Współczynnik zmienności resztowej V
S
V = Å"100%
y
" Im współczynnik V mniejszy, tym model lepiej pasuje do danych.
a1 = 1.05 , r = 0.97 a1 = 1.2 , r = 0.78
Współczynnik determinacji = 0.94 Współczynnik determinacji = 0.6
S = 0.73 S = 2.79
V = 10.9 % V = 38.73 %
2 4 6 8 10 2 4 6 8 10
x x
a1 = 0.52 , r = 0.92 a1 = 0.49 , r = 0.56
Współczynnik determinacji = 0.85 Współczynnik determinacji = 0.31
S = 0.63 S = 2.11
V = 17.45 % V = 43.62 %
2 4 6 8 10 2 4 6 8 10
x x
y
y
0
2
4
6
8
10
12
0
2
4
6
8
10
12
y
y
0
2
4
6
8
10
12
0
2
4
6
8
10
12
" Poniewa\
oszacowania nieznanych współczynników a0 i a1
(nazywanych tak\
e parametrami modelu) dokonuje siÄ™ na podstawie
danych pochodzących z próby losowej, to przy ich szacowaniu równie\

popełnia się błędy. Wartości tych błędów mo\
na oszacować
wyznaczając średnie błędy szacunku parametrów:
1 x2
SE(a0) = S +
n
n
"(x - x)2
i
i=1
1
SE(a1) = S
n
"(x - x)2
i
i=1
> summary(lm(tlen~temperatura))
Call:
lm(formula = tlen ~ temperatura)
Residuals:
Min 1Q Median 3Q Max
-1.1481 -0.6898 0.3034 0.6316 0.8314
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.8224 1.8176 9.806 9.83e-06 ***
temperatura -0.8012 0.1341 -5.974 0.000333 ***
---
Signif. codes: 0  *** 0.001  ** 0.01  * 0.05  . 0.1   1
Residual standard error: 0.7977 on 8 degrees of freedom
Multiple R-squared: 0.8169, Adjusted R-squared: 0.794
F-statistic: 35.68 on 1 and 8 DF, p-value: 0.000333
Wybrane nazwy pól wartości funkcji lm()
coefficients współczynniki modelu
coefficients
coefficients
coefficients
residuals reszty modelu
residuals
residuals
residuals
fitted.values wartości teoretyczne (y  z daszkiem )
fitted.values
fitted.values
fitted.values
df.residual liczba stopni swobody dla reszt (liczba obserwacji
df.residual
df.residual
df.residual
minus liczba estymowanych parametrów modelu)
> model<-lm(tlen~temperatura)
#Reszty modelu
> model$residuals
1 2 3 4 5 6 7
0.3711669 0.8314111 0.6717775 -1.1481004 -0.9468792 -0.7655360 0.2356852
8 9 10
0.5561737 0.6567843 -0.4624830
> as.matrix(model$residuals)
[,1]
1 0.3711669
2 0.8314111
Residuals:
3 0.6717775
Min 1Q Median 3Q Max
4 -1.1481004
-1.1481 -0.6898 0.3034 0.6316 0.8314
5 -0.9468792
6 -0.7655360
7 0.2356852
Minimum, I kwartyl, mediana, III kwartyl i maksimum
8 0.5561737
dla reszt modelu
9 0.6567843
10 -0.4624830
> reszty<-model$residuals
> y_z_daszkiem<-model$fitted.values
> tabelka<-data.frame(empiryczne=tlen,
teoretyczne=y_z_daszkiem,
reszty=reszty)
> tabelka
empiryczne teoretyczne reszty
1 9.3 8.928833 0.3711669
2 9.6 8.768589 0.8314111
3 9.2 8.528223 0.6717775
4 7.3 8.448100 -1.1481004
5 6.7 7.646879 -0.9468792
6 6.0 6.765536 -0.7655360
7 6.2 5.964315 0.2356852
8 6.2 5.643826 0.5561737
9 5.9 5.243216 0.6567843
10 4.3 4.762483 -0.4624830
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.8224 1.8176 9.806 9.83e-06 ***
temperatura -0.8012 0.1341 -5.974 0.000333 ***
Oszacowania Średnie błędy
Wartości
p-wartości w teście
parametrów szacunku
statystyki
istotności parametrów
modelu parametrów
testowej t
---
Signif. codes: 0  *** 0.001  ** 0.01  * 0.05  . 0.1   1
Oznaczenia ułatwiające szybką ocenę poziomu p-wartości
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.8224 1.8176 9.806 9.83e-06 ***
temperatura -0.8012 0.1341 -5.974 0.000333 ***
Oszacowania Średnie błędy
parametrów szacunku
modelu parametrów
Oszacowanie parametru a0 (czyli a0 ) jest równe 17,8, zaś średni błąd
szacunku tego parametru jest równy 1,82
Oszacowanie parametru a1 (czyli a1 ) jest równe  0,8, zaś średni błąd
szacunku tego parametru jest równy 0,13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.8224 1.8176 9.806 9.83e-06 ***
temperatura -0.8012 0.1341 -5.974 0.000333 ***
p-wartości w teście
istotności parametrów
H0: a0 = 0 vs. H1: a0 +" 0
Poniewa\
p-wartość = 9.83e-06, to na poziomie istotności 0,05 nale\
y odrzucić H0
H0: a1 = 0 vs. H1: a1 +" 0
Poniewa\
p-wartość = 0.000333, to na poziomie istotności 0,05 nale\
y odrzucić H0
Wniosek:
oba parametry modelu są statystycznie istotne na poziomie istotności 0,05.
Residual standard error: 0.7977 on 8 degrees of freedom
Średni błąd szacunku modelu (pierwiastek z błędu średniokwadratowego) wynosi
0.7977 (przy 8 stopniach swobody)
Multiple R-squared: 0.8169, Adjusted R-squared: 0.794
Współczynnik determinacji R2 jest równy 0.8169. Obserwowana zmienność
zawartości tlenu rozpuszczonego w wodzie jest w około 82% wyjaśniona przez
liniowÄ… zale\
ność od temperatury wody.
F-statistic: 35.68 on 1 and 8 DF, p-value: 0.000333
W przypadku, gdy w modelu jest tylko jedna zmienna objaśniająca (w tym przypadku
jest niÄ… temperatura), wynik tzw. testu F jest to\
samy z wynikiem testu istotności
parametru a1)
Diagnostyka modelu  badanie własności reszt
Czy reszty mają rozkład normalny?
> shapiro.test(reszty)
Shapiro-Wilk normality test
data: reszty
W = 0.8691, p-value = 0.09767
Poniewa\
p-wartość > 0,05,
to na poziomie istotności 0,05
nie ma podstaw do odrzucenia hipotezy
o normalności rozkładu reszt
Przedziały ufności dla parametrów modelu
#95% przedziały ufności
> confint(object=model,level=0.95)
2.5 % 97.5 %
(Intercept) 13.631107 22.013669
temperatura -1.110514 -0.491928
#90% przedziały ufności
> confint(object=model,level=0.90)
5 % 95 %
(Intercept) 14.442564 21.202212
temperatura -1.050633 -0.551809
" Mo\
na przewidywać średnią wartość zmiennej tlen dla interesującej nas wartości
zmiennej temperatura (pamiętając o tym, \
e wartość ta powinna być z zakresu
zmienności zmiennej temperatura, tzn. w naszym przypadku dla wartości z zakresu
range(temperatura)).
" Np. gdyby interesowało nas przewidywanie średniej zawartości tlenu dla temperatur: 11.5,
12.0, 12.4, 15.7 stopni Celsjusza, to mo\
na u\
yć funkcji predict.lm, przy czym jako
jednego z argumentów tej funkcji nale\
y u\
yć ramki danych, w której jedna ze zmiennych ma
nazwę taką, jak zmienna objaśniająca modelu (w naszym przypadku  temperatura )
i zawiera wartości, dla których zamierzamy dokonać przewidywania.
> predict.lm(model,
data.frame(temperatura=c(11.5, 12.0, 12.4, 15.7 )))
1 2 3 4
8.608345 8.207734 7.887246 5.243216
" Wyznaczamy przedziały ufności dla powy\
szych przewidywań
(interval= confidence ) na wybranym przez siebie poziomie
ufności (level=)
> predict.lm(model,
data.frame(temperatura=c(11.5, 12.0, 12.4, 15.7 )),
interval="confidence",level=0.95)
fit lwr upr
1 8.608345 7.777042 9.439647
Górne (upper) granice
2 8.207734 7.478825 8.936643
przedziałów ufności
3 7.887246 7.225474 8.549017
4 5.243216 4.329047 6.157385
Wartości Dolne (lower) granice
przewidywane przedziałów ufności
dla poszczególnych
temperatur
Zawartość tlenu w zale\
ności od temperatury
(wykres rozrzutu wraz z 95% przedziałem ufności dla przewidywanych średnich)
11 12 13 14 15 16
temperatura
9
8
7
tlen
6
5
> iksy<- seq(from=min(temperatura),to=max(temperatura),by=0.1)
> dane.do.przewidywania<-data.frame(temperatura=iksy)
> średnie.wartości<-predict.lm(model,
dane.do.przewidywania,
interval="confidence",
level=0.95)
> plot(tlen~temperatura,
main="Zawartość tlenu w zale\
ności od temperatury \n (wykres
rozrzutu wraz z 95% przedziałem ufności dla przewidywanych
średnich)",
cex.main=0.8)
> abline(model)
> lines(x=iksy,y=średnie.wartości[,2],lty=2)
> lines(x=iksy,y=średnie.wartości[,3],lty=2)
" Mo\
na przewidywać przyszłą wartość zmiennej tlen.
" W odró\
nieniu od poprzedniej sytuacji estymujemy wartość
zmiennej losowej, a nie wartość oczekiwaną tej zmiennej.
" Z tego powodu, choć oszacowania punktowe będą identyczne, jak
w przypadku poprzednim, to przedziały ufności (nazywane w tym
przypadku przedziałami predykcji albo tolerancji) są dłu\
sze).
> predict.lm(model,
data.frame(temperatura=c(11.5, 12.0, 12.4, 15.7 )),
interval="prediction",level=0.95)
fit lwr upr
1 8.608345 6.589629 10.627060
2 8.207734 6.228983 10.186485
3 7.887246 5.932229 9.842262
4 5.243216 3.188988 7.297444
Zawartość tlenu w zale\
ności od temperatury
oraz prosta regresji wraz z 95% przedziałami ufności i tolerancji
prosta regresji
przedział ufności
przedział tolerancji
11 12 13 14 15 16
temperatura
9
8
7
tlen
6
5
> iksy<- seq(from=min(temperatura),to=max(temperatura),by=0.1)
> dane.do.przewidywania<-data.frame(temperatura=iksy)
> średnie.wartości<-predict.lm(model, dane.do.przewidywania,
interval="confidence",level=0.95)
> predykcja.wartości<-predict.lm(model, dane.do.przewidywania,
interval="prediction",level=0.95)
> plot(tlen~temperatura, main="Zawartość tlenu w zale\
ności od
temperatury \n oraz prosta regresji wraz z 95% przedziałami
ufności i tolerancji",cex.main=0.75)
> abline(model)
> lines(x=iksy,y=średnie.wartości[,2],lty=2,col="red",lwd=2)
> lines(x=iksy,y=średnie.wartości[,3],lty=2,col="red",lwd=2)
> lines(x=iksy,y=predykcja.wartości[,2],lty=3,col="blue",lwd=2)
> lines(x=iksy,y=predykcja.wartości[,3],lty=3, col="blue",lwd=2)
> legend("topright",legend=c("prosta regresji","przedział
Å‚
Å‚
Å‚
ufnoÅ› Å‚
ści","przedział tolerancji"),
Å› Å‚
Å› Å‚
text.col= c("black","red","blue"), lty=c(1,2,3),
lwd=c(1,2,2),col=c("black","red","blue"))
Regresja wielokrotna (wieloraka)
" Rzadko zdarza się sytuacja, gdy zmienna objaśniana zale\
y tylko od
jednej zmiennej objaśniającej.
" Załó\
my, \
e mamy k zmiennych objaśniających
y1 x11 K x1k
M M M
yn xn1 K xnk
Yi = Ä…0 + Ä…1xi1 + ... + Ä…k xik + µi
i = 1, 2,..., n
1 x11 K x1k
îÅ‚ Å‚Å‚
y1
îÅ‚ Å‚Å‚
Ä…0
îÅ‚ Å‚Å‚
µ1
îÅ‚ Å‚Å‚
ïÅ‚1
ïÅ‚y śł
ïÅ‚Ä… śł
ïÅ‚µ śł
x21 K x2k śł
2 1
ïÅ‚ śł 2
ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł
X = Ä… =
y =
µ =
ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł M ïÅ‚ śł
M M O M M
ïÅ‚ śł
ïÅ‚ śł
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚Ä…k ûÅ‚ ðÅ‚µn ûÅ‚
xnk K xnk ûÅ‚
n
ðÅ‚y ûÅ‚
ðÅ‚1
Yi = Ä…0 + Ä…1xi1 + ... + Ä…k xik + µi
Zapis skalarny
i = 1, 2,..., n
y = XÄ… + µ
Zapis macierzowy
Metoda najmniejszych kwadratów w przypadku wielowymiarowym:
poszukujemy wektora a, gdzie
a0
îÅ‚ Å‚Å‚
ïÅ‚a śł
1
ïÅ‚ śł
a =
ïÅ‚ śł
M
ïÅ‚ śł
ðÅ‚ak ûÅ‚
takiego, \
e
n n
S(a0,a1,..., ak ) =
"( yi - w
i )2 = "( yi - (a0 + a1xi1 + ... + ak xik ))2 = (y - Xa)T (y - Xa)
i=1 i=1
przyjmie wartość minimalną.
Poszukiwany wektor a jest określony wzorem:
a = (XTX)-1 XTy
> airquality
Ozone Solar.R Wind Temp Month Day
1 41 190 7.4 67 5 1
2 36 118 8.0 72 5 2
3 12 149 12.6 74 5 3
4 18 313 11.5 62 5 4
5 NA NA 14.3 56 5 5
6 28 NA 14.9 66 5 6
7 23 299 8.6 65 5 7
8 19 99 13.8 59 5 8
9 8 19 20.1 61 5 9
10 NA 194 8.6 69 5 10
.
.
.
> str(airquality)
'data.frame': 153 obs. of 6 variables:
$ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
$ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 5 6 7 8 9 10 ...
NA  not available (brak danych)
> help(airquality)
# wykluczenie braków danych
> ozon<-airquality[complete.cases(airquality),]
> str(ozon)
'data.frame': 111 obs. of 6 variables:
$ Ozone : int 41 36 12 18 23 19 8 16 11 14 ...
$ Solar.R: int 190 118 149 313 299 99 19 256 290 274 ...
$ Wind : num 7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...
$ Temp : int 67 72 74 62 65 59 61 69 66 68 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 7 8 9 12 13 14 ...
> ozon<-read.csv2("F:/Moje dokumenty/dane/ozon_dane.csv")
> plot(ozon)
0 150 300 60 80 0 10 20 30
Ozone
Solar.R
Wind
Temp
Month
Day
0 50 150 5 15 5 6 7 8 9
0
50
150
0
150
300
5
15
60
80
5
6
7
8
9
0
10 20 30
> names(ozon)
[1] "Ozone" "Solar.R" "Wind" "Temp" "Month"
"Day"
> dane<-ozon[,1:4]
> names(dane)
[1] "Ozone" "Solar.R" "Wind" "Temp"
#wyznaczamy macierz korelacji pomiędzy zmiennymi
> cor(dane)
Ozone Solar.R Wind Temp
Ozone 1.0000000 0.3483417 -0.6124966 0.6985414
Solar.R 0.3483417 1.0000000 -0.1271835 0.2940876
Wind -0.6124966 -0.1271835 1.0000000 -0.4971897
Temp 0.6985414 0.2940876 -0.4971897 1.0000000
Chcemy zbadać, czy zawartość ozonu w powietrzu mo\
na
przewidywać w oparciu o zmienne "Solar.R","Wind","Temp"
Konstruujemy model liniowy postaci
Ozone = Ä…0 + Ä…1Solar.R + Ä…2Wind + Ä…3Temp + µ
Przyłączamy ramkę dane komendą attach(). Dzięki temu będziemy
mieć bezpośredni dostęp do zmiennych z tej ramki poprzez nazwy tych
zmiennych, co ułatwia pracę w programie.
> attach(dane)
> #konstruujemy model liniowy
> model<-lm(Ozone~Solar.R+Wind+Temp)
> model
Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)
Coefficients:
(Intercept) Solar.R Wind Temp
-64.34208 0.05982 -3.33359 1.65209
Ozone = Ä…0 + Ä…1Solar.R + Ä…2Wind + Ä…3Temp + µ
Ozone = -64,34 + 0,06Solar.R - 3,33Wind +1,65Temp
> summary(model)
Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)
Residuals:
Min 1Q Median 3Q Max
-40.485 -14.219 -3.551 10.097 95.619
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -64.34208 23.05472 -2.791 0.00623 **
Solar.R 0.05982 0.02319 2.580 0.01124 *
Wind -3.33359 0.65441 -5.094 1.52e-06 ***
Temp 1.65209 0.25353 6.516 2.42e-09 ***
---
Signif. codes: 0  *** 0.001  ** 0.01  * 0.05  . 0.1   1
Residual standard error: 21.18 on 107 degrees of freedom
Multiple R-squared: 0.6059, Adjusted R-squared: 0.5948
F-statistic: 54.83 on 3 and 107 DF, p-value: < 2.2e-16
Diagnostyka modelu  badanie własności reszt
Czy reszty mają rozkład normalny?
> reszty<-model$resid
> shapiro.test(reszty)
Shapiro-Wilk normality test
data: reszty
W = 0.9171, p-value = 3.618e-06
Poniewa\
p-wartość < 0,05,
to na poziomie istotności 0,05
nale\
y odrzucić hipotezę
o normalności rozkładu reszt
> par(mfrow=c(2,2))
Histogram of reszty
> hist(reszty)
> boxplot(reszty)
> qqnorm(reszty)
> qqline(reszty)
> par(mfrow=c(1,1))
-50 0 50 100
reszty
Normal Q-Q Plot
-2 -1 0 1 2
Theoretical Quantiles
Frequency
0
10
30
50
-40
0
40
80
Sample Quantiles
-40
0
40
80
" Budowanie modelu zale\
ności między zmiennymi odbywa się często metodą
prób i błędów.
" W kolejnych etapach do modelu dodaje siÄ™ (lub usuwa) zmienne lub transformacje
(przekształcenia) zmiennych. Opieramy się tu niekiedy na merytorycznej znajomości
modelowanego zjawiska.
" Próbujemy poprawić model, tak aby był lepiej dopasowany do danych
empirycznych, a tak\
e aby reszty modelu miały rozkład normalny.
" Pierwsza próba poprawienia dopasowania modelu do danych empirycznych:
dodamy do modelu kwadraty wartości poszczególnych zmiennych.
> model.2<-
lm(Ozone~Solar.R+I(Solar.R^2)+Wind+I(Wind^2)+Temp+I(Temp^2))
I() to funkcja identycznościowa
> summary(model.2)
Call:
lm(formula = Ozone ~ Solar.R + I(Solar.R^2) + Wind + I(Wind^2) + Temp +
I(Temp^2))
Residuals:
Min 1Q Median 3Q Max
-48.698 -10.926 -3.786 9.201 79.932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.984e+02 1.016e+02 2.937 0.00408 **
Solar.R 1.347e-01 8.806e-02 1.529 0.12919
I(Solar.R^2) -2.044e-04 2.549e-04 -0.802 0.42444
Wind -1.334e+01 2.308e+00 -5.783 7.79e-08 ***
I(Wind^2) 4.642e-01 1.010e-01 4.594 1.22e-05 ***
Temp -6.585e+00 2.742e+00 -2.402 0.01808 *
I(Temp^2) 5.223e-02 1.786e-02 2.925 0.00423 **
---
Signif. codes: 0  *** 0.001  ** 0.01  * 0.05  . 0.1   1
Residual standard error: 18.3 on 104 degrees of freedom
Multiple R-squared: 0.7141, Adjusted R-squared: 0.6976
F-statistic: 43.29 on 6 and 104 DF, p-value: < 2.2e-16
" model.2 ma wy\
szy współczynnik determinacji R2
" Zmienne Solar.R oraz I(Solar.R^2)sÄ… nieistotne statystycznie.
" Mo\
na spróbować usunąć najpierw zmienną I(Solar.R^2)
> model.3<-lm(Ozone~Solar.R+Wind+I(Wind^2)+Temp+I(Temp^2))
> summary(model.3)
Call:
lm(formula = Ozone ~ Solar.R + Wind + I(Wind^2) + Temp + I(Temp^2))
Residuals:
Min 1Q Median 3Q Max
-48.017 -10.810 -4.144 8.120 80.125
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 291.09564 101.00727 2.882 0.00479 **
Solar.R 0.06593 0.02007 3.285 0.00139 **
Wind -13.37647 2.30330 -5.808 6.83e-08 ***
I(Wind^2) 0.46372 0.10087 4.597 1.20e-05 ***
Temp -6.34116 2.72014 -2.331 0.02165 *
I(Temp^2) 0.05104 0.01777 2.873 0.00492 **
---
Signif. codes: 0  *** 0.001  ** 0.01  * 0.05  . 0.1   1
Residual standard error: 18.27 on 105 degrees of freedom
Multiple R-squared: 0.7123, Adjusted R-squared: 0.6986
F-statistic: 51.99 on 5 and 105 DF, p-value: < 2.2e-16