Rozdział 8 BADANIE WSPÓŁZALEŻNOŚCI ZMIENNYCH LOSOWYCH
Dane są dwie zmienne losowe X i Y oraz dwuwymiarowa prosta próba losowa {xi,yi} = {(x1,y1), (x2,y2), ..., (xn,yn)}. Interesują nas dwa problemy:
(A) czy zmienne X i Y są niezależne?
(B) jeśli zmienne X, Y są zależne, to jaki jest kształt tej zależności?
Rozwiązania tych zadań muszą być oczywiście uzyskane na podstawie tej próby losowej. Najlepiej jest, aby próba pochodziła z dwuwymiarowej populacji normalnej, gdyż do takiego przypadku odnosi się większość teoretycznych rozwiązań obu powyższych problemów. Również w niniejszym tekście założenie to jest przyjęte – w przeciwnym przypadku będzie to wyraźnie zaznaczone.
8.1 TESTOWANIE NIEZALEŻNOŚCI ZMIENNYCH X I Y
Problem (A) sprowadza się do testowania hipotezy statystycznej H0(zmienne X i Y są niezależne).
Istnieje szereg sposobów testowania tej hipotezy.
1. Test niezależności χ2 zmiennych X i Y: H0[FX,Y(x,y)/FX(x)@FY(y)];
2. Testowanie istotności
współczynnika korelacji ρ zmiennych X i Y: H0(ρ=0);
3. Test istotności współczynnika α regresji zmiennej Y względem zmiennej X: H0(α=0)
4. Testowanie istotności współczynnika korelacji rangowej zmiennych X i Y (τ Kendalla lub ρ Spearmanna)
Rys. 8.1. Przypomnienie: niezależność zmiennych X,Y pociąga za sobą ich nieskorelowanie ale (niestety) nie na odwrót.
nieskorelowanie
( ρ ρ ρ ρ = 0) skorelowanie
( ρ ρ ρ ρ ≠ ≠ ≠ ≠ 0)
Najczęściej zamieniamy słowo "niezależność" na "nieskorelowanie" (jakkolwiek nie są to synonimy, zob. rys. 8.1) i hipoteza H0 przybiera kształt H0(ρ=0), gdzie ρ jest współczynnikiem korelacji liniowej badanych zmiennych. Dalsze postępowanie jest już znane – analogiczne do opisanego w rozdziale dotyczącym testowania hipotez parametrycznych.
8.2 ANALIZA REGRESJI
8.2.1 Sformułowanie problemu i metody jego rozwiązania
Sformułowanie problemu. Mamy powody przypuszczać, że zmienne X i Y są zależne i interesuje nas kształt tej zależności wyrażony przez funkcję regresji freg(x)
E( | )Y x = freg( ; , , ,...) ?xα β γ = (8.1) gdzie parametry α, β, γ,..., rzeczywistej funkcji regresji freg(x) nazywane są rzeczywistymi współczynnikami regresji. Jeśli funkcja freg(x) jest liniowa regresję nazywamy liniową – w wypadku przeciwnym regresja jest krzywoliniowa lub nieliniowa.
Informację podaną w równaniu zapisuje się dość często nieco inaczej, wyróżniając pewne wielkości, na jakie można rozłożyć zmienną warunkową Y*x oznaczoną w poniższym równaniu symbolem Y|x.
|x E( | ) x
Y = Y x + ∈ (8.2)
Wielkość 0x jest zmienną losową nazywaną resztą zmiennej Y względem średniej warunkowej E(Y*x) (używane są też terminy zmienna resztowa lub zmienna reszt), a wyrażenie (8.2) wraz z ewentualnymi założeniami nazywane jest modelem regresji zmiennej Y względem zmiennej X. Wielkość ∈x – pomimo iż jest zmienną losową – jest przeważnie oznaczana małą literą; tak będzie również w dalszym ciągu niniejszego tekstu.
Z definicji (8.2) reszty∈x mamy, że
E∈ =x 0 (8.3)
dla każdego x z zakresu zmienności zmiennej X. Tak więc zmienność zmiennej Y została rozłożona na wielkość, która jest funkcją (tj. E(Y*x)) oraz wielkość, która jest zmienną losową (tj. ∈x). Można też spotkać określenia składnik systematyczny i składnik losowy modelu. Przy ustalonej wartości x zmiennej X zmienność zmiennej Y|x jest zawarta wyłącznie w zmiennej resztowej:
varY|x =var∈ x (8.4)
bo E(Y*x) jest tutaj stałą. Korzystając z oznaczeń (8.2) i (8.4) możemy przedstawić rozdział wariancji zmiennej Y poprzez regresję w następującej postaci
|
var var E( | ) E var ( | ) var E X E var X
Y Y X Y X
Y
= +
= + ∈ (8.5)
Wielkość Evar∈x nazywa się wariancją resztową. Poszukiwanie kształtu funkcji regresji może oznaczać dwa różne zadania:
(a1) znaleźć postać funkcji freg(x: α,β,γ,...) i wartości jej parametrów α, β, γ,..., albo
(a2) postać funkcji freg(x; α,β,γ,..., ) jest znana – znaleźć wartości jej parametrów α, β, γ,...,
Na ogół postać funkcji regresji jest postulowana – np. funkcja liniowa, funkcja kwadratowa, eksponencjalna (czyli typu y(x)=aebx), potęgowa itp. – i zadanie regresji to zadanie (a2).
Po rozwiązaniu zadania (a1) lub (a2) – a więc gdy znany jest kształt funkcji regresji – zachodzi naturalne pytanie:
(b) jaka jest jakość otrzymanego rozwiązania problemu (a1) lub (a2)?
Zadania (a1), (a2) i (b) to zadania estymacji krzywej regresji. Jeśli zostały one wykonane, możemy krzywą tę wykorzystywać do zadań predykcji wartości zmiennej losowej Y
Dwuwymiarowa próba losowa jest podstawową informacją o regresji. Podstawową informacją o funkcji regresji i jej współczynnikach jest n łącznych pomiarów wielkości X i Y.
Pomiary te tworzą dwuwymiarową prostą próbę losową ZXY = {xi,yi} = {(x1,y1), (x2,y2), ..., (xn,yn)}, np. taką, jak pokazana na rys. 8.2. Próba ta – wraz z pewnymi założeniami – służy do rozwiązania problemów (a1), (a2) i (b).
W języku statystyki zadanie regresji można sformułować następująco: na podstawie danej próby losowej ZXY
znaleźć estymator
ˆE( | )Y x = fˆreg( ; , , ,...)x a b c ozn= yśr( ; , , ,...)x a b c (8.6) krzywej regresji (8.1) i ocenić jego jakość przez podanie na ile otrzymana (empiryczna) krzywa regresji yśr(x) (8.6) różni się od rzeczywistej krzywej regresji freg(x). Musi więc zostać podane kryterium, na podstawie którego będziemy poszukiwać empirycznej krzywej regresji. Najczęściej do tego celu wykorzystuje się metodę najmniejszych kwadratów.
Rys. 8.2. 50-elementowa dwuwymiarowa prosta próba losowa wylosowana z pewnej populacji
Empiryczna postać modelu (8.2), która posłuży do rozwiązywania problemów estymacji w regresji zmiennej Y względem zmiennej X, jest następująca:
( ; , , ,...) ˆ
i śr i i
y = y x a b c + ∈ (8.7)
gdzie yi jest wartością zmiennej Y*xi, a ∈ˆi jest resztą zmiennej Y*xi, i=1,2,...,n, (rys.
8.3). Zakładając dalej, że sama postać funkcji yśr(x) jest znana (a więc rozwiązujemy zadanie (a2) i, ewentualnie, (b)), a nie znane są tylko wartości parametrów α, ß, γ,..., zadanie (a2) regresji sprowadza się do znalezienia wartości a, b, c,... estymatorów współczynników ", ß, γ,..., funkcji .
Znalezienie współczynników funkcji regresji nie powinno kończyć zadania estymacji, gdyż istotną sprawą jest ocena jakości uzyskanego rozwiązania. Aby to móc uczynić, zauważmy, że otrzymane wartości a, b, c,... estymatorów
współczynników α, β, γ, ..., nie są jedynymi możliwymi wartościami parametrów.
Gdyby z tej samej populacji (X,Y) wylosować drugą, trzecią, ..., n-elementową próbę losową, to za każdym razem obliczone wartości a, b, c,... byłyby inne. Fakt ten zaznaczamy pisząc
( ; , , ,...) ˆ
i śr i i
Y =Y x A B C + ∈ (8.8)
Wszystkie trzy składniki powyższego równania są dla ustalonej wartości xi zmiennymi losowymi. Jeśli potrafimy obliczyć np. varYśr(x), to potrafimy powiedzieć o jakości uzyskanej regresji.
Metoda najmniejszych kwadratów jako metoda estymacji współczyn- ników regresji. Istota tej metody polega na minimalizacji sumy kwadratów różnic wartości pomierzonych yi i modelowych yśr(xi) we wszystkich punktach xi, i=1,2,...,n:
[ ]
2 21 1
( , , ,...) ( ; , , ,...) min
n n
i śr i i
i i
F a b c y y x a b c
= =
=
∑
− =∑
∈ = (8.9)czyli znalezieniu takich wartości a, b, c, ..., dla których funkcja F(a,b,c,...) przyjmuje najmniejszą z możliwych wartości. Jeśli F(a,b,c,...)/n nazwiemy średnim błędem kwadratowym, to opisaną metodę możemy nazwać metodą minimalizacji średniego błędu kwadratowego.
Rys. 8.3. Metoda najmniejszych kwadratów polega na znalezieniu takiej krzywej yśr(x),która minimalizuje sumę kwadratów
∑
∈2i reszt ∈ ixi
yśr (xi) yi
∈i
Ponieważ F(a,b,c,...) jest sumą kwadratów, jest więc funkcją nieujemną. Wynika stąd, że posiada ona (co najmniej jedno) minimum. Oznacza to, że pierwsze pochodne funkcji F zerują się w punkcie a, b, c, ... – wystarczy więc znaleźć rozwiązania następującego układu równań
[ ]
[ ]
2
1
2
1
( , , ,...)
( ; , , ,...) 0 ( , , ,...)
( ; , , ,...) 0 ...
n
i śr i
i n
i śr i
i
F a b c
y y x a b c
a a
F a b c
y y x a b c
b b
=
=
∂ ∂
= − =
∂ ∂
∂ ∂
= − =
∂ ∂
∑
∑
(8.10)Układ ten nie zawsze daje rozwiązania w sposób analityczny i trzeba wtedy uzyskiwać je za pomocą metod numerycznych.
W kolejnych rozdziałach rozpatrywać będziemy szczegółowo szczególny przypadek regresji – regresję liniową – który posłuży nam jako przypadek wzorcowy ujawniający kilka ważnych zagadnień wchodzących w skład analizy regresji.
8.2.2 Regresja liniowa zmiennej Y względem zmiennej X
Klasyczny model regresji liniowej. Zakładamy, że regresja Y względem X jest regresją liniową i ma postać:
E( | )Y x =αx+β (8.11)
gdzie α i β są rzeczywistymi (nieznanymi) współczynnikami prostej regresji (parametrami populacji (X,Y)). Rozpatrywany model (8.2) zależności zmiennych X i Y ma więc konkretną postać
|x E( | ) x x
Y = Y x + ∈ =αx+ + ∈ β (8.12)
Dana jest n-elementowa dwuwymiarowa próba losowa {xi,yi} = {(x1,y1), (x2,y2), ..., (xn,yn)}.
E( | )
i i i i i
Y = Y x + ∈ =αx + + ∈β (8.13)
gdzie Yi oznacza zmienną warunkową Y|xi.
Na razie nic nie zakładamy o postaci rozkładów prawdopodobieństwa zmiennych występujących w równaniu (8.13). (Nie musimy tego robić dopóki nie jest nam potrzebne prawdopodobieństwo związane z regresją). Musimy natomiast przyjąć pewne założenia o własnościach zmiennych losowych Yi. Zwykle przyjmuje się, że
2
(i) E
(ii) var var
(iii) cov( , ) 0 dla
i i
i i
i j
Y x
Y const
Y Y i j
α β
σ
= +
= ∈ = =
= ≠
(8.14)
Model (8.13) wraz założeniami (8.14) nazywany jest klasycznym modelem regresji liniowej.
Założenie (8.14)(i) w jest niezbędne, gdyż oznacza, że wszystkie wartości próby {xi, yi} pochodzą z populacji o tej samej (prostoliniowej) średniej warunkowej.
Założenie (8.14)(ii) wyraża bardzo wygodny dla obliczeń warunek, że zmienność warunkowa zmiennej losowej Y nie zależy od wartości x zmiennej X, co powoduje, że wariancja resztowa 2
Yr
σ zmiennej Y wyraża się bardzo prosto:
2 E var( | ) E var E 2 2
r X
Y Y X
σ = = ∈ = σ =σ (8.15)
Warunek (8.14)(iii) w jest zapisem o braku korelacji pomiędzy zmiennymi warunkowymi Yi i Yj i jest związany ze sposobem, w jaki pobierana jest próba losowa.
Oczywiście, na mocy własności kowariancji, mamy cov(Yi,Yi) = varYi.
Sformułowanie problemu estymacji. Zadanie (a2) regresji wygląda nastę- pująco: na podstawie danej próby losowej ZXY znaleźć estymator prostej regresji (8.11) w postaci
śr( )
y x = ⋅ +a x b (8.16)
W tym celu należy obliczyć estymatory a i b współczynników α i β równania (8.11).
Przykład 8.1. Korelacja lipcowych sum opadów w Rabce i Krynicy. Analiza wstępna. W tabeli 8.1 dany jest zbiór sum opadu deszczu (w mm), jakie zanotowano w lipcu w Rabce i Krynicy w latach 1961-1990. Zadanie Przeprowadzić wstępną analizę graficzną.
Tabela 8.1. Lipcowe sumy opadów w Rabce (x) i Krynicy (y)
rok 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 x [mm] 133.0 157.5 88.7 68.4 117.2 152.7 31.2 210.5 81.8 232.4 y [mm] 97.4 136.1 79.0 66.7 174.2 88.8 74.8 182.2 58.1 213.1 rok 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 x [mm] 81.0 100.8 276.3 191.1 194.1 73.4 125.3 64.3 139.2 222.6 y [mm] 81.9 63.7 126.4 144.2 132.9 75.2 100.1 60.9 76.7 211.5 rok 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 x [mm] 107.8 129.0 116.0 145.3 94.4 62.6 49.5 120.6 110.8 81.8 y [mm] 140.2 147.6 132.1 101.3 100.5 51.5 48.9 78.2 80.7 73.3
Rys. 8.4 zawiera dane z tabeli 8.1. Wy- daje się, że sumy opadów w obu stacjach są dość silnie skorelowane dodatnio. Tego zresztą (tzn. dodatniej korelacji) należa- łoby oczekiwać ze względu warunki klimatyczne rejonu, w którym stacje te się znajdują i przedział czasu, w jakim sumowana jest wysokość opadu. Należy też zwrócić uwagę na fakt, że w miarę wzrostu wartości sumy opadu wzrasta rozrzut wartości (xi,yi). Może to oznaczać, że nie tylko wartość średnia ale i wariancja warunkowa zmiennej y zależy od wartości zmiennej X.
Rozwiązanie zadania regresji będzie polegało na:
(1) znalezieniu (estymatorów) współczynników regresji oraz (2) ocenie jej jakości, co można uczynić
(2a) podając wariancje współczynników regresji
(2b) podając przedziały ufności dla estymowanych współczynników regresji lub (2c) wykreślając tzw. (1-α)100% obszar ufności dla prostej regresji.
Wyznaczenie współczynników regresji. Celem wyznaczenia wartości współczynników a i b posłużymy się warunkiem (8.9), gdzie funkcja freg jest funkcją liniową:
[ ]
21
( , ) n i ( i ) min
i
F a b y ax b
=
=
∑
− + = (8.17)Dalszy ciąg obliczeń sprowadza się do obliczenia pochodnych cząstkowych funkcji F po a i b i przyrównaniu ich do zera. Dostaniemy
[ ]
[ ]
1
1
( , )
2 ( ) ( ) 0
( , )
2 ( ) ( 1) 0
n
i i i
i n
i i
i
F a b
y ax b x
a F a b
y ax b b
=
=
∂ = − + − =
∂
∂ = − + − =
∂
∑
∑
(8.18)Dzieląc oba równania przez -2n i używając symbolu "kreska ponad" jako znaku średniej arytmetycznej, można napisać w bardziej zwartej formie:
2 0
0 yx a x bx
y ax
− − =
− = (8.19)
Rys. 8.4. „Chmura punktów” – realizacje (xi,yi),i=1,2,...,30, zmiennych losowych X i Y – lipcowych sum opadu w Krynicy i Rabce.
skąd po krótkich rachunkach otrzymujemy równania definiujące empiryczne współczynniki regresji
( )
22
yx x y a
x x
b y ax
= − ⋅
−
= −
(8.20)
W liczniku wzoru na wartość a rozpoznajemy empiryczny (tj. z próby) odpowiednik znanej nam już kowariancji; w mianowniku znajduje się wariancja sX2 z próby zmiennej X. Wstawiając do wzoru (8.20) na a wzór (7.16) definiujący współczynnik korelacji z próby rXY, dostajemy inną postać współczynnika kierunkowego regresji liniowej:
Y XY
X
a r s
= s (8.21)
Przykład 8.2. Obliczanie prostej regresji lipcowych sum opadów w Krynicy względem lipcowych sum opadów w Rabce.
Zadanie: Na podstawie danych zawartych w tabeli 8.1 obliczyć współczynniki regresji (8.20) i wykreślić prostą regresji yśr(x).
Rozwiązanie. Do znalezienia współ- czynników regresji a i b (do obliczenia wartości a użyjemy równania (8.21)) wystarczy obliczyć pięć momentów:
dwie wartości średnie x i y, dwa od- chylenia standardowe sX i sY oraz współczynnik korelacji r zmiennych X i Y. Obliczono, że:
x= 125.31 mm, y = 106.61 mm, sX = 57.74 mm, sY = 45.11 mm oraz r = 0.754.
Podstawiając te wartości do wzoru (8.21) i drugiego ze wzorów (8.20), dostajemy: a = 0.589 mm/mm, b = 32.83 mm. Otrzymana prosta regresji wykreślona została na rys. 8.5.
Przypominamy, że prosta ta jest oceną
(estymatorem) prawdziwej średniej warunkowej E(Y*x).
Niepewność uzyskanej oceny prostej regresji. Mając obliczone a i b możemy już narysować estymowaną prostą regresji yśr(x), pokazaną przykładowo na rys. 8.5. Prosta ta jest odpowiednikiem wartości estymatora w estymacji punktowej – jest to jedna z możliwych prostych regresji , empirycznych ocen rzeczywistej prostej
Rys. 8.5. Estymowana prosta regresji lipcowych sum opadu w Krynicy względem lipcowych sum opadu w Rabce, yśr(x) = ax+b, o współczynnikach obliczonych wzorem (8.20).
y(x) = αx + β. Zachodzi więc naturalne pytanie: w jakim stopniu uzyskana realizacja yśr(x) = ax + b rzeczywistej prostej y(x) = αx + β odzwierciedla tę prostą. Odpowiedź na to pytanie będziemy wyrażać liczbowo w dwojaki sposób: przez wariancję resztową z próby (i ewentualnie przez związany z nią współczynnik determinacji R2) oraz przez obszar ufności prostej regresji.
8.2.3 Ocena jakości estymacji : wariancja resztowa z próby i jej związek z jakością predykcji wartości y zmiennej Y
Wariancja resztowa. Jeśli do równania wstawimy obliczone wartości współ- czynników a i b minimalizujących sumę kwadratów odchyleń od prostej yśr(x)=ax+b, otrzymamy empiryczny odpowiednik (czyli estymator) wariancji resztowej E[var(Y|x)] (8.15):
2 2
[ ]
21 1
1 1
E var( | ) rY n i n i ( i )
i i
Y X s y ax b
n = n =
= =
∑
∈ =∑
− + (8.22)Wariancja resztowa (8.22) może być interpretowana jako średni kwadrat odchylenia zmiennej Y od uzyskanej średniej warunkowej, przypadający na jedną obserwację.
Nie jest to więc bezpośrednia ocena niepewności uzyskanej prostej regresji yśr, a ocena niepewności związanej z pojedynczą obserwacją y zmiennej Y. Inaczej możemy powiedzieć, że wariancja resztowa jest miarą rozproszenia zmiennej Y wokół prostej regresji.
Estymator (8.22) jest estymatorem obciążonym, tzn. że dla ustalonej liczebności n próby zachodzi nierówność
ESrY2 ≠E var ( | )Y X (8.23) Nieobciążonym estymatorem wariancji resztowej jest
[ ]
22
1
1 ( )
2
n
rY i i
i
s y ax b
n =
= − +
−
∑
(8.24)
Pierwiastek z wariancji (8.22) lub (8.24), srY, jest odchyleniem resztowym. Obie te wielkości (tzn. wariancja resztowa i odchylenie resztowe) pełnią rolę miary rozrzutu punktów eksperymentalnych(wartości z próby) wokół prostej regresji, a tym samym miary jakości dopasowania prostej do punktów eksperymentalnych (do wartości z próby): im są mniejsze, tym dopasowanie jest lepsze.
Przykład 8.3. Odchylenie resztowe w regresji lipcowych sum opadów w Krynicy względem lipcowych sum opadów w Rabce.
Zadanie. Obliczyć wartość odchylenia resztowego (8.22) i wykreślić proste yśr(x)+srY i yśr(x)- srY.
Rozwiązanie. Wstawiając wartości (x,y) z tabeli 8.1 i wartości współczynników a i b regresji yśr(x) z przykładu 8.2: a = 0.589 mm/mm, b =
32.83 mm, do wzoru (8.22) znajdujemy, że wartość odchylenia resztowego srY = 29.7 mm (rys. 8.6 ).
Interpretacja. Fakt ten można praktycznie wykorzystać np. w ten sposób: jeśli suma opadu w lipcu w Rabce (np. w roku 2001, chociaż rok nie ma znaczenia) jest znana i wynosi x mm i znana jest postać prostej regresji yśr(x), to możemy powiedzieć, że suma opadu w lipcu w Krynicy (np. w roku 2001, chociaż rok nie ma znaczenia) wynosi yśr(x) ± srY. Na przykład, jeśli x=200 mm, to z rys. 8.6 znajdujemy, że yśr(x) = 150 mm, a niepewność tej predykcji wynosi srY
= 29.7 mm. 30 mm.
Szczególna postać wzoru na wariancję resztową w regresji prostolinio- wej. Podstawiając do (8.22) wartości a z (8.21) i b z (8.20), dostajemy, że
[ ] [ ]
( )
2 2
2
2 2 2
2
2 2
2 2
( ) ( ) ( )
( ) 2 ( )( ) ( )
2 1
rY
Y Y
Y X Y X
X X
Y
s y ax b y y a x x
y y a x x y y a x x
s s
s r rs s s
s s
s r
= − + = − − −
= − − − − + −
= − +
= −
(8.25)
Jest to postać już znana dla przypadku teoretycznego (3.107)
Różnica pomiędzy wariancją resztową a wariancją zmiennej Y pokazuje ile zmien- ności zmiennej Y zostało, jak się to mówi, wyjaśnione wprowadzeniem prostej re- gresji, czyli o ile poprawi się (średnio!) dokładność przewidywania wartości zmiennej Y, gdy znać będziemy wartość x zmiennej X. Fakt ten ilustruje poniższy przykład.
Przykład 8.4. Różnica pomiędzy odchyleniem resztowym srY a odchyleniem sY w regresji lipcowych sum opadów w Krynicy względem lipcowych sum opadów w Rabce.
Zadanie. Obliczyć wartość odchylenia sY i wykreślić – oprócz prostych yśr(x)+srY i yśr(x)-srY,
jak w przykładzie 8.3 – proste y = y, y = y+ sY i y = y– sY.
Rozwiązanie. Wykorzystując rys. 8.6 i uzyskane już w przykładzie 8.2 wyniki: y = 106.61 mm i sY = 45.11 mm uzyskujemy wynik pokazany na rys. 8.7.
Rys. 8.6. Prosta yśr(x) regresji sumy opadów lipcowych w Krynicy względem sumy opadów lipcowych w Rabce wraz z odchyleniem resztowym srY = 29.7 mm
Interpretacja. Jeśli nic nie wiemy o związku regresyjnym zmiennej Y z inną zmienną, to przewidywanie (predykcja) wartości y tej zmiennej jest realizowane najczęściej w postaci y = y± sY. Jeśli na- tomiast znamy związek regresyjny y=yśr(x) i wartość x to przewidywanie wartości y ma postać y = yśr(x) ± srY. Wprowadzenie regresji poprawia więc predykcję na dwa sposoby: (1) wartość yśr(x) jest lepszą oceną ewentualnego y dla ustalonej wartości x niż wartość y oraz (2) niepewność sY związana z oceną yśr(x) jest mniejsza niż niepewność sY
związana z oceną y. Na przykład, jeśli interesuje nas odpowiedź na pytanie jaka będzie wartość y zmiennej Y dla x=200 mm, to z rys. 8.7 znajdujemy, że yśr(x) = 150 mm a niepewność tej predykcji wy-
nosi srY = 29.7 mm. 30 mm. Gdybyśmy nie znali postaci regresji, to odpowiedź mogłaby być tylko taka: y = 106.6 ± 45.1 mm.
Przykład powyższy pokazuje, w jaki sposób można wykorzystać odchylenie resztowe i jaki sens jest z nim związany. Nie jest to jednak ocena jakości otrzymanej prostej regresji, tzn. odchylenie resztowe (czy wariancja resztowa) nic nie mówi o niepewnoś- ci sYśr2 czy sYśr związanej z wartością yśr(x) otrzymaną dla ustalonej wartości x. Otrzy- manie i wykorzystanie takiej oceny jest przedmiotem następnych podrozdziałów.
8.2.4 Ocena niepewności oszacowania współczynników regresji liniowej Prosta regresji yśr(x) = ax+b otrzymana metodą najmniejszych kwadratów (czy inną metodą) jest estymatorem (czyli oceną) nieznanej rzeczywistej prostej regresji E(Y|x)
= αx + β. Estymator jest funkcją próby – oznacza to, że dla ustalonej wartości x zmiennej losowej X możliwych jest wiele wartości yśr(x). Dla każdej próby losowej pobranej z tej samej populacji dostaniemy trochę inne współczynniki regresji a i b. A więc – ze względu na te współczynniki – yśr(x) jest zmienną losową, co możemy zapisać następująco (por. (8.16))
śr( )
Y x = A x⋅ +B (8.26)
W zapisie tym zaznaczone zostało, że współczynniki regresji są zmiennymi losowymi (dlatego duże litery A i B) i dlatego również otrzymana wartość średnia Yśr(x) jest też
Rys. 8.7. Zysk w stosowaniu regresji sumy opadów lipcowych w Krynicy względem sumy opadów lipcowych w Rabce polega na lepszej predykcji wartości y sumy opadu w Krynicy (np. y=147 mm dla x=185 mm, zamiast y=106.6 mm dla wszystkich x) oraz mniejszym średnim błędzie tej predykcji (30 mm zamiast 45 mm)
zmienną losową. Ponieważ Yśr(x) jest zależna od x, jest to wartość średnia warun- kowa.
Przykład 8.5.. Prosta regresji z próby jest wielkością losową. Z po- pulacji normalnej N(µX,σX,µY, σY,ρ), gdzie µX = µY = 0, σX = 2 ,σY = 1 i ρ=0.8 wylosowano dwie 30-elementowe próby losowe (x1i,y1i) i(x2i,y2i), i=1,2 ,...,30; na podstawie każdej z nich obliczono współczynniki regresji a1, b1 oraz a2, b2 i wykreślono proste regresji, odpowied- nio yśr1(x) = a1 x + b1 i yśr2(x) = a2 x + b2. Wyniki przedstawione zostały na rys.
8.8. Wszystkie proste są najbliżej siebie mniej więcej w środku ciężkości danych.
W miarę oddalania się od miejsca, gdzie znajdują się dane zwiększa się odstęp pomiędzy prostą znalezioną na podsta- wie próby a rzeczywistą prostą regresji.
Aby ocenić niepewność związaną z uzyskaną wartością yśr(x) zmiennej Yśr(x) znajdziemy wariancje varA i varB. W tym celu zapiszmy nieco inaczej estymatory (8.20)współczynników regresji α i β:
2 2
1
1 ( )( )
n
i i
X X i
Yx x Y
A x x Y Y
s ns
B Y Ax
=
= − ⋅ = − −
= −
∑
(8.27)Zapis ten uwidacznia zwłaszcza te wielkości, które traktujemy jak zmienne losowe, co ma znaczenie dla dalszych obliczeń.
Aby lepiej zrozumieć sens sformułowań (8.27), pokażemy najpierw, że estymatory A i B są nieobciążonymi estymatorami parametrów α i β, potem znajdziemy ich wariancje, a na końcu obliczymy wariancję varYśr(x) oceny średniej warunkowej (8.26).
Nieobciążoność estymatorów A i B współczynników regresji αααα i ββββ. Przypomnijmy, że nieobciążoność estymatora pewnej wielkości oznacza, iż jego wartość średnia z próby jest równa prawdziwej wartości danej wielkości. Znajdźmy tę średnią EA najpierw dla współczynnika kierunkowego prostej regresji A.
Rys. 8.8. Prosta regresji yśr(x) jest obliczana na podstawie próby losowej i dlatego jest wielkością losową – realizacją rzeczywistej prostej regresji y(x)=αx + β
2 2
1 1 1
1 1
E E n ( i )( i ) E n ( i ) i n ( i )
i i i
X X
A x x Y Y x x Y Y x x
ns = ns = =
= − − = − − −
∑
∑ ∑
( )
2 2 2
1 1 1
2 2 2
2 2
1 1
1 1
( )E ( )( ) ( )
n n n
i i i i i i
i i i
X X X
n n
i i
i i
X X
x x Y x x x x x x
ns ns ns
x x x x x
ns s
α β α
α α
α
= = =
= =
= − = − + = −
= − = − =
∑ ∑ ∑
∑ ∑
(8.28)
W powyższych obliczeniach dwukrotnie wykorzystano przekształcony nieco wzór definicyjny średniej arytmetycznej: Σ(xi- x )=0. Podobnie, choć znacznie krócej, obliczamy średnią wartość EB:
1 1
1 1
E E E E ( )
( )
n n
i i
i i
B Y x A Y x x x
n n
x x
α α β α
α β α β
= =
= − ⋅ = − = + −
= + − =
∑ ∑
(8.29) Nieobciążoność estymatorów A i B pociąga za sobą nieobciążoność estymatora (8.26) prostej regresji E(Y|x):
EYśr( ) Ex = A x⋅ +EB=α⋅ +x β (8.30) Wariancje varA i varB estymatorów A i B współczynników regresji α i β. Teraz obliczymy wariancję varA:
( )
2 2
2
2 2
1 1
2 2 2
2 2
2 4
2 2
2
2 2
1
1 1
var var ( ) ( ) var
1 1
( )
n n
i i i i
i i
X X
n Y
i X
X i X X X
A x x Y x x Y
ns ns
x x s
ns ns ns ns
σ ρ
σ σ
σ
= =
=
= − = −
−
= − = = =
∑ ∑
∑
(8.31)W praktyce nie znamy prawdziwych wartości wariancji σY2 zmiennej Y i współczynnika korelacji ρ i wielkości te zastępujemy ich wartościami z próby, odpowiednio: sY2 i r. Zamiast (8.31) dostajemy
( )
2 2
2
( 2) 2 Y 1
A
X
s r
s s n
−
= − (8.32)
Obliczenie wariancji współczynnika B jest nieco bardziej skomplikowane
( ) ( )
2varB=var Y −Ax =varY −2 cov( , )x A Y + x varA (8.33) Trzeci składnik sumy jest już znany; pozostaje obliczyć pierwszy i drugi. Pierwszy składnik liczy się łatwo:
2 2
2 2
1 1
1 1 1
var var n i n var i
i i
Y Y Y n
n n n n
σ σ
= =
= = = =
∑
∑
(8.34)Aby obliczyć wartość drugiego zapiszemy A dla wygody w skróconej postaci:
2
1 1
( )
n n
i
i i i
i X i
x x
A Y c Y
= ns =
=
∑
− =∑
(8.35)gdzie Eci=0 (dlaczego?). Na mocy własności kowariancji mamy dalej
1 1 1 1
2
1 1
1 1
cov( , ) cov , cov( , )
1 var 0
n n n n
i i j i i j
i j i j
n n
i i i
i i
A Y c Y Y c Y Y
n n
c Y c
n n
σ
= = = =
= =
= =
= = =
∑ ∑ ∑ ∑
∑ ∑
(8.36)Suma po j została zniknęła wskutek drugiego założenia (8.14) – założenia o nieza- leżności zmiennych Yi. Uzyskane wyniki i wraz z pozwalają na kontynuowanie obli- czeń i uzyskanie końcowego wyniku na varB:
( ) ( ( ) )
2 2 2 2
2 2 2 2 2
2 2 2
var X var( )
X X X
B x s x x A x
n ns ns ns
σ σ σ σ
= + = + = = (8.37)
Oczywiście, podobnie jak w dla varA (8.31), nie znamy prawdziwych wartości wa- riancji σY2 zmiennej Y i współczynnika korelacji ρ i wielkość varA zastępujemy jej wartością z próby sA2 . Dostajemy więc
2 2 2
B A
s =s x (8.38)
Jesteśmy teraz przygotowani do obliczenia niepewności varYśr(x) związanej z oceną Yśr(x) nieznanej prostej regresji E(Y|x) = αx+β. Wielkości: varYśr(x) lub, częściej odchylenie standardowe DYśr(x), mogą służyć jako miary niepewności oceny prostej regresji Yśr(x).
Wariancja varYśr(x) średniej warunkowej z próby Yśr(x). Biorąc wariancję obu stron równania (8.26) dostajemy
varYśr( ) var(x = A x⋅ +B)=x2varA+2 cov( , ) varx A B + B=... (8.39) Aby dopisać dalszy ciąg równania niezbędne jest obliczenie tylko jednej wielkości:
cov(A,B), gdyż pozostałe składniki są już obliczone. Korzystając z drugiego ze wzorów (8.27) dostajemy:
cov( , ) cov( ,A B = A Y −Ax) cov( , )= A Y −xcov( , )A A = −xvarA (8.40) Po drodze wykorzystaliśmy udowodnioną wcześniej własność (8.36). Możemy więc pisać dalej:
( ) ( )
( ) ( )
2 2 2 2 2
2 2 2 2 2 2 2
... var 2 var var 2
( ) ( )
A
A A X
x A xx A x A x xx x
x x x x x x s
σ
σ σ
= + + = + +
= − + − = − + (8.41)
Jeśli wykorzystamy przedostatnią równość we wzorze (8.31), otrzymamy końcowy wzór na niepewność prostej regresji z próby w punkcie x
2 2 2
2 2
2 2
( )
var śr ( ) X 1
X X
x x
Y x x s
ns n s
σ σ −
= − + = +
(8.42)
Wyrażenie to jest wyrażeniem teoretycznym; w praktyce musimy stosować estymator
2
|
sY x wariancji resztowej σY x2| , co prowadzi do nieco innej wersji wzoru (8.42):
2 2
2 |
2
( )
1
Y x Yśr
X
s x x
s n s
−
= +
(8.43)
Uzyskane wyrażenia (8.42) i (8.43) pokazują, że:
1. niepewność uzyskanej oceny Yśr(x) prostej regresji y(x)=αx+β nie jest stała: jest najmniejsza w punkcie x = x i zwiększa się symetrycznie (i nieograniczenie) względem tego punktu;
2. niepewność oceny uzyskanej prostej regresji maleje do zera wraz z n zmierzającym do nieskończoności, co oznacza, że im bardziej liczebna próba tym bardziej ocena prostej regresji jest zbliżona do rzeczywistej prostej regresji.
Przykład 8.6 ilustruje stwierdzenie 1.
Niepewność pojedynczej przewidywanej wartości zmiennej Y. Jeśli chce- my przewidzieć, jaka wartość y zostanie zrealizowana w doświadczeniu losowym (x,Y) posłużymy się liniową wersją wzoru (8.8) (zob. też (8.13)):
ˆi śr( )i i i i
Y =Y x + ∈ = Ax +B+ ∈ (8.44)
Uwzględniając (8.42) i (8.14)(ii) możemy obliczyć błąd predykcji pojedynczej wartości zmiennej Y, liczony jako wariancja
( )
2 2 2 2
2
2 2
var ˆ var ( ) var var
( ) ( )
1 1
i śr i i śr
X X
Y Y x Y
x x x x
n s n n s
σ σ
σ
= + ∈ = + ∈
− −
= + + = + +
(8.45)
Druga równość w pierwszym wierszu (8.45) jest prawdziwa z powodu braku korelacji pomiędzy Yśr a ∈.
8.2.5 Ocena jakości estymacji: (1-αααα)100% obszar ufności dla prostej regresji
Zakładamy, że zmienne losowe X i Y są zależne i podlegają dwuwymiarowemu rozkła- dowi normalnemu. Wyznaczenie prostej regresji yśr(x) (8.16) (tj. wyznaczenie jej parametrów (8.20)) nie daje informacji o dokładności, z jaką ta prosta odwzorowuje rzeczywistą prostą regresji y(x) = αx+β. Pewną informację na ten temat daje – jak to jest pokazane w poprzednim podrozdziale – wariancja sYśr2 (8.43) charakteryzująca średni błąd oceny średniej warunkowej E(Y|x) przy znanej wartości x zmiennej X.
Ocenę jakości dopasowania znalezionej prostej regresji możemy też przeprowadzić w sposób charakterystyczny dla estymacji przedziałowej: ustalamy dla każdej wartości x warunkowy przedział ufności dla zmiennej Y. W ten sposób, przy założonym poziomie ufności (1-α), określamy pewien obszar płaszczyzny (x,y), zwany obszarem ufności (dokładniej: (1-α)@100% obszarem ufności), jako obszar zawarty pomiędzy dwoma krzywymi: górną y+(x) i dolną y-(x), zawierającymi na danym poziomie ufności rzeczywistą prostą regresji E(Y|x) = αx+β:
2
2
( )
( ) śr( ) Yśr( ) p( , 2) rY 1
X
s x x
y x y x s x ax b t n
n s
± α
= ± = + ± − + − (8.46)
gdzie tp(α,n-2) jest odpowiednim kwantylem wziętym z rozkładu Studenta o n-2 stopniach swobody. Funkcje te są wykreślone na rys. 8.9.
Przykład 8.6. Niepewność sY x2| uzyskanej oceny yśr(x) prostej regresji y(x)=αx+β lipcowych sum opadów w Krynicy względem
lipcowych sum opadów w Rabce.
Zadanie. Na podstawie wzoru (8.43) i uzyskanych już wcześniej (na podsta- wie danych z tabeli 8.1) wyników: sY*x = 29.7 mm, x= 125.3 mm i sX = 57.7 mm znaleźć odchylenie standardowe sYśr(x) i wykreślić funkcje y+(x) = yśr(x) + sYśr(x) i y-(x) = yśr(x) - sYśr(x).
Rozwiązanie zadania pokazane jest na rys. 8.9. Zgodnie z niepewność oce- ny yśr(x) rzeczywistej prostej regresji jest najmniejsza w punkcie = 125.3 mm i zwiększa się w miarę oddalania się wzdłuż osi x od tego punktu.
Interpretacja. Korzystając z rys. 8.9 możemy odpowiedzieć na następujące
y+(x)
y-(x)
yśr(x)
sY****x=29.7 mm
y
x
Rys. 8.9. Krzywe y+(x) i y-(x) ograniczające pas o szero- kości pionowej 2sYśr(x)ilustrujący wielkość niepewności związanej z oceną yśr(x) rzeczywistej prostej regresji
pytanie: Jaka jest wartość prawdziwej prostej regresji E(Y*x)=αx+β w punkcie x? Odpowiedź brzmi: E(Y*x) . yśr(x) ± sYśr(x) = ax+b ± sYśr(x).
(1-α)100% obszar ufności jest miarą dokładności prostej regresji w tym sensie, w jakim (1-α)100% przedział ufności jest miarą dokładności parametru g. Posługując się tą analogią możemy
a) przy spełnionych założeniach dotyczących tworzenia obszaru ufności, b) przy ustalonym poziomie ufności 1-α
c) i przed wyznaczeniem tego obszaru
powiedzieć, że obszar ten będzie z prawdopodobieństwem 1-α zawierał rzeczywistą prostą regresji. Oczywiście, po wyznaczeniu obszaru ufności dla prostej regresji, zawiera on albo nie zawiera rzeczywistej prostej regresji i podkreślone stwierdzenie jest nieprawdziwe (nie spełniony warunek c).
Uwagi dotyczące regresji nieliniowej. Często spotykanym przypadkiem w analizie regresyjnej jest przypadek, gdy postulowana zależność y(x) nie może mieć postaci linii prostej. Mamy wtedy do czynienia z regresją nieliniową. Istnieją przy tym dwa najczęściej spotykane podejścia:
(i) Metoda bezpośrednia: obliczyć współczynniki krzywej y(x) korzystając bezpośrednio z metody najmniejszych kwadratów , gdzie zamiast ax+b wpiszemy inną postać postulowanej zależności y(x), np. y=a@xb.
(ii) Linearyzacja zagadnienia. Znaczy to, że należy przetransformować zmienne (x,y) na takie zmienne (x',y'), które będą zależały liniowo od siebie. Jeśli np.
przypuszczamy, że zależność y(x) ma postać y=a@xb, to logarytmując to równanie dostajemy lny=lna+b@lnx. Dokonując formalnego podstawienia y'=lny, a'=lna i x'=lnx mamy: y'=a'+bx', a więc zagadnienie liniowe.
Linearyzacja wprowadza pewnego rodzaju obciążenie, np. poprzez transformację logarytmiczną mniejsze wartości dostają większą wagę, a większe – mniejszą.
Jednakże takie podejście jest często spotykane w praktyce. Argumentem decydującym, jak to często bywa, jest prostota dalszych obliczeń.
PRZYKŁAD 8.7. Mamy daną 50-elementową dwuwymiarową prostą próbę losową przedstawioną na rys. 8.10A w postaci chmury punktów. Widać wyraźną nieliniowość zależności y(x). Załóżmy, że zależność tę będziemy przybliżać równaniem
y= ⋅a xb (8.47)
Rys. 8.10. Przykładowe porównanie metody bezpośredniej i linearyzacji w regresji nieliniowej; r2 jest wartością współczynnika determinacji będącego miarą dopasowania
Analogon równania (8.17) dla tego przypadku będzie następujący:
2
1
( , ) ) min
n
b
i i
i
F a b y ax
=
=
∑
− = (8.48)Pochodne cząstkowe funkcji F po a i b mają postać
( )
( )
1
1
( , )
2 ) 0
( , )
2 ) ln 0
n
b b
i i i
i n
b b
i i i i
i
F a b
y ax x
a F a b
y ax ax x
b
=
=
∂ = − − =
∂
∂ = − − =
∂
∑
∑
(8.49)Używając symbolu "kreska ponad" jako znaku średniej arytmetycznej, możemy powyższe równanie zapisać następująco
2
2
0
ln ln 0
b b
b b
yx a x
yx x a x x
− =
− = (8.50)
Po wyrugowaniu a dostaniemy dosyć skomplikowane równanie na parametr b:
2 2
ln 0 ln
b b
b b
yx yx x
x x x
− = (8.51)
Można je rozwiązać tylko numerycznie, co wykonuje się najczęściej za pomocą komputera.
Zostało ono rozwiązane, a wynik w postaci krzywej przedstawiony jest na rys. 8.10B. Ocena jakości dopasowania nie może być w tym przypadku wyrażona przez współczynnik korelacji liniowej, gdyż zagadnienie jest nieliniowe; do tego celu użyto współczynnika determinacji R2, zdefiniowanego równaniem
( )
( )
( )
( )
2 2
2
2 1 1
2 2 2
1 1
1
1 1 1
1
n n
b b
i i i i
i i
rY
n n
Y
i i
i i
y ax y ax
s n
R s
y y y y
n
= =
= =
− −
= − = − = −
− −
∑ ∑
∑ ∑
(8.52)Jego wartość jest podana na rys. 8.10B.
Drugi sposób znajdowania krzywej regresji – linearyzację zagadnienia – ilustrują kolejno rys. 8.10C i D (linie przerywane). Jeśli bowiem zlogarytmujemy stronami równanie (8.47),
lny=lna+blnx (8.53)
to traktując lnx i lny jako nowe zmienne, odpowiednio x′ i y′ i korzystając ze wzorów (8.20) dostaniemy równanie liniowe. Rys. 8.10C jest właśnie ilustracją tego podejścia.
Podana wartość r2 jest w tym przypadku zarówno wartością współczynnika korelacji liniowej, jak i wartością współczynnika determinacji. Jeśli teraz przejdziemy do zwykłych współrzędnych (x,y), to prosta regresji stanie się krzywą regresji, co pokazane jest linią przerywaną na rys. 8.10D. Dla porównania wykreślona jest tam również krzywa z rys. 8.10B. Jak widać różnice są niewielkie.
