# Rozdział 3 Model liniowy

Definicja 3.1 Hiperpłaszczyzna \(H\) wymiaru \(d \geq 0\) w przestrzeni \(\mathbf{R}^p\) o kierunku podprzestrzeni \(U\) i przechodząca przez punkt \(A \in \mathbf{R}^p\) \[ H=H(U,A) \stackrel{def}{=}U+A,\\ dim(U) = d \]

Przykład 3.1 Gdy \(U=\{\mathbf{0}\}\) to hiperpłaszczyzna przechodząca przez \(A\) jest punktem \(A\). Przyjmiemy, ze wymiar takiej hiperpłaszczyzny jest równy 0.
Gdy \(U\) jest jednowymiarową podprzestrzenią \(\mathbf{R}^p\) o wektorze bazowym \(B\) to hiperpłaszczyzna wymiaru 1 przechodząca przez \(A\) jest prostą o kierunku wektora \(B\) przechodzącą przez punkt \(A\)

Dla macierzy \(X\) i hiperpłaszczyzny \(H\) Oznaczenie \(X \subset H\) jest równoważne spełnieniu warunku \[ X_i \in H \;\; i=1,2,\ldots,n \]

Definicja 3.2 Model liniowy
Niech \(X\) macierz danych \(n \times p\). Modelem liniowym macierzy \(X\) o wymiarze \(d\) jest macierz \(Y \subset H\) wymiaru \(n \times p\), dla pewnej hiperpłaszczyzny \(H\) wymiaru \(0 \leq d \leq p\).
Definicja 3.3 Najlepszy model liniowy wymiaru \(d\)
Niech \(Y\) będzie dowolnym modelem liniowym macierzy \(X\) o wymiarze \(d\).
Model \(Y^*\) jest najlepszym modelem liniowym macierzy \(X\) o wymiarze \(d\) gdy spełnia warunek \[ d^2(X,Y^*)\leq d^2(X,Y) \;\; \forall Y \]
Twierdzenie 3.1 Twierdzenie Pitagorasa dla punktu
\[ \forall A \in \mathbf{R}^p\,\, d^2(X,A^n) = d^2(X,g(X)^n)+d^2(A^n,g(X)^n)=\\ J(X)+||g(X)-A||^2 \]
Wniosek 3.1 \[ \forall A \in \mathbf{R}^p \,\, d^2(X,g(X)^n)\leq d^2(X,A^n)\\ J(X) \leq d^2(X,A^n) \]
Wniosek 3.2 Najlepszym modelem 0-wymiarowym macierzy \(X\) jest \(g(X)^n\)
Najlepszym modelem \(p\)-wymiarowym macierzy \(X\) jest \(X\)
Definicja 3.4 Kwadratowy błąd względny modelu
Niech \(Y\) będzie modelem liniowym macierzy \(X\). Kwadratowy błąd względny modelu to liczba \[ RSE(Y)=\frac{d^2(X,Y)}{J(X)} \]
Twierdzenie 3.2 Niech \(Y_1\) i \(Y_2\) będą najlepszymi modelami liniowymi macierzy \(X\) o wymiarach \(0\leq d_1 \leq d_2 \leq p\). Wtedy \[ 0\leq RSE(Y_2) \leq RSE(Y_1) \leq 1 \]
Definicja 3.5 Rzut prostopadły na hiperpłaszczyznę
Niech \(H=H(U,A)\) będzie hiperpłaszczyzną o kierunku podprzestrzeni \(U\) o wymiarze \(d>0\) i przechodzącą przez punkt \(A \in \mathbf{R}^p\). Rzutem punktu \(Q \in \mathbf{R}^p\) na \(H\) jest punkt \[ P_H(Q)=P_U(Q-A)+A=P_U(Q)+(A-P_U(A)) \] gdzie \(P_U()\) jest rzutem prostopadłym na podprzestrzeń \(U\)
Wniosek 3.3 \(P_H\) jest operatorem liniowym \(\iff\) \(H\) jest podprzestrzenią liniową
Wniosek 3.4 Dla \(n \times p\) macierzy \(W\) \[ P_H(W)=W \iff W \subset H \]
Twierdzenie 3.3 Twierdzenie Pitagorasa dla hiperpłaszczyzn
Niech H będzie hiperpłaszczyzną wymiaru \(0<d \leq p\), \(X,Y\) macierzami danych wymiaru \(n \times p\), \(Y \subset H\) \[ d^2(X,Y) = d^2(X,P_H(X))+d^2(P_H(X),Y) \]
Wniosek 3.5 Niech \(Y^*\) będzie najlepszym modelem liniowym wymiaru \(d\) dla macierzy \(X\), \(Y^* \subset H^*\).
Wtedy \[ Y^*=P_{H^*}(X) \]

> Z ostatniego wniosku wynika, że budowanie najlepszego modelu \(Y^*\) jest równoważne szukaniu najlepszej hiperpłaszczyzny rzutu \(H^*\)
> Dlatego wymiennie można używać jako synonimu modelu macierzy danych \(Y^*\) i hiperpłaszczyzny rzutu \(H^*\)

Wniosek 3.6 \[ RSE(X,P_{H^*}(X))=1-\frac{J(P_{H^*}(X))}{J(X)} \]
Definicja 3.6 Współczynnik determinacji najlepszego modelu liniowego wymiaru \(d\) \[ r^2(P_{H^*}(X)) \stackrel{def}{=}\frac{J(P_{H^*}(X))}{J(X)} \]
Twierdzenie 3.4 \[ H^*=U^*+g(X) \] dla pewnej podprzestrzeni \(U^* \subset \mathbf{R}^p\)

> Budowanie najlepszego modelu wymiaru \(d\) mozna ograniczyć do szukania najlepszej podprzestrzeni wymiaru \(d\).

Twierdzenie 3.5 Najlepszy model \(d\) wymiarowy \(H^*=U^*+g(X)\) dla macierzy danych \(X\) spełnia jeden z dwóch warunków; \[ J(P_{U^*+g(X)}(X)) \geq J(P_{U+g(X)}(X)) \\ Tr(V(P_{U^*}(X^0))) \geq J(P_U(X^0)) \] dla dowolnej podprzestrzeni \(U\) wymiaru \(d\)

Dla zadanego układu współrzędnych operator rzutu na podprzestrzeń \(d\)-wymiarową można utożsamić z macierzą \(C\) wymiaru \(d \times p\), której wiersze są wektorami bazy (ortonormalnej) przestrzeni rzutu \(U\): \[ C= \begin{bmatrix} C_1^T \\ C_2^T \\ \vdots \\ C_d^T \end{bmatrix} \]

Twierdzenie 3.6 Macierz \(C^*\) odpowiadająca rzutowi na optymalną podprzestrzeń dla macierzy \(X\) o macierzy wariancji \(V=V(X)\) spełnia warunek: \[ C^*=argmax(Tr(CVC^T)) \] z warunkiem ubocznym \[ ||C_i||=1,\;i=1,2,\dots,d \]

Korzystając z metody mnożników Lagrange’a zadanie maksymalizacji z poprzedniego twierdzenia może być przedstawione jako zagadnienie maksymalizacji funkcji, gdzie zmienną jest macierz \(C\) wymiaru \(d \times p\) \[ Tr(CVC^T)- Tr(C^T \Lambda C) \] Macierz \(\Lambda = diag(\lambda_1,\lambda_2,\dots,\lambda_d)\) jest diagonalna. Wartości \(\lambda_i\) muszą być tak dobrane, aby spełniony był warunek \(||C_i||=1,\;i=1,2,\dots,d\)

Skorzystamy tu z warunku dostatecznego na istnienie ekstremum funkcji rzeczywistej z argumentem macierzowym (uogólnienie pojecia gradientu).

Definicja 3.7 Pochodna macierzowa funkcji rzeczywistej macierzy \(D\) rozmiaru \(d \times p\) jest macierzą rozmiaru \(d \times p\): \[ \left [ \frac{\partial f}{\partial D}(D) \right ]_{ij} = \frac{\partial f(D)}{\partial D_{ij}} \]
Twierdzenie 3.7 Jeżeli funkcja \(f(D)\) ma ekstremum dla macierzy \(D_0\) to \[ \left.\frac{\partial f}{\partial D}(D)\right|_{D=D_0}=0 \]
Lemat 3.1 Dla symetrycznej macierzy \(W\) \[ \frac{\partial Tr(A^TWA)}{\partial A}=2WA\\ \frac{\partial Tr(BWB^T)}{\partial B}=2BW\\ \]
Twierdzenie 3.8 Macierz \(C\) odpowiadająca rzutowi na optymalną podprzestrzeń dla macierzy \(X\) o macierzy wariancji \(V=V(X)\) (twierdzenie 3.6) spełnia równanie \[ CV=\Lambda C \]
Wniosek 3.7 Wiersze macierzy \(C\) (składowe główne) są wektorami własnymi macierzy \(V\), odpowiadającymi jej \(d\) największym wartościom własnym \(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_d\)
Wniosek 3.8 Współczynnik determinacji najlepszego modelu \(d\)-wymiarowego wyraża się wzorem; \[ r^2=\frac{\sum_{i=1}^d \lambda_i}{\sum_{i=1}^p \lambda_i} \]

## 3.1 Reifikacja modelu

Składowe główne są zbiorem nowych cech w \(d\)-wymiarowej przestrzeni najbliższej w sensie odległości Frobeniusa danym z macierzy \(X\). Interesujące jest uzyskanie praktycznej interpretacji tych nowych cech. Taki proces nazywa się reifikacją modelu.
Zazwyczaj wszystkie analizy prowadzi się dla danych scentrowanych \(X^0\), to znaczy w układzie współrzędnych, którego początek umieszcza się w środku ciężkości \(X\).
Model danych \(X\) w przestrzeni składowych głównych (ang. PCA scores) jest macierzą \(Y=X^0C^T\).
O związku między kolumnami \(X^0\) i \(Y\) świadczy macierz korelacji \(R(X,Y)\)

Propozycja 3.1 Niech \(V=V(X)=[v_{ij}]\) będzie macierzą wariancji \(X\), \(C=[c_{ij}]\) macierzą, której wiersze są składowymi głównymi \[ r_{ij} \stackrel{def}{=} [R(X,Y)]_{ij}=\sqrt{\frac{\lambda_j}{v_{ii}}}c_{ij} \]
Propozycja 3.2 \[ \sum_{j=1}^p r_{ij}^2=1 \]

Liczba \(r_{ij}\) jest korelacją między \(X^i\) a \(Y^j\), czyli cosinusem kąta miedzy nimi. Im ten kąt mniejszy (a więc korelacja większa) tym bardziej nowa zmienna \(Y^j\) jest związana z \(X^i\). Propozycja 3.2, mówi że wektor \([r_{i1}, r_{i2},\dots ,r_{ip}]\), reprezentujący zmienną \(X^i\) leży na sferze jednostkowej w układzie współrzędnych korelacyjnych między \(X^i\) a \(Y^1,Y^2,\dots,Y^p\).
Wykres w układzie dwóch pierwszych składowych korelacyjnych nazywa się kołem korelacyjnym. Pokazuje on jak “stare” zmienne są związane z dwiema najważniejszymi składowymi głównymi.

## 3.2 Przykład

Dane są związane z wynikami zawodów olimpijskich w Seulu (1988) w siedmioboju kobiet.
System punktowy przelicza wyniki 7 dyscyplin na punkty. Porównamy ten system punktów z modelem składowych głównych.

Słowniczek
hurdles = 100 m płotki
shot = pchnięcie kulą
javelin = oszczep

    load("siedmioboj.Rda")
    hm <- as.matrix(siedmioboj[,1:7])
    knitr::kable(
      hm[1:5,], booktabs = TRUE,
      caption = 'Wyniki 5 najlepszych zawodniczek.'
    )

Tabela 3.1: Wyniki 5 najlepszych zawodniczek.
| hurdles | highjump | shot | run200m | longjump | javelin | run800m
--- | --- | --- | --- | --- | --- | --- | ---
Joyner-Kersee (USA) | 12.69 | 1.86 | 15.80 | 22.56 | 7.27 | 45.66 | 128.51
John (GDR) | 12.85 | 1.80 | 16.23 | 23.65 | 6.71 | 42.56 | 126.12
Behmer (GDR) | 13.20 | 1.83 | 14.20 | 23.10 | 6.68 | 44.54 | 124.20
Sablovskaite (URS) | 13.61 | 1.80 | 15.23 | 23.92 | 6.25 | 42.78 | 132.24
Choubenkova (URS) | 13.51 | 1.74 | 14.76 | 23.93 | 6.32 | 47.46 | 127.90

    kor <- round(cor(hm), 2)
    knitr::kable(
      kor, booktabs = TRUE,
      caption = 'Korelacje między zmiennymi'
    )

Tabela 3.2: Korelacje między zmiennymi
| hurdles | highjump | shot | run200m | longjump | javelin | run800m
--- | --- | --- | --- | --- | --- | --- | ---
hurdles | 1.00 | -0.81 | -0.65 | 0.77 | -0.91 | -0.01 | 0.78
highjump | -0.81 | 1.00 | 0.44 | -0.49 | 0.78 | 0.00 | -0.59
shot | -0.65 | 0.44 | 1.00 | -0.68 | 0.74 | 0.27 | -0.42
run200m | 0.77 | -0.49 | -0.68 | 1.00 | -0.82 | -0.33 | 0.62
longjump | -0.91 | 0.78 | 0.74 | -0.82 | 1.00 | 0.07 | -0.70
javelin | -0.01 | 0.00 | 0.27 | -0.33 | 0.07 | 1.00 | 0.02
run800m | 0.78 | -0.59 | -0.42 | 0.62 | -0.70 | 0.02 | 1.00

    siedmioboj_pca <- prcomp(hm, scale = TRUE)
    knitr::kable(
      siedmioboj_pca$rotation, booktabs = TRUE, digits=4,
      caption = 'Składowe główne'
    )

Tabela 3.3: Składowe główne
| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7
--- | --- | --- | --- | --- | --- | --- | ---
hurdles | 0.4529 | -0.1579 | -0.0451 | 0.0265 | -0.0949 | -0.7833 | -0.3802
highjump | -0.3772 | 0.2481 | 0.3678 | -0.6800 | -0.0188 | -0.0994 | -0.4339
shot | -0.3631 | -0.2894 | -0.6762 | -0.1243 | -0.5117 | 0.0509 | -0.2176
run200m | 0.4079 | 0.2604 | 0.0836 | -0.3611 | -0.6498 | 0.0250 | 0.4534
longjump | -0.4562 | 0.0559 | -0.1393 | -0.1113 | 0.1843 | -0.5902 | 0.6121
javelin | -0.0754 | -0.8417 | 0.4716 | -0.1208 | -0.1351 | 0.0272 | 0.1729
run800m | 0.3750 | -0.2245 | -0.3959 | -0.6034 | 0.5043 | 0.1556 | 0.0983
> scale=TRU | E oznacza, | że używam | y danych s | tandaryzow | anych |  |

    lambdy <- siedmioboj_pca$sdev^2
    lmb_tab <- data.frame(lambda=round(lambdy,4),r2=round(cumsum(lambdy)/7,2))
    lmb_tab <- t(lmb_tab)
    colnames(lmb_tab) <- 1:7
    knitr::kable(
      lmb_tab, booktabs = TRUE,
      caption = 'Lambdy i współczynniki determinacji jako funkcja wymiaru modelu d'
    )

Tabela 3.4: Lambdy i współczynniki determinacji jako funkcja wymiaru modelu d
| 1 | 2 | 3 | 4 | 5 | 6 | 7
--- | --- | --- | --- | --- | --- | --- | ---
lambda | 4.4603 | 1.1943 | 0.521 | 0.4572 | 0.2453 | 0.073 | 0.049
r2 | 0.6400 | 0.8100 | 0.880 | 0.9500 | 0.9800 | 0.990 | 1.000

    punkty_pca <- -predict(siedmioboj_pca)[,1]
    pca_sco <- data.frame(score=siedmioboj$score,pca=punkty_pca)
    rownames(pca_sco) <- rownames(siedmioboj)
    knitr::kable(
      t(pca_sco[1:5,]), booktabs = TRUE,digits = 2,
      caption = 'Porównanie punktów siedmioboju i punktów z pca dla 5 zawodniczek'
    )

Tabela 3.5: Porównanie punktów siedmioboju i punktów z pca dla 5 zawodniczek
| Joyner-Kersee (USA) | John (GDR) | Behmer (GDR) | Sablovskaite (URS) | Choubenkova (URS)
--- | --- | --- | --- | --- | ---
score | 7291.00 | 6897.00 | 6858.00 | 6540.00 | 6540.00
pca | 4.12 | 2.88 | 2.65 | 1.34 | 1.36

> Znak “-” przy predict(siedmioboj_pca)[,1] został wybrany tak, aby zwrot wektora pca i score był taki sam

    require(FactoMineR)
    pca.FM  <- PCA(hm, graph = FALSE)
    knitr::kable(
      pca.FM$var$coord, booktabs = TRUE,digits = 2,
      caption = 'Korelacje zmiennych oryginalnych ze składowymi głównymi'
    )

Tabela 3.6: Korelacje zmiennych oryginalnych ze składowymi głównymi
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | Dim.5
--- | --- | --- | --- | --- | ---
hurdles | -0.96 | 0.17 | 0.03 | -0.02 | 0.05
highjump | 0.80 | -0.27 | -0.27 | 0.46 | 0.01
shot | 0.77 | 0.32 | 0.49 | 0.08 | 0.25
run200m | -0.86 | -0.28 | -0.06 | 0.24 | 0.32
longjump | 0.96 | -0.06 | 0.10 | 0.08 | -0.09
javelin | 0.16 | 0.92 | -0.34 | 0.08 | 0.07
run800m | -0.79 | 0.25 | 0.29 | 0.41 | -0.25

Składowa główna pc1 jest dodatnio i silnie skorelowana ze sportami siłowymi i z prędkościami (z czasem biegu ujemnie!). Jedynym niepasującym sportem jest rzut oszczepem, związanym z drugą składową główną. Czyli punktacja, silnie związana ze składową 1 nie uwzględnia w istocie rzutu oszczepem.