Zagadnienia

8.1 Predykcja cechy ciągłej
8.2 Metoda najmniejszych kwadratów (mnk)
8.3 Inne wyprowadzenie estymatora najmniejszych kwadratów
8.4 Estymatory metody największej wiarygodności prametrów modelu liniowego
8.5 Kolejne własności estymatorów mnk
8.6 Model liniowy przy założeniu normalności
8.7 Test ilorazu wiarygodności (Likelihood Ratio Test) hipotez liniowych
8.8 Popularne kryteria wyboru modelu – kryteria informacyjne
8.9 Model logistyczny – przykład uogólnionego modelu liniowego
8.10 Przykłady w programie R

8. Modele liniowe

8.1. Predykcja cechy ciągłej

Będziemy obserwować ciągłą zmienną objaśnianą $y_{i}$ oraz zmienne objaśniające $x_{i}=(x_{{i1}},\ldots,x_{{ip}})$ , $i=1,\ldots,n$ . Na ich podstawie będziemy chcieli znaleźć funkcję zależącą od $x$ , która będzie najlepiej przybliżać cechę $y$ . Ograniczymy się przy tym tylko do zależności liniowej. Na podstawie znalezionej funkcji, dla nowo zaobserwoawanych $x_{{n+1}}$ będziemy mogli znależć predykcję $\widehat{y}_{{n+1}}$ . Jeżeli rozpatrzymy jednowymiarowy $x_{i}$ ( $p=1$ ), szukanie funkcji liniowej najlepiej przybliżającej dane obrazuje rysunek 8.1.

$\par$

Rys. 8.1. Regresja liniowa jako dopasowanie prostej do danych.

Dane są postaci:

$\underbrace{\left(\begin{array}[]{c}y_{1}\\ \ldots\\ y_{n}\\ \end{array}\right)}_{{\text{losowe,obserwowane}}}=\underbrace{\left(\begin{array}[]{ccc}x_{{11}}&\ldots&x_{{1p}}\\ \ldots&\ldots&\ldots\\ x_{{n1}}&\ldots&x_{{np}}\\ \end{array}\right)}_{{\text{deterministyczne,obserwowane}}}\cdot\underbrace{\left(\begin{array}[]{c}\beta _{1}\\ \ldots\\ \beta _{p}\\ \end{array}\right)}_{{\text{szukane}}}+\underbrace{\left(\begin{array}[]{c}\varepsilon _{1}\\ \ldots\\ \varepsilon _{n}\\ \end{array}\right)}_{{\text{losowe,nieobserwowane}}}.$

W zapisie macierzowym:

$y=X\beta+\varepsilon,$

gdzie:

wektor $y$ będziemy nazywać zmienną objaśnianą;
macierz $X$ macierzą planu;
$\beta$ estymowanymi parametrami;
$\varepsilon$ to wektor efektów losowych (wektor realizacji zmiennej losowej).

Dla tak sformułowanego problemu przyjmiemy następujące założenia:

$\mathbb{E}(\varepsilon)=0$ ;
$\text{Var}(\varepsilon)=\sigma^{2}I_{n}$ ;
rząd macierzy $X$ jest pełny: $\text{rank}(X)=p$ .

Dla tak sformułowanych danych, problem szukania estymatora parametru $\beta$ będziemy nazywać problemem liniowym.

Twierdzenie 8.1

Rozkład QR macierzy

Szeroki rozkład QR: Każdą rzeczywistą macierz $A$ wymiaru $n\times m$ $(m\leq n)$ można zapisać jako iloczyn maierzy ortogonalnej $Q$ wymiaru $n\times n$ $(Q^{T}Q=I_{n})$ oraz górnotrójkątnej macierzy $R$ wymiaru $n\times m$ :

$A=QR=\left(\begin{array}[]{ccc}&&\\ &Q&\\ &&\\ \end{array}\right)\left(\begin{array}[]{cc}\circ&\circ\\ &\circ\\ &\\ \end{array}\right).$

Wąski rozkład QR: Ponieważ $(n-m)$ dolnych wierzy macierzy $R$ jest zerowa, można skrócić zapis do

$A=QR=Q\left(\begin{array}[]{c}R\\ 0\\ \end{array}\right)=\left(\begin{array}[]{c}Q_{1}\quad Q_{2}\\ \end{array}\right)\cdot\left(\begin{array}[]{c}R_{1}\\ 0\\ \end{array}\right)=Q_{1}R_{1}=,$

$=\left(\begin{array}[]{cc}&\\ &Q_{1}\\ &\\ \end{array}\right)\left(\begin{array}[]{cc}\circ&\circ\\ &\circ\\ \end{array}\right).$

gdzie $Q_{1}$ jest macierzą wymairu $n\times m$ o ortogonalnych kolumnach a $R_{1}$ jest macierzą górnotrójkątną wymairu $m\times m$ .

Wąski rozkład QR jest zapisem macierzowym ortogonalizacji Gramma-Schmidta układu wektorów będących kolumnami macierzy $A$ . Szeroki rozkład otrzymujemy dopełniając macierz $Q_{1}$ do bazy przestrzeni $\mathbb{R}^{n}$ .

∎

Problem 8.1

Przy założeniu modelu liniowego, będziemy chcieli wyestymować nieznane parametry: $\beta$ i $\sigma^{2}$ .

8.2. Metoda najmniejszych kwadratów (mnk)

Zauważmy, że

$\mathbb{E}(y)=\mathbb{E}(X\beta+\varepsilon)=X\beta+\mathbb{E}(\varepsilon)=X\beta.$

Estymator najmnieszych kwadratów parametru $\beta$ to taka jego wartość, dla której odległości euklidesowe przybliżanych danych od prostej je przybliżających jest najmniejsza:

$\widehat{\beta}=\min _{{\beta}}||y-X\beta||^{2}.$

Twierdzenie 8.2

Estymator najmniejszych kwadratów wyraża się wzorem

$\widehat{\beta}=R_{1}^{{-1}}Q_{1}^{T}y,$

gdzie $R_{1}$ i $Q_{1}$ pochodzą z wąskiego rozkładu QR macierzy planu $X$ .

Skorzystajmy z szerokiego rozkładu QR macierzy $X$ : $X=QR$ . Ponieważ mnożenie wektora przez macierz ortogonalną nie zmienia jego normy, możemy zapisać:

$||y-X\beta||^{2}=||Q^{T}(y-QR\beta)||^{2}=$

$=||Q_{1}^{T}(y-Q_{1}R_{1}\beta)||^{2}+||Q_{2}^{T}(y-Q_{1}R_{1}\beta)||^{2}=$

$=||Q_{1}^{T}y-R_{1}\beta||^{2}+||Q_{2}^{T}y||^{2}.$

Wyrażenie to osiąga minimum ze względu na parametr $\beta$ , jeżeli wyzerujemy pierwszy składnik sumy:

$Q_{1}^{T}y=R_{1}\widehat{\beta}.$

Ponieważ macierz $R_{1}$ jest kwadratowa i pełnego rzędu (rank $(X)=p$ ), możemy ją odwrócić:

$\widehat{\beta}=R_{1}^{{-1}}Q_{1}^{T}y.$

∎

Wniosek 8.1

Zauważmy, że:

Predykcja dla $y$ jest równa $\widehat{y}=X\widehat{\beta}$ ;
$||y-\widehat{y}||^{2}=||Q_{2}^{T}y||^{2}.$ (8.1)

Przyjrzyjmy się własnościom metody najmniejszych kwadratów (zostaną one udowodnione w dalszej części wykładu):

$\widehat{y}=X\widehat{\beta}$ jest rzutem ortogonalnym $y$ na przestrzeń rozpiętą przez kolumny macierzy planu $X$ .
Nieobciążonym estymatorem parametru $\sigma^{2}$ jest $\widehat{\sigma}^{2}=\frac{||y-\widehat{y}||^{2}}{n-p}$ .
Twierdzenie Gaussa-Markowa: estymator $\widehat{\beta}$ jest liniowym, nieobciążonym estymatorem o najmniejszej wariancji parametru $\beta$ (BLUE- Best Linear Unbiased Estimator).
Przy założeniu $\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$ , zachodzi twierdzenie Fishera:
- $\widehat{\beta}\sim\mathcal{N}(\beta,\sigma^{2}(X^{T}X)^{{-1}})$ ;
- $\frac{(n-p)\widehat{\sigma}^{2}}{\sigma^{2}}\sim\chi^{2}(n-p)$ ;
- $\widehat{\beta}$ i $\widehat{\sigma}^{2}$ są niezależne.

8.3. Inne wyprowadzenie estymatora najmniejszych kwadratów

Wyprowadzimy estymator mnk jako rozwiązanie zadania BLUE – liniowy, nieobciążony estymator o najmniejszej wariancji. Rozumowanie będzie jednocześnie dowodem twierdzenia Gaussa-Markowa.

Twierdzenie 8.3

Dla problemu liniowego estymator postaci $\widehat{\beta}=(X^{T}X)^{{-1}}X^{T}y$ jest liniowym, nieobciążonym estymatorem o najmniejszej wariancji parametru $\beta$ .

Żeby łatwiej mówić o nieobciążoności, czy minimalnej wariancji, zredukujemy wymiar problemu do jednowymiarowego zakładając, że własności będą zachodzić dla wszystkich możliwych kombinacji liniowych zmiennej objaśnianej: Dla danego wektora $l$ , konstruujemy kombinację liniową $\mu=l^{T}\beta$ i szukamy dla niej estymatora zależącego liniowo od $y$ :

$\widehat{\mu}=c^{T}y$

przy założeniu nieobciążoności:

$\mathbb{E}(\widehat{\mu})=\mu=l^{T}\beta.$

Jednocześnie wiemy, że:

$\mathbb{E}(\widehat{\mu})=c^{T}\mathbb{E}(y)=c^{T}X\beta.$

Stąd:

$c^{T}X=l^{T}.$

Będziemy minimalizować wariancję estymatora $\widehat{\mu}$ :

$\text{Var}(\widehat{\mu})=\text{Var}(c^{T}y)=c^{T}\text{Var}(y)c=\sigma^{2}c^{T}c.$

Zadanie optymalizacyjne wygląda następująco:

$\min _{{c}}\left[\sigma^{2}c^{T}c\right]\text{ pod warunkiem }c^{T}X=l^{T}.$

Skorzystajmy z metody mnożników Lagrange'a:

$\min _{c}\left[F(c,\lambda)=\sigma^{2}c^{T}c+(c^{T}X-l^{T})\underbrace{\lambda}_{{\text{wektor}}}\right].$

Szukamy estymatora wektora $c$ , spełniającego dwa równania:

$\frac{\partial F}{\partial c}=2\sigma^{2}c+X\lambda=0;$

(8.2)

$\frac{\partial F}{\partial\lambda}=X^{T}c-l=0.$

(8.3)

Z równania 8.2 otrzymujemy: $c=-\frac{1}{2\sigma^{2}}X\lambda$ , wstawiamy do równania 8.3:

$-\frac{1}{2\sigma^{2}}X^{T}X=l,$

skąd:

$\lambda=-2\sigma^{2}(X^{T}X)^{{-1}}l.$

Macierz $X$ jest pełnego rzędu, więc macierz $(X^{T}X)$ jest odwracalna. Wstawiając $\lambda$ do wzoru na $c$ , otrzymujemy:

$c=X(X^{T}X)^{{-1}}l\text{ , inaczej: }c^{T}=l^{T}(X^{T}X)^{{-1}}X^{T}.$

Estymator $\widehat{\mu}$ jest więc postaci:

$\widehat{\mu}=c^{T}y=l^{T}(X^{T}X)^{{-1}}X^{T}y;$

podstawiając za $l$ kolejne wektory bazy kanonicznej $e_{i}=(0,\ldots,0,1,0,\ldots,0)$ , znajdujemy kolejne estymatory kombinacji liniowych $l^{T}\beta=e_{i}^{T}\beta=\widehat{\beta _{i}}$ , co łącznie możemy zapisać jako:

$\widehat{\beta}=(X^{T}X)^{{-1}}X^{T}y.$

∎

Stwierdzenie 8.1

Liniowy, nieobciążony estymator o najmniejszej wariancji parametru $\beta$ w modelu liniowym jest równy estymatorowi najmniejszych kwadratów.

$\hat{\beta}=(X^{T}X)^{{-1}}X^{T}y=$

korzystając z wąskiego rozkładu QR: $X=Q_{1}R_{1}$ ,

$=(R_{1}^{T}\underbrace{Q_{1}^{T}Q_{1}}_{{=I_{p}}}R_{1})^{{-1}}R_{1}^{T}Q_{1}^{T}y=$

$=R_{1}^{{-1}}\underbrace{(R_{1}^{T})^{{-1}}R_{1}^{T}}_{{=I_{p}}}Q_{1}^{T}y=R_{1}^{{-1}}Q_{1}^{T}y.$

∎

8.4. Estymatory metody największej wiarygodności prametrów modelu liniowego

Estymatory największej wiarygodności to takie wartości parametrów, których prawdopodobieństwo zaobserwowania danych jest największe. Żeby skorzystać z tej metody estymacji, potrzebna jest funkcja wiarygodności, niezbędne więc będzie założenie na temat rozkładu danych:

Zamiast zakładać:

$\mathbb{E}\varepsilon=0\quad,\quad\text{Var}(\varepsilon)=\sigma^{2}I_{n}$

założymy:

$\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n}),$

skąd mamy $y\sim\mathcal{N}(X\beta,\sigma^{2}I_{n})$ oraz funkcję wiarygodności:

$f_{{\beta,\sigma^{2}}}(y_{1},\ldots,y_{n})=f_{{\beta,\sigma^{2}}}(y)=\frac{1}{(2\pi)^{{\frac{n}{2}}}|\sigma^{2}I_{n}|^{{\frac{1}{2}}}}\exp\left[-\frac{1}{2}(y-X\beta)\sigma^{{-2}}I_{n}(y-X\beta)\right]=$

$=\frac{1}{(2\pi)^{{\frac{n}{2}}}(\sigma^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{1}{2\sigma^{{2}}}||y-X\beta||^{2}\right].$

Funkcję wiarygodności będziemy chcieli zmaksymalizować po parametrach $\beta$ i $\sigma^{2}$ . Ponieważ logarytm jest funkcją rosnącą, jest to równoważne z maksymalizacją logarytmu funkcji wiarygodności:

$F_{{\beta,\sigma^{2}}}(y)=\log f_{{\beta,\sigma^{2}}}(y)=C-\frac{n}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{{2}}}||y-X\beta||^{2},$

gdzie $C$ jest stałą niezależną od szukanych parametrów. Zadanie maksymalizacji logwiarygodności $F_{{\beta,\sigma^{2}}}$ jest rónoważne minimalizacji $-2F_{{\beta,\sigma^{2}}}$ :

$-2F_{{\beta,\sigma^{2}}}(y)=C^{{\prime}}+n\log(\sigma^{2})+\frac{1}{\sigma^{{2}}}||y-X\beta||^{2}.$

Część sumy zależąca od parametru $\beta$ to $\frac{1}{\sigma^{{2}}}||y-X\beta||^{2}$ . Wartością parametru $\beta$ minimalizującą to wyrażenie jest:

$\widehat{\beta}_{{mnw}}=(X^{T}X)^{{-1}}X^{T}y=R_{1}^{{-1}}Q_{1}^{T}y,$

co udowodniliśmy już w twierdzeniu 8.2.

Ponieważ $\widehat{\beta}$ nie zależy od parametru $\sigma^{2}$ , mogę wstawić estymator do funkcji wiarygodności przy szukaniu optymalnego parametru $\sigma^{2}$ . Oznaczmy także $\tau=\sigma^{2}$ żeby nie mylił się nam kwadrat przy parametrze:

$-2F_{{\beta,\sigma^{2}}}(y)=C^{{\prime}}+n\log(\tau)+\frac{1}{\tau}||y-X\widehat{\beta}||^{2}.$

$\frac{\partial(-2F_{{\beta,\sigma^{2}}})}{\partial\tau}=\frac{n}{\tau}-\frac{1}{\tau^{2}}||y-X\widehat{\beta}||^{2}=0,$

skąd otrymujemy:

$n\tau=||y-X\widehat{\beta}||^{2};$

$\widehat{\tau}=\widehat{\sigma}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n}.$

Wniosek 8.2

Przy założeniu rozkładu normalnego:

estymatory parametru $\beta$ dla metody największej wiarygodności i metody najmniejszych kwadratów są równe:

$\widehat{\beta}_{{mnk}}=\widehat{\beta}_{{mnw}}=(X^{T}X)^{{-1}}X^{T}y;$
estymatory parametru $\sigma^{2}$ dla metody największej wiarygodności i metody najmniejszych kwadratów są równe z dokładnością do stałej:

$\widehat{\sigma}_{{mnk}}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n-p},\quad\widehat{\sigma}_{{mnw}}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n}.$

8.5. Kolejne własności estymatorów mnk

8.5.1. Wartość oczekiwana i wariancja estymatora $\widehat{\beta}$

Wartość oczekiwana:

$\mathbb{E}\widehat{\beta}=\mathbb{E}(X^{T}X)^{{-1}}X^{T}y=(X^{T}X)^{{-1}}X^{T}\mathbb{E}y=(X^{T}X)^{{-1}}X^{T}X\beta=\beta.$ (8.4)

Estymator jest niebciążony.
Macierz wariancji:

$\text{Var}(\widehat{\beta})=\mathbb{E}(\widehat{\beta}-\beta)(\widehat{\beta}-\beta)^{T}=$

$=\mathbb{E}((X^{T}X)^{{-1}}X^{T}y-(X^{T}X)^{{-1}}X^{T}\mathbb{E}y)((X^{T}X)^{{-1}}X^{T}y-(X^{T}X)^{{-1}}X^{T}\mathbb{E}y)^{T}=$

$=(X^{T}X)^{{-1}}X^{T}\mathbb{E}\left[(y-\mathbb{E}y)(y-\mathbb{E}y)^{T}\right]((X^{T}X)^{{-1}}X^{T})^{T}=$

$=(X^{T}X)^{{-1}}X^{T}\sigma^{2}I_{n}((X^{T}X)^{{-1}}X^{T})^{T}=\sigma^{2}(X^{T}X)^{{-1}}\underbrace{(X^{T}X)(X^{T}X)^{{-1}}}_{{=I_{p}}}=$

$=\sigma^{2}(X^{T}X)^{{-1}}.$ (8.5)

8.5.2. Dopasowanie $\widehat{y}$ jako rzut ortogoanlny $y$ na przestrzeń rozpiętą przez kolumny macierzy $X$

$\par$

Rys. 8.2. Dopasowanie $\widehat{y}$ jako rzut ortogonalny $y$ na $\text{lin}(X)$ .

Prypomnijmy

$\widehat{\beta}=(X^{T}X)^{{-1}}X^{T}y\quad,\quad\widehat{y}=X\widehat{\beta}.$

Definicja 8.1

Macierzą daszkową $H$ nazwiemy taką macierz, że:

$\widehat{y}=X(X^{T}X)^{{-1}}X^{T}y=Hy.$

Stąd:

$H=X(X^{T}X)^{{-1}}X^{T}.$

Stwierdzenie 8.2

Zauważmy, że $\widehat{y}$ jest nieobciążonym estymatorem $y$ :

$\mathbb{E}\widehat{y}=X(X^{T}X)^{{-1}}X^{T}X\beta=X\beta.$

Stwierdzenie 8.3

Własności macierzy daszkowej $H$ :

$HX=X$ :

$X\underbrace{(X^{T}X)^{{-1}}X^{T}X}_{{=I_{p}}}=X;$
macierz $H$ jest idempotentna, czyli $HH=H$ :

$HH=X\underbrace{(X^{T}X)^{{-1}}X^{T}X}_{{=I_{p}}}(X^{T}X)^{{-1}}X^{T}=H;$
symetryczna, czyli $H^{T}=H$ :

$H^{T}=(X(X^{T}X)^{{-1}}X^{T})^{T}=X(X^{T}X)^{{-1}}X^{T}=H.$
korzystając z wąskiego rozkładu QR macierzy $X$ , $H=Q_{1}Q_{1}^{T}$ :

$H=Q_{1}R_{1}(R_{1}^{T}\underbrace{Q_{1}^{T}Q_{1}}_{{=I_{p}}}R_{1})^{{-1}}R_{1}^{T}\underbrace{Q_{1}^{T}Q_{1}}_{{=I_{p}}}\underbrace{R_{1}R_{1}^{{-1}}}_{{=I_{p}}}\underbrace{(R_{1}^{T})^{{-1}}R_{1}^{T}}_{{=I_{p}}}Q_{1}^{T}=$

$=Q_{1}\underbrace{R_{1}R_{1}^{{-1}}}_{{=I_{p}}}\underbrace{(R_{1}^{T})^{{-1}}R_{1}^{T}}_{{=I_{p}}}Q_{1}^{T}=Q_{1}Q_{1}^{T}.$
korzystając z szerokiego rozkładu QR macierzy $X$ , możemy przyjrzeć się rozkładowi spektralnemu macierzy daszkowej:

$H=QR(R^{T}\underbrace{Q^{T}Q}_{{=I_{n}}}R)^{{-1}}R^{T}Q^{T}=$

ponieważ $R^{T}R=R_{1}^{T}R_{1}$ ,

$=QR(R_{1}^{T}R_{1})^{{-1}}R^{T}Q^{T}=QRR_{1}^{{-1}}(R_{1}^{T})^{{-1}}R^{T}Q^{T}=$

$=\underbrace{Q\left(\begin{array}[]{cccc}1&&&\\ &\ldots&&\\ &&1&\\ &&&0\\ \end{array}\right)Q^{T}}_{{\text{rozkład spektralny macierzy }H}}=$

dla $Q=[q_{1},\ldots,q_{p},\ldots,q_{n}]$ :

$=Q_{1}Q_{1}^{T}=\sum _{{i=1}}^{p}1\cdot q_{i}q_{i}^{T}+\sum _{{i=p+1}}^{n}0\cdot q_{i}q_{i}^{T}$

Wniosek 8.3

Macierz daszkowa $H$ jest macierzą rzutu ortogonalnego na przestrzeń rozpiętą przez kolumny macierzy $X$ .
Jeżeli $\widehat{y}=X\widehat{\beta}$ minimalizuje wyrażenie $||y-\widehat{y}||^{2}$ , to jest rzutem ortogonalnym $y$ na $\text{lin}(X)$ .

8.5.3. Nieobciążony estymator parametru $\sigma^{2}$

Stwierdzenie 8.4

Macierz $(I_{n}-H)$ jest macierzą rzutu ortogonalnego na przestrzeń prostopadłą do przestrzeni rozpiętej przez kolumny macierzy $X$ , jest więc w szczególności symetryczna i idempotentna.

$I_{n}-H=I_{n}-Q\left(\begin{array}[]{cc}I_{p}&\\ &0\\ \end{array}\right)Q^{T}=$

$=QQ^{T}-Q\left(\begin{array}[]{cc}I_{p}&\\ &0\\ \end{array}\right)Q^{T}=Q\left(\begin{array}[]{cc}0&\\ &I_{{n-p}}\\ \end{array}\right)Q^{T}.$

∎

Wniosek 8.4

Ponieważ ślad macierzy równy jest sumie jego wartości własnych, ślady macierzy daszkowej $H$ i macierzy $(I-H)$ to:

$\text{tr}(H)=p;$

$\text{tr}(I-H)=n-p.$

Stwierdzenie 8.5

Twierdzenie Pitgorasa w postaci macierzowej:

$\text{Var}(y)=\text{Var}(\widehat{y})+\text{Var}(y-\widehat{y}).$

$\sigma^{2}I_{n}=\text{Var}(Hy)+\text{Var}((I-H)y);$

$\sigma^{2}I_{n}=H\text{Var}(y)H^{T}+(I-H)\text{Var}(y)(I-H)^{T};$

ponieważ macierze $H$ i $I-H$ są symetryczne i idempotentne, zachodzi:

$\sigma^{2}I_{n}=\sigma^{2}H+\sigma^{2}(I_{n}-H).$

(8.6)

∎

Stwierdzenie 8.6

Niebciążonym estymatorem parametru $\sigma^{2}$ w modelu liniowym jest:

$\widehat{\sigma}^{2}=\frac{||y-\widehat{y}||^{2}}{n-p}.$

Ponieważ $\widehat{y}$ jest nieobciążonym estymatorem $y$ , możemy zapisać:

$\mathbb{E}||y-\widehat{y}||^{2}=\sum _{{i=1}}^{n}\left(\text{Var}(y_{i}-\widehat{y}_{i})\right)=\text{tr}\left(\text{Var}(y-\widehat{y})\right)=$

$=\text{tr}\left(\sigma^{2}(I-H)\right)=\sigma^{2}(n-p).$

Stąd:

$\frac{\mathbb{E}||y-\widehat{y}||^{2}}{n-p}=\sigma^{2}.$

∎

8.5.4. Model z większą liczbą parametrów nie musi być lepiej dopasowany dla nowych danych

Błąd predykcji $y$ za pomocą $\widehat{y}$ na tej samej próbie, korzystając ze wzoru 8.6, można zapisać w postaci:

$\mathbb{E}(y_{i}-\widehat{y}_{i})^{2}=\text{Var}(y_{i}-\widehat{y}_{i})=\sigma^{2}(1-h_{{ii}}),$

gdzie $h_{{ii}}$ jest elementem macierzy daszkowej: $H=(h_{{ij}})_{{i,j=1}}^{n}$ .

Definicja 8.2

Elementy przekątnej macierzy daszkowej $H$ : $h_{{ii}}$ będziemy nazywać ładunkami obserwsacji $i$ -tej i oznaczać $h_{{i}}$ .
Dla nowych obserwacji mamy:

Zakładamy niezależność nowych obserwacji zmiennej objaśnianej i $p$ -wymiarowego wektora zmiennych objaśniających: $(y_{{n+1}},x_{{n+1}}^{T})$ od $(y,X)$ . Będziemy estymować parametry używając danych treningowych $(y,X)$ , a obliczać błąd dla nowych danych testowych:

$y_{{n+1}}=x_{{n+1}}^{T}\beta+\varepsilon _{{n+1}};$

$\widehat{\beta}=(X^{T}X)^{{-1}}X^{T}y;$

$\widehat{y}_{{n+1}}=x_{{n+1}}^{T}\widehat{\beta}.$

Błąd predykcji jest równy:

$\mathbb{E}(y_{{n+1}}-\widehat{y}_{{n+1}})^{2}=\text{Var}(y_{{n+1}}-\widehat{y}_{{n+1}})\underbrace{=}_{{\text{z niezależności}}}\text{Var}(y_{{n+1}})+\text{Var}(\widehat{y}_{{n+1}})=$

$=\sigma^{2}+\text{Var}(x_{{n+1}}^{T}\widehat{\beta})=\sigma^{2}+x_{{n+1}}^{T}\text{Var}(\widehat{\beta})x_{{n+1}}=$

$=\sigma^{2}+x_{{n+1}}^{T}\sigma^{2}(X^{T}X)^{{-1}}x_{{n+1}}=\sigma^{2}(1+h_{{n+1}}),$

gdzie $h_{{n+1}}=x_{{n+1}}^{T}(X^{T}X)^{{-1}}x_{{n+1}}$ , analogicznie do ładunków obserwacji dla $i=1,\ldots,n$ : $h_{i}=x_{i}^{T}(X^{T}X)^{{-1}}x_{i}$ .
Porównanie obu błędów predykcji dla tej samej macierzy planu:

Dane treningowe, dla których będziemy estymować parametr $\widehat{\beta}$ to $(y,X)$ gdzie $y=(y_{1},\ldots,y_{n})$ . Dane testowe, dla których będziemy liczyć błąd predykcji to w pierwszym przypadku ten sam zbiór $(y,X)$ , a w drugim $(y^{{te}},X)$ gdzie $y=(y_{{n+1}},\ldots,y_{{2n}})$ są nowymi obserwacjami, a macierz planu $X$ pozostaje niezmieniona. Porównajmy uśrednione oba błędy predykcji:

$\mathbb{E}\left(\frac{1}{n}\sum _{{i=1}}^{n}(y_{i}-\widehat{y}_{i})^{2}\right)=\frac{1}{n}\mathbb{E}||y-\widehat{y}||^{2}=\frac{(n-p)\sigma^{2}}{n};$

$\mathbb{E}\left(\frac{1}{n}\sum _{{i=n+1}}^{{2n}}(y_{i}-\widehat{y}_{i})^{2}\right)=\frac{1}{n}\sum _{{i=n+1}}^{{2n}}\sigma^{2}(1+h_{i})=\frac{(n+p)\sigma^{2}}{n},$

gdzie korzystamy z równości $h_{{n+i}}=h_{i}$ , co zachodzi dzięki użyciu tej samej macierzy planu $X$ w zbiorze treningowym i testowym oraz własnści macierzy daszkowej $\sum _{{i=1}}^{n}h_{{i}}=p$ .
Wniosek 8.5

Na podstawie obliczonych błędów predykcji możemy wywnioskować:
1. Większy model nie zawsze oznacza lepsze dopasowanie.
2. Różnica pomiędzy błędami predykcji wynosi:
  
  $\frac{1}{n}\mathbb{E}||y-\widehat{y}||^{2}-\frac{1}{n}\mathbb{E}||y^{{tr}}-\widehat{y}||^{2}=\frac{2p\sigma^{2}}{n}$

8.5.5. Kroswalidacja leave-one-out

Estymację błędu predykcji można oprzeć na kroswalidacji leave-one-out. Dla każdej obserwacji będziemy estymować model za pomocą wszystkich obserwacji oprócz niej samej i obliczać błąd predykcji na nowych danych dla tej pominiętej obserwacji. W ten sposób dostaniemy $n$ błędów predykcji, które następnie uśrednimy.

Niech $X_{{(i)}}$ , $i=1,\ldots,n$ oznacza macierz $X$ z usuniętą $i$ -tą obserwacją ( $i$ -tym wierszem), $y_{{(i)}}$ wektor obserwacji z usuniętą $i$ -tą obserwacją. Estymator $\widehat{\beta}_{{(i)}}$ będzie oznaczać estymator mnk na podstawie danych $(y_{{(i)}},X_{{(i)}})$ :

$\widehat{\beta}_{{(i)}}=(X_{{(i)}}^{T}X_{{(i)}})^{{-1}}X_{{(i)}}^{T}y_{{(i)}};$

Predykcja dla pominiętej obserwacji wyraża się wzorem:

$\widehat{y}_{{(i)}}=x_{i}^{T}\widehat{\beta}_{{(i)}},$

gdzie tak jak przy liczeniu błędu predykcji na nowych danych, $(y_{i},x_{i}^{T})$ jest niezależne od $(y_{{(i)}},X_{{(i)}})$ .

Korzystając z tego, że $\mathbb{E}(\widehat{y}_{{(i)}})=\mathbb{E}(y_{i})=x_{i}^{T}\beta$ , otrzymujemy:

$\mathbb{E}(\widehat{y}_{{(i)}}-y_{i})^{2}=\text{Var}(\widehat{y}_{{(i)}}-y_{i})=\sigma^{2}(1+x_{i}^{T}(X_{{(i)}}^{T}X_{{(i)}})^{{-1}}x_{i})=\frac{\sigma^{2}}{1-h_{{ii}}},$

gdzie $h_{{ii}}$ to $i$ -ty wyraz na przekątnej macierzy daszkowej dla pełnej macierzy $X$ : $H=X(X^{T}X)^{{-1}}X^{T}$ . Fakt ostatniej równości w powyższym wzorze przyjmiemy bez dowodu.

Wniosek 8.6

Estymator błędu predykcji przy użyciu kroswalidacji leave-one-out można uprościć do wzoru:

$\frac{||y-\widehat{y}_{{()}}||^{2}}{n}=\frac{1}{n}\sum _{{i=1}}^{n}(y_{i}-\widehat{y}_{{(i)}})^{2}=\sum _{{i=1}}^{n}\frac{\sigma^{2}}{1-h_{{ii}}}.$

(8.7)

8.6. Model liniowy przy założeniu normalności

Zamiast w modelu liniowym zakładać:

$\mathbb{E}\varepsilon=0\quad,\quad\text{Var}(\varepsilon)=\sigma^{2}I_{n},$

założymy:

$\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n}).$

Dzięki takiemu sformułowaniu zadania, będziemy mogli znaleźć rozkłady estymatorów $\widehat{\beta}$ i $\widehat{\sigma}^{2}$ , co umożliwi wnioskowanie statystyczne na ich temat, na przykład kondtrukcję przedziałów ufności. Udowodnimy:

Twierdzenie 8.4 (Fishera)

Przy założeniu $\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$ , estymatory modelu liniowego spełniają:

$\widehat{\beta}\sim\mathcal{N}(\beta,\sigma^{2}(X^{T}X)^{{-1}})$ ;
$\widehat{\beta}$ i $\widehat{\sigma}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n-p}$ są niezależne;
$\frac{(n-p)\widehat{\sigma}^{2}}{\sigma^{2}}\sim\chi^{2}(n-p)$ ;

Ponieważ $y=X\beta+\varepsilon$ , mamy:

$y\sim\mathcal{N}(X\beta,\sigma^{2}I_{n}).$

(8.8)

Wiemy, że nieobciążonymi estymatorami parametrów modelu liniowego są:

$\widehat{\beta}=(X^{T}X)^{{-1}}X^{T}y\quad,\quad\widehat{\sigma}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n-p}.$

Korzystając z twierdzenia 7.1 oraz wzoru 8.8, możemy wywnioskować, że $\widehat{\beta}$ ma rozkład normalny. Wartość oczekiwana estymatora (8.4) i wariancja (8.5) wyznaczają jednoznacznie rozkład.
Przypomnijmy wzory 8.2 i 8.1:

$\widehat{\beta}=R_{1}^{{-1}}Q_{1}^{T}y\quad,\quad||y-X\widehat{\beta}||^{2}=||Q_{2}y||^{2}.$

Wiemy z twierdzenia 7.1, że $Qy$ ma rozkład normalny. Policzmy macierz wariancji tego wektora losowego:

$\text{Var}(Q^{T}y)=Q^{T}\text{Var}(y)Q=\sigma^{2}I_{n}.$ (8.9)

Wiemy zatem, że wektory losowe $Q_{1}^{T}y$ i $Q_{2}^{T}y$ są nieskorelowane. Dla rozkładu normalnego brak korelacji jest równoważny niezależności.
Rozkład $\chi^{2}$ z $n-p$ stopniami swobody to suma $n-p$ kwadratów niezależnych zmiennych losowych o rozkładzie standardowym normalnym. Udowodnimy, że $\frac{\widehat{\sigma}^{2}(n-p)}{\sigma^{2}}=\frac{||y-X\widehat{\beta}||^{2}}{\sigma^{2}}=\frac{||Q_{2}^{T}y||^{2}}{\sigma^{2}}$ ma rozkład $\chi^{2}(n-p)$ .

Z rozkładu QR macierzy $X$ znamy wymiary macierzy $Q_{2}$ , długość wektora $Q_{2}^{T}y$ to $n-p$ . Oznaczmy:

$||Q_{2}^{T}y||^{2}=(q_{{p+1}}^{T}y)^{2}+\ldots+(q_{{n}}^{T}y)^{2}=z_{1}^{2}+\ldots+z_{{n-p}}^{2}.$

Udowodnimy, że $z_{i}$ , $i=1,\ldots,n-p$ są niezależne i mają rozkład $\mathcal{N}(0,\sigma^{2})$ .

Współrzędne wektora $Q_{2}^{T}y$ są niezależnymi zmiennymi losowymi o rozkładzie normalnym. Normalność wynika z twierdzenia 8.2, niezależność z braku korelacji (8.9). Ze wzoru 8.9 widzimy także, że wariancje $z_{i}$ są równe $\sigma^{2}$ .

Współrzędne wektora $Q_{2}^{T}y$ mają wartość oczekiwaną równą zero:

$\mathbb{E}(Q_{2}^{T}y)=Q_{2}^{T}\mathbb{E}(y)=Q_{2}^{T}X\beta=$

z wąskiego rozkładu QR macierzy $X$ ,

$=\underbrace{Q_{2}^{T}Q_{1}}_{{=0}}R_{1}=0$

z ortogonalności kolumn macierzy $Q$ .

Otrzymujemy więc:

$\left(\frac{z_{1}}{\sigma}\right)^{2}+\ldots+\left(\frac{z_{{n-p}}}{\sigma}\right)^{2}=\frac{||Q_{2}^{T}y||^{2}}{\sigma^{2}}\sim\chi^{2}(n-p),$

gdzie $\frac{z_{i}}{\sigma}$ są niezależnymi zmiennymi losowymi o rozkładzie $\mathcal{N}(0,1)$ .

∎

8.7. Test ilorazu wiarygodności (Likelihood Ratio Test) hipotez liniowych

Hipotezy liniowe przy założeniach modelu liniowego można ogólnie sformułować jako:

$\left\{\begin{array}[]{l}y=X\beta+\varepsilon\\ A\beta=0\end{array}\right.$

gdzie macierz $X$ jest wymiaru $n\times p$ , a macierz $A$ wymiaru $(p-q)\times p$ .

Przykład 8.1

Jeżeli wektor współczynników jest postaci:

$\beta=\left(\begin{array}[]{c}\beta _{1}\\ \beta _{2}\\ \beta _{3}\\ \beta _{4}\\ \end{array}\right)$

i chcemy nałożyć ograniczenie liniowe na parametry: $\beta _{2}=\beta _{3}$ , to można go zapisać postaci:

$\underbrace{\left(\begin{array}[]{cccc}0&1&-1&0\\ \end{array}\right)}_{{A}}\beta=0.$

8.7.1. LRT ogólnie

Ogólnie test ilorazu wiarygodności dotyczący parametru $\theta$ rozkładu zmiennej losowej $X$ można zapisać jako:

$H_{0}:\text{ }X\sim f_{{\theta}}(x)\text{ , }\theta\in\Theta _{0};$

$H_{1}:\text{ }X\sim f_{{\theta}}(x)\text{ , }\theta\in\Theta _{1}=\Theta\setminus\Theta _{0};$

gdzie $f_{{\theta}}(x)$ oznacza gęstość rozkładu zmiennej $X$ zależącą od parametru $\theta$ .

Statystyka testowa wyraża się wzorem:

$\lambda(x)=\frac{\sup _{{\theta\in\Theta _{1}}}f_{{\theta}}(x)}{\sup _{{\theta\in\Theta _{0}}}f_{{\theta}}(x)}=\frac{f_{{\widehat{\theta}(x)}}(x)}{f_{{\dot{\theta}(x)}}(x)},$

gdzie:

$\widehat{\theta}\text{ to estymator największej wiarygodności dla }\theta\in\Theta _{1};$

$\dot{\theta}\text{ to estymator największej wiarygodności dla }\theta\in\Theta _{0}.$

Uwaga 8.1

Jeżeli $f_{{\theta}}(x)=\prod _{{i=1}}^{n}f_{{\theta}}(x_{i})$ to:

$\log\lambda(x)=\sum _{{i=1}}^{n}\log f_{{\widehat{\theta}(x)}}(x)-\sum _{{i=1}}^{n}\log f_{{\dot{\theta}(x)}}(x).$

8.7.2. Modele zagnieżdżone

Z modelem zagnieżdżonym mamy do czynienia gdy $\Theta _{0}\subseteq\Theta$ .

Rozpatrzmy następujący problem:

$\Theta\subseteq\mathbb{R}^{p},\quad h:\mathbb{R}^{p}\rightarrow\mathbb{R}^{{p-q}},\quad\Theta _{0}=\{\theta:h(\theta)=0\}.$

Dla hipotez liniowych mamy:

$h(\theta)=A\theta=0,$

wtedy typowo $\Theta _{0}\subseteq\mathbb{R}^{q}$ , skąd możemy zapisać:

$\sup _{{\theta\in\Theta _{1}}}f_{{\theta}}(x)=\sup _{{\theta\in\Theta}}f_{{\theta}}(x).$

Dzięki takiemu zapisowi upraszcza się wzór na statystykę testową LRT:

$\lambda(x)=\frac{\sup _{{\theta\in\Theta}}f_{{\theta}}(x)}{\sup _{{\theta\in\Theta _{0}}}f_{{\theta}}(x)}.$

Twierdzenie 8.5 (Asymptotyczny rozkład LRT)

Przy założeniach: $\Theta\subseteq\mathbb{R}^{p}$ otwarty, $f_{{\theta}}(x)$ regularna rodzina gęstości, $h:\mathbb{R}^{p}\rightarrow\mathbb{R}^{{p-q}}$ funkcja gładka, $\Theta _{0}=\{\theta:h(\theta)=0\}$ :

$\forall\theta\in\Theta _{0}\quad\forall t\quad\mathbb{P}_{{\theta}}\left(2\log\lambda(x)\leq t\right)\xrightarrow[n\rightarrow\infty]{}F_{{\chi^{2}(p-q)}}(t),$

gdzie $F_{{\chi^{2}(p-q)}}$ oznacza dystrybuantę rozkładu $\chi^{2}$ o $p-q$ stopniach swobody.

8.7.3. LRT w modelu liniowym

Wracamy teraz do modelu linowego i zakładamy normalność rozkładu $\varepsilon$ :

$\left\{\begin{array}[]{l}y=X\beta+\varepsilon\\ A\beta=0\end{array}\right.$

gdzie $X$ ma wymiary $n\times p$ , $A$ wymiary $(p-q)\times p$ ;

$\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n}).$

Dla tak sformułowanego zadania wiemy, że rozkład danych $y$ jest normalny i wyraża się wzorem:

$f_{{\theta}}(y)=\frac{1}{(2\pi\sigma^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{||y-X\beta||^{2}}{2\sigma^{2}}\right],$

gdzie $\theta=(\beta _{1},\ldots,\beta _{p},\sigma^{2})$ .

Stwierdzenie 8.7

Statystyka testowa testu ilorazu wiarygodności dla $H_{0}:A\beta=0$ jest równa:

$\lambda(y)=\left[\frac{\dot{\sigma}^{2}}{\widehat{\sigma}^{2}}\right]^{{\frac{n}{2}}}=\left[\frac{||y-X\dot{\beta}||^{2}}{||y-X\widehat{\beta}||^{2}}\right]^{{\frac{n}{2}}},$

(8.10)

gdzie:

$\widehat{\theta}=(\widehat{\beta},\widehat{\sigma}^{2})\text{ estymatory największej wiarygodności bez ograniczeń (na zbiorze $\Theta$);}$

$\dot{\theta}=(\dot{\beta},\dot{\sigma}^{2})\text{ estymatory największej wiarygodności z ograniczeniami}$

$\text{($A\beta=0$, na zbiorze $\Theta _{0}$)}.$

$\lambda(y)=\frac{f_{{\widehat{\theta}(y)}}(y)}{f_{{\dot{\theta}(y)}}(y)}=\frac{\frac{1}{(\sqrt{2\pi})^{n}}\frac{1}{(\widehat{\sigma}^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{1}{2}\frac{||y-X\widehat{\beta}||^{2}}{\widehat{\sigma}^{2}}\right]}{\frac{1}{(\sqrt{2\pi})^{n}}\frac{1}{(\dot{\sigma}^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{1}{2}\frac{||y-X\dot{\beta}||^{2}}{\dot{\sigma}^{2}}\right]}=\clubsuit$

korzystając z postać estymatora największej wiarygodności dla parametru $\sigma^{2}$ w modelu liniowym, możemy zapisać:

$\widehat{\sigma}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n}\quad,\quad\dot{\sigma}^{2}=\frac{||y-X\dot{\beta}||^{2}}{n};$

pdstawiając otrzymujemy:

$\clubsuit=\frac{\frac{1}{(\widehat{\sigma}^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{1}{2}\frac{\widehat{\sigma}^{2}}{n\widehat{\sigma}^{2}}\right]}{\frac{1}{(\dot{\sigma}^{2})^{{\frac{n}{2}}}}\exp\left[-\frac{1}{2}\frac{\dot{\sigma}^{2}}{n\dot{\sigma}^{2}}\right]}=\left[\frac{\dot{\sigma}^{2}}{\widehat{\sigma}^{2}}\right]^{{\frac{n}{2}}}.$

∎

Stwierdzenie 8.8

Statystyka testowa

$\lambda(y)=\left[\frac{\dot{\sigma}^{2}}{\widehat{\sigma}^{2}}\right]^{{\frac{n}{2}}}$

jest równoważna statystyce:

$F=\frac{(R_{0}-R)/(p-q)}{R/(n-p)},$

gdzie $R_{0}=||y-X\dot{\beta}||^{2}$ , $R=||y-X\widehat{\beta}||^{2}$ .

Ze wzoru 8.10 widzimy, że:

$\lambda(y)=\left(\frac{R_{0}}{R}\right)^{{\frac{n}{2}}}.$

Statystyka $\lambda$ jako iloraz norm dwóch wektorów, jest nieujemna, a dla $p<n$ dodatnia z dodatniości $R$ i $R_{0}$ . Istnieje rosnące przekształcenie $\lambda$ w $F$ dla $\lambda>0$ , więc statystyki są rónoważne.

∎

Twierdzenie 8.6

Statystyka $F$ przy $p\ll n$ ma rozkład $\mathcal{F}$ -Snedecora:

$F=\frac{(R_{0}-R)/(p-q)}{R/(n-p)}\sim\mathcal{F}(p-q,n-p).$

Zmieńmy oznaczenia dotyczące macierzy planu. Macierz $X=(x_{1},\ldots,x_{p})$ gdzie $x_{i}$ będą oznaczać kolumny macierzy, zwane predyktorami. Możemy wtedy zapisac:

$y=X\beta+\varepsilon=x_{1}\beta _{1}+\ldots+x_{p}\beta _{p}+\varepsilon.$

Wiemy, że:

$X\widehat{\beta}\in\mathcal{L}=\{ X\beta=x_{1}\beta _{1}+\ldots+x_{p}\beta _{p}\text{, }\beta\in\mathbb{R}^{p}\};$

$X\dot{\beta}\in\mathcal{L}_{0}=\{ X\beta\text{ : }A\beta=0\};$

gdzie

$\mathcal{L}_{0}\subseteq\mathcal{L}\subseteq\mathbb{R}^{n},$

przestrzenie $\mathcal{L}_{0}$ i $\mathcal{L}_{0}$ są przestrzeniami liniowymi o wymiarach:

$\underbrace{\text{dim}(\mathcal{L}_{0})}_{{=q}}<\underbrace{\text{dim}(\mathcal{L})}_{{=p}}\ll n.$

Ortogonalizujemy bazę przestrzeni $\mathcal{L}_{0}$ , uzupełniamy do bazy $\mathcal{L}$ , a następnie do bazy $\mathbb{R}^{n}$ . Oznaczmy:

$V=(\underbrace{\underbrace{v_{1},\ldots,v_{q}}_{{\text{baza }\mathcal{L}_{0}}},v_{{q+1}},\ldots,v_{p}}_{{\text{baza }\mathcal{L}}},v_{{p+1}},\ldots,v_{n}).$

oraz:

$Z=V^{T}y,\quad\widehat{Z}=V^{T}X\widehat{\beta},\quad\dot{Z}=V^{T}X\dot{\beta}.$

Zauważmy, że wektory te są postaci:

$Z=\left(\begin{array}[]{c}z_{1}\\ \ldots\\ z_{q}\\ \ldots\\ z_{p}\\ \ldots\\ \ldots\\ z_{n}\\ \end{array}\right)\quad\widehat{Z}=\left(\begin{array}[]{c}\widehat{z}_{1}\\ \ldots\\ \widehat{z}_{q}\\ \ldots\\ \widehat{z}_{p}\\ 0\\ \ldots\\ 0\\ \end{array}\right)\quad\dot{Z}=\left(\begin{array}[]{c}\dot{z}_{1}\\ \ldots\\ \dot{z}_{q}\\ 0\\ \ldots\\ \ldots\\ \ldots\\ 0\\ \end{array}\right).$

Możemy wtedy zapisać:

$R=||y-X\widehat{\beta}||^{2}=$

ponieważ mnożenie wektora przez macierz ortogonalną nie zmienia jego normy,

$=||V^{T}(y-X\widehat{\beta})||^{2}=||Z-\widehat{Z}||^{2};$

$R=||y-X\dot{\beta}||^{2}=||V^{T}(y-X\dot{\beta})||^{2}=||Z-\dot{Z}||^{2}.$

Najlepszymi dopasowaniami $\dot{Z}$ do $Z$ oraz $\widehat{Z}$ do $Z$ minimalizującymi błędy średniokwadratowe $R_{0}$ i $R$ są:

$\dot{Z}=\left(\begin{array}[]{c}\dot{z}_{1}=z_{1}\\ \ldots\\ \dot{z}_{q}=z_{q}\\ 0\\ \ldots\\ \ldots\\ 0\\ \end{array}\right)\quad\widehat{Z}=\left(\begin{array}[]{c}\widehat{z}_{1}=z_{1}\\ \ldots\\ \ldots\\ \widehat{z}_{p}=z_{p}\\ 0\\ \ldots\\ 0\\ \end{array}\right)$

Stąd:

$R_{0}=z_{{q+1}}^{2}+\ldots+z_{p}^{2}+\ldots+z_{n}^{2};$

$R=z_{{p+1}}^{2}+\ldots+z_{n}^{2}.$

Ponieważ założyliśmy rozkład normalny dla $\varepsilon$ , możemy zapisać:

$y\sim\mathcal{N}(X\beta,\sigma^{2}I_{n}),$

a także, ponieważ $V$ jest macierzą ortogonalną:

$Z=V^{T}y\sim\mathcal{N}(V^{T}X\beta,\sigma^{2}I_{n}).$

(8.11)

Współrzędne wektora $Z$ : $z_{1},\ldots,z_{n}$ mają więc rozkłady normalne i są niezależne (bo nieskorelowane). Co więcej, przy założeniu hipotezy zerowej, $X\beta\in\mathcal{L}_{0}$ , czyli jest postaci:

$X\beta=\left(\begin{array}[]{c}w_{1}\\ \ldots\\ w_{q}\\ 0\\ \ldots\\ 0\\ \end{array}\right)$

w bazie $V$ . Ze wzoru 8.11, $\mathbb{E}(Z)=V^{T}X\beta$ , czyli:

$\mathbb{E}(z_{{q+1}})=\ldots=\mathbb{E}(z_{{n}})=0.$

Widzimy teraz, że wyrażenie:

$\frac{R_{0}-R}{\sigma^{2}}=\left(\frac{z_{{q+1}}}{\sigma}\right)^{2}+\ldots+\left(\frac{z_{{p}}}{\sigma}\right)^{2}$

ma rozkład $\chi^{(}p-q)$ , a wyrażenie:

$\frac{R}{\sigma^{2}}=\left(\frac{z_{{p+1}}}{\sigma}\right)^{2}+\ldots+\left(\frac{z_{{n}}}{\sigma}\right)^{2}$

rozkład $\chi^{(}n-p)$ oraz oba wyrażenia są od siebie ziezależne. Wróćmy do postaci statystyki F:

$F=\frac{(R_{0}-R)/(p-q)}{R/(n-p)}=\frac{(z_{{q+1}}^{2}+\ldots+z_{p}^{2})/(p-q)}{(z_{{p+1}}^{2}+\ldots+z_{n}^{2})/(n-p)}$

ma więc rozkład $\mathcal{F}(p-q,n-p)$ .

∎

Uwaga 8.2

Zauważmy ciekawą własność bazującą na dowodzie twierdzenia: dla $H_{0}:$ $\beta _{1}=\ldots=\beta _{p}=0$ przy modelu postaci:

$y=\left(\begin{array}[]{cccc}1&&&\\ 1&&\star&\\ \ldots&&&\\ 1&&&\\ \end{array}\right)\left(\begin{array}[]{c}\beta _{0}\\ \beta _{1}\\ \ldots\\ \beta _{p}\\ \end{array}\right)+\varepsilon,$

możemy zapisać:

$||y-\underbrace{\overline{y}1_{n}}_{{=X\dot{\beta}}}||^{2}=||y-\underbrace{\widehat{y}}_{{=X\widehat{\beta}}}||^{2}+||\widehat{y}-\overline{y}1_{n}||^{2},$

gdzie $\overline{y}$ jest średnią arytmetyczną z obserwacji w wektorze $y$ .

Wniosek 8.7

Testowanie hipotez o istotności współczynników (testowanie hipotez, czy kolejne grupy $\beta _{i}$ są równe zeru) służy wyborowi modelu (podzbioru zmiennych objaśniających $x_{1},\ldots,x_{p}$ ).

8.8. Popularne kryteria wyboru modelu – kryteria informacyjne

W poprzednim rozdziale zostało opisane testowanie hipotez o istotności współczynników jako sposób wyboru modelu. Wybór predyktorów można także oprzeć na minimalizacji estymatora błędu predykcji wyliczonego na podstawie kroswalidacji leave-one-out (8.7). Opiszemy teraz jeszcze inną metodę wyboru zmiennych objaśniających bazującą na tak zwanych kryteriach informacyjnych postaci:

$\text{kryterium }=\text{ }-2\cdot\text{logwiarygodność }+\text{ kara za złożoność modelu},$

które obliczane są dla każdego modelu (dla każdego podzbioru $p$ predyktorów) i wybierany jest ten minimalnej wartości kryterium. Dwa popularne kryteria informacyjne:

Akaike Information Criterion (AIC):

$AIC=\text{ }-2\cdot\text{loglik }+2\cdot p,$

gdzie $p$ to liczba zmiennych objaśniających w modelu.
Bayes Information Criterion (BIC):

$BIC=\text{ }-2\cdot\text{loglik }+\log(n)\cdot p,$

gdzie $n$ to liczba obserwacji w modelu.

Przy założeniach modelu liniowego i normalności rozkładu $\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{n})$ , kryteria przyjmują łatwiejszą postać:

Przy znanym $\sigma^{2}$ :

$AIC=\frac{||y-\widehat{y}||^{2}}{\sigma^{2}}+2\cdot p;$

$BIC=\frac{||y-\widehat{y}||^{2}}{\sigma^{2}}+\log(n)\cdot p;$
Przy nieznanym $\sigma^{2}$ :

$AIC=n\log\left(\frac{||y-\widehat{y}||^{2}}{n}\right)+2\cdot p;$

$BIC=n\log\left(\frac{||y-\widehat{y}||^{2}}{n}\right)+\log(n)\cdot p;$

8.9. Model logistyczny – przykład uogólnionego modelu liniowego

Modelu logistycznego używa się do objaśniania zmiennej binarnej, czyli przyjmującej wartości ze zbioru $\{ 0,1\}$ . Poprzednio zakładaliśmy:

$y_{i}\sim\mathcal{N}(x_{i}^{T}\beta,\sigma^{2}),\quad y_{i}\text{ niezależne }i=1,\ldots,n,$

gdzie wektor $x_{i}$ oznacza wiersz macerzy planu.

Teraz będziemy zakładać rozkład:

$y_{i}\sim\underbrace{\mathcal{B}(1,p(x_{i}))}_{{\text{rozkład Bernoulliego}}}=\mathcal{B}\left(1,\frac{\exp(x_{i}^{T}\beta)}{1+\exp(x_{i}^{T}\beta)}\right),\quad y_{i}\text{ niezależne }i=1,\ldots,n,$

gdzie postać funkcji $p(x_{i})$ można tłumaczyć tym, że prawdopodobieństwo powinno przyjmować wartości z przedziału $[0,1]$ . Parametry modelu ( $\beta$ ) estymuje się metodą największej wiarygodności, gdzie funkcja wiarygodności jest równa:

$L(\beta)=\prod _{{i=1}}^{n}\left(p(x_{i})\right)^{{y_{i}}}\left(1-p(x_{i})\right)^{{1-y_{i}}}=\prod _{{i=1}}^{n}\frac{\exp(y_{i}x_{i}^{T}\beta)}{1+\exp(x_{i}^{T}\beta)}.$

Logarytm funkcji wiarygodności maksymalizuje się numerycznie aby otrzymać estymatory $\widehat{\beta}$ . Predykcję w modelu można oprzeć na klasyfikatorze:

$t=x_{{n+1}}^{T}\widehat{\beta},$

gdzie $x_{{n+1}}$ jest wektorem nowych obserwacji. Przewidywany na podstawie modelu $\widehat{y}_{{n+1}}$ to wtedy:

$\widehat{y}_{{n+1}}=\left\{\begin{array}[]{ll}1,&\hbox{$t\geq 0$ (wtedy $p(x_{{n+1}})\geq\frac{1}{2}$);}\\ 0,&\hbox{$t<0$ (wtedy $p(x_{{n+1}})<\frac{1}{2}$).}\end{array}\right.$

8.10. Przykłady w programie R

Model liniowy:

regresja liniowa z diagnostyką dla danych Bodyfat: http://www.mimuw.edu.pl/~pokar/StatystykaII/DANE/bodyfat.R
regresja liniowa z diagnostyką dla danych Samochody: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/samochodyNowe.R
regresja liniowa dla danych Iris: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/lm.R
regresja liniowa dla danych Samochody: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/samochody.R
porównanie metody najmniejszych kwadratów i sieci neuronowych: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/crossValRegr.R

Logit (model logistyczny):

estymacja parametrów: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/logit.R
estymacja parametrów i rysowanie wyników: http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/Orings.R

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Statystyka II wykłady

Zagadnienia

8. Modele liniowe

8.1. Predykcja cechy ciągłej

Twierdzenie 8.1

Problem 8.1

8.2. Metoda najmniejszych kwadratów (mnk)

Twierdzenie 8.2

Wniosek 8.1

8.3. Inne wyprowadzenie estymatora najmniejszych kwadratów

Twierdzenie 8.3

Stwierdzenie 8.1

8.4. Estymatory metody największej wiarygodności prametrów modelu liniowego

Wniosek 8.2

8.5. Kolejne własności estymatorów mnk

8.5.1. Wartość oczekiwana i wariancja estymatora

8.5.2. Dopasowanie jako rzut ortogoanlny na przestrzeń rozpiętą przez kolumny macierzy

Definicja 8.1

Stwierdzenie 8.2

Stwierdzenie 8.3

Wniosek 8.3

8.5.3. Nieobciążony estymator parametru

Stwierdzenie 8.4

Wniosek 8.4

Stwierdzenie 8.5

Stwierdzenie 8.6

8.5.4. Model z większą liczbą parametrów nie musi być lepiej dopasowany dla nowych danych

Definicja 8.2

Wniosek 8.5

8.5.5. Kroswalidacja leave-one-out

Wniosek 8.6

8.6. Model liniowy przy założeniu normalności

Twierdzenie 8.4 (Fishera)

8.7. Test ilorazu wiarygodności (Likelihood Ratio Test) hipotez liniowych

Przykład 8.1

8.7.1. LRT ogólnie

Uwaga 8.1

8.7.2. Modele zagnieżdżone

Twierdzenie 8.5 (Asymptotyczny rozkład LRT)

8.7.3. LRT w modelu liniowym

Stwierdzenie 8.7

Stwierdzenie 8.8

Twierdzenie 8.6

Uwaga 8.2

Wniosek 8.7

8.8. Popularne kryteria wyboru modelu – kryteria informacyjne

8.9. Model logistyczny – przykład uogólnionego modelu liniowego

8.10. Przykłady w programie R

8.5.1. Wartość oczekiwana i wariancja estymatora $\widehat{\beta}$

8.5.2. Dopasowanie $\widehat{y}$ jako rzut ortogoanlny $y$ na przestrzeń rozpiętą przez kolumny macierzy $X$

8.5.3. Nieobciążony estymator parametru $\sigma^{2}$