Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Statystyka II – 8. Modele liniowe – MIM UW

Zagadnienia

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:

  1. \mathbb{E}(\varepsilon)=0;

  2. \text{Var}(\varepsilon)=\sigma^{2}I_{n};

  3. 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:

  1. Predykcja dla y jest równa \widehat{y}=X\widehat{\beta};

  2. ||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):

  1. \widehat{y}=X\widehat{\beta} jest rzutem ortogonalnym y na przestrzeń rozpiętą przez kolumny macierzy planu X.

  2. Nieobciążonym estymatorem parametru \sigma^{2} jest \widehat{\sigma}^{2}=\frac{||y-\widehat{y}||^{2}}{n-p}.

  3. Twierdzenie Gaussa-Markowa: estymator \widehat{\beta} jest liniowym, nieobciążonym estymatorem o najmniejszej wariancji parametru \beta (BLUE- Best Linear Unbiased Estimator).

  4. 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:

  1. 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;
  2. 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:

  1. HX=X:

    X\underbrace{(X^{T}X)^{{-1}}X^{T}X}_{{=I_{p}}}=X;
  2. 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;
  3. symetryczna, czyli H^{T}=H:

    H^{T}=(X(X^{T}X)^{{-1}}X^{T})^{T}=X(X^{T}X)^{{-1}}X^{T}=H.
  4. 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}.
  5. 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
  1. Macierz daszkowa H jest macierzą rzutu ortogonalnego na przestrzeń rozpiętą przez kolumny macierzy X.

  2. 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ą:

  1. \widehat{\beta}\sim\mathcal{N}(\beta,\sigma^{2}(X^{T}X)^{{-1}});

  2. \widehat{\beta} i \widehat{\sigma}^{2}=\frac{||y-X\widehat{\beta}||^{2}}{n-p} są niezależne;

  3. \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}.
  1. 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.

  2. 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.

  3. 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:

  1. Akaike Information Criterion (AIC):

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

    gdzie p to liczba zmiennych objaśniających w modelu.

  2. 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ć:

  1. 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;
  2. 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:

Logit (model logistyczny):

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

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.