Jak przeprowadzić regresję ortogonalną (suma najmniejszych kwadratów) za pomocą PCA?

Zawsze używam lm()w R do regresji liniowej na . Ta funkcja zwraca współczynnik taki, że $y$ $x$ $\beta$

y = β x .

$y = \beta x.$

Dzisiaj dowiedziałem się o całkowitej liczbie najmniejszych kwadratów i tej princomp()funkcji (analiza głównego składnika, PCA) można użyć do jej wykonania. To powinno być dla mnie dobre (dokładniejsze). Zrobiłem kilka testów przy użyciu princomp(), takich jak:

r <- princomp( ~ x + y)

Mój problem brzmi: jak interpretować jego wyniki? Jak mogę uzyskać współczynnik regresji? Przez „współczynnik” rozumiem liczbę , której muszę użyć do pomnożenia wartości celu uzyskania liczby zbliżonej do . $\beta$ $x$ $y$

— Dail
źródło

Chwilkę, jestem trochę zmieszany. spójrz na: zoonek2.free.fr/UNIX/48_R/09.html Nazywa się to PCA (Principal Component Analysis, inaczej „regresja ortogonalna” lub „prostopadłe sumy kwadratów” lub „suma najmniejszych kwadratów”), więc myślę, że mówimy o TLS z princomp () Nie?

— Dail

Nie; to dwie różne rzeczy, zobacz artykuł w Wikipedii o PCA. Fakt, że jest tu używany, to hack (nie wiem jak dokładny, ale zamierzam to sprawdzić); dlatego złożona ekstrakcja współczynników.

Powiązane pytanie: stats.stackexchange.com/questions/2691/…, a na blogu znajduje się jedna z odpowiedzi: cerebralmastication.com/2010/09/...

— Jonathan

Zwykłe najmniejsze kwadraty vs. suma najmniejszych kwadratów

Rozważmy najpierw najprostszy przypadek tylko jednej zmiennej predykcyjnej (niezależnej) . Dla uproszczenia, niech i są wyśrodkowane, tzn. Przecięcie jest zawsze zerowe. Różnica między standardową regresją OLS i „ortogonalną” regresją TLS jest wyraźnie pokazana na tej (dostosowanej przeze mnie) liczbie z najpopularniejszej odpowiedzi w najpopularniejszym wątku na PCA: $x$ $x$ $y$

OLS vs TLS

OLS dopasowuje się do równania przez minimalizowanie kwadratów odległości pomiędzy obserwowanymi wartościami i przewidywanych wartości . TLS pasuje do tego samego równania, minimalizując kwadratowe odległości między punktami i ich rzut na linię. W tym najprostszym przypadku linia TLS jest po prostu pierwszym głównym składnikiem danych 2D. Aby znaleźć , wykonaj PCA w punktach , tj. Skonstruuj macierz kowariancji i znajdź pierwszy wektor własny $y=\beta x$ $y$ $\hat y$ $(x,y)$ $\beta$ $(x,y)$ $2\times 2$ $\boldsymbol \Sigma$ $\mathbf v = (v_x, v_y)$ ; następnie . $\beta = v_y/v_x$

W Matlabie:

 v = pca([x y]);    //# x and y are centered column vectors
 beta = v(2,1)/v(1,1);

W R:

 v <- prcomp(cbind(x,y))$rotation
 beta <- v[2,1]/v[1,1]

Nawiasem mówiąc, to wydajność prawidłowe nachylenie nawet jeśli i nie wyśrodkowany (ponieważ wbudowanej funkcji automatycznego wykonywania centrowania PCA). Aby odzyskać przechwycenie, oblicz . $x$ $y$ $\beta_0 = \bar y - \beta \bar x$

OLS vs. TLS, regresja wielokrotna

Biorąc pod uwagę zmienną zależną i wiele zmiennych niezależnych (ponownie wszystkie wyśrodkowane dla uproszczenia), regresja pasuje do równaniaOLS dopasowuje się, minimalizując do kwadratu błędy między obserwowanymi wartościami a wartościami przewidywanymi . TLS dopasowuje się, minimalizując kwadratowe odległości między zaobserwowanymi punktami a najbliższymi punktami na płaszczyźnie regresji / hiperpłaszczyźnie. $y$ $x_i$

y = β_{1} x_{1} + \dots + β_{p} x_{p} .

$y= \beta_1 x_1 + \ldots + \beta_p x_p.$

y

$y$

\hat{y}

$\hat y$

(x, y) \in R^{p + 1}

$(\mathbf x, y)\in\mathbb R^{p+1}$

Zauważ, że nie ma już „linii regresji”! Powyższe równanie określa hiperpłaszczyznę : jest to płaszczyzna 2D, jeśli istnieją dwa predyktory, hiperpłaszczyzna 3D, jeśli istnieją trzy predyktory itp. Tak więc powyższe rozwiązanie nie działa: nie możemy uzyskać rozwiązania TLS, biorąc tylko pierwszy komputer (który jest linia). Mimo to rozwiązanie można łatwo uzyskać za pomocą PCA.

Tak jak poprzednio, PCA wykonuje się na punktach . Daje to wektory w kolumnach . Pierwsze wektory zdefiniować wymiarową hiperpłaszczyznę że musi; ostatni (liczba ) wektor własny jest do niego ortogonalny. Pytanie brzmi, jak przekształcić podstawę podanej przez pierwszych wektorów własnych do współczynników. $(\mathbf x, y)$ $p+1$ $\mathbf V$ $p$ $p$ $\mathcal H$ $p+1$ $\mathbf v_{p+1}$ $\mathcal H$ $p$ $\boldsymbol \beta$

Zauważ, że jeśli ustawimy dla wszystkich i tylko , wtedy , tj. Wektor leży w hiperpłaszczyzna . Z drugiej strony wiemy, że jest do niego ortogonalny. Czyli ich iloczyn iloczynu musi wynosić zero: $x_i=0$ $i \ne k$ $x_k=1$ $\hat y=\beta_k$

(0, \dots, 1, \dots, β_{k}) \in H

$(0,\ldots, 1, \ldots, \beta_k) \in \mathcal H$

H

$\mathcal H$

v_{p + 1} = (v_{1}, \dots, v_{p + 1}) ⊥ H

$\mathbf v_{p+1}=(v_1, \ldots, v_{p+1}) \:\bot\: \mathcal H$

v_{k} + β_{k} v_{p + 1} = 0 \Rightarrow β_{k} = - v_{k} / v_{p + 1} .

$v_k + \beta_k v_{p+1}=0 \Rightarrow \beta_k = -v_k/v_{p+1}.$

W Matlabie:

 v = pca([X y]);    //# X is a centered n-times-p matrix, y is n-times-1 column vector
 beta = -v(1:end-1,end)/v(end,end);

W R:

 v <- prcomp(cbind(X,y))$rotation
 beta <- -v[-ncol(v),ncol(v)] / v[ncol(v),ncol(v)]

Ponownie, to otrzymując odpowiednie terenie, nawet w przypadku i nie wycentrowany (ponieważ wbudowanej funkcji automatycznego wykonywania centrowania PCA). Aby odzyskać przechwytywanie, oblicz . $x$ $y$ $\beta_0 = \bar y - \bar {\mathbf x} \boldsymbol \beta$

W ramach kontroli poczytalności zauważ, że to rozwiązanie pokrywa się z poprzednim w przypadku tylko jednego predyktora . Rzeczywiście, wówczas przestrzeń jest 2D, a zatem biorąc pod uwagę, że pierwszy wektor własny PCA jest ortogonalny do drugiego (ostatniego), . $x$ $(x,y)$ $v^{(1)}_y/v^{(1)}_x=-v^{(2)}_x/v^{(2)}_y$

Rozwiązanie w formie zamkniętej dla TLS

Nieoczekiwanie okazuje się, że istnieje równanie o zamkniętej formie dla . Poniższy argument pochodzi z książki Sabine van Huffel „Suma najmniejszych kwadratów” (sekcja 2.3.2). $\boldsymbol \beta$

Niech i będą wyśrodkowanymi macierzami danych. Ostatni wektor własny PCA jest wektorem własnym macierzy kowariancji o wartości własnej . Jeśli jest to wektor własny, to tak samo jest . Zapisywanie równania wektora własnego: $\mathbf X$ $\mathbf y$ $\mathbf v_{p+1}$ $[\mathbf X\: \mathbf y]$ $\sigma^2_{p+1}$ $-\mathbf v_{p+1}/v_{p+1} = (\boldsymbol \beta\:\: -1)^\top$

(\begin{matrix} X^{⊤} X & X^{⊤} y \\ y^{⊤} X & y^{⊤} y \end{matrix}) (\begin{matrix} β \\ - 1 \end{matrix}) = σ_{p + 1}^{2} (\begin{matrix} β \\ - 1 \end{matrix}),

$\left(\begin{array}{c}\mathbf X^\top \mathbf X & \mathbf X^\top \mathbf y\\ \mathbf y^\top \mathbf X & \mathbf y^\top \mathbf y\end{array}\right) \left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right) = \sigma^2_{p+1}\left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right),$ i obliczając produkt po lewej, natychmiast otrzymujemy, że co mocno przypomina znane wyrażenie OLS

β_{T L S} = (X^{⊤} X - σ_{p + 1}^{2} I)^{- 1} X^{⊤} y,

$\boldsymbol \beta_\mathrm{TLS} = (\mathbf X^\top \mathbf X - \sigma^2_{p+1}\mathbf I)^{-1} \mathbf X^\top \mathbf y,$

β_{O L S} = (X^{⊤} X)^{- 1} X^{⊤} y .

$\boldsymbol \beta_\mathrm{OLS} = (\mathbf X^\top \mathbf X)^{-1} \mathbf X^\top \mathbf y.$

Wieloczynnikowa regresja wielokrotna

Tę samą formułę można uogólnić na przypadek wielowymiarowy, ale nawet zdefiniowanie działania TLS wielowymiarowego wymagałoby pewnej algebry. Zobacz Wikipedia na TLS . Wielowymiarowa regresja OLS jest równoważna wiązce jednoczynnikowych regresji OLS dla każdej zmiennej zależnej, ale w przypadku TLS tak nie jest.

— ameba mówi Przywróć Monikę
źródło

Nie znam R, ale nadal chciałem udostępnić fragmenty R do wykorzystania w przyszłości. Jest tu wielu ludzi biegle posługujących się językiem R. W razie potrzeby możesz edytować moje fragmenty! Dziękuję Ci.

— ameba mówi Przywróć Monikę

(0, \dots, 1, \dots, β_{k})

$(0,\ldots, 1, \ldots, \beta_k)$

x_{i}

$x_i$

x_{k} = 1

$x_k=1$

y = \sum β_{j} x_{j}

$y=\sum \beta_j x_j$

y = β_{k} \cdot 1 = β_{k}

$y=\beta_k\cdot 1 = \beta_k$

(0, \dots, 1, \dots β_{k})

$(0,\ldots, 1, \ldots \beta_k)$

y = \sum β_{j} x_{j}

$y=\sum \beta_j x_j$

— ameba mówi Przywróć Monikę

Wydaje mi się, że źle odczytałem tę część, ale teraz jest jasne. Dziękuję również za wyjaśnienie.

— JohnK

W R możesz preferować „ wektory własne (cov (cbind (x, y)))) $ ” niż „prcomp (cbind (x, y)) $ rotacja”, ponieważ to pierwsze jest znacznie szybsze dla większych wektorów.

— Thomas Browne,

Oparta na naiwnej realizacji GNU Octave znaleźć tutaj , coś takiego może (przymrużeniem oka, że jest późno) praca.

tls <- function(A, b){

  n <- ncol(A)
  C <- cbind(A, b)

  V <- svd(C)$v
  VAB <- V[1:n, (n+1):ncol(V)]
  VBB <- V[(n+1):nrow(V), (n+1):ncol(V)]
  return(-VAB/VBB)
}

— kaszmir
źródło

princompuruchamia analizę głównych składników zamiast regresji sumy najmniejszych kwadratów. O ile mi wiadomo, nie ma funkcji R ani pakietu, który obsługuje TLS; w MethComp występuje co najwyżej regresja Deminga .
Proszę jednak potraktować to jako sugestię, że najprawdopodobniej nie jest tego warte.

Myślałem, że Deming w pakiecie MethComp to TLS - jaka jest różnica?

— mark999

Musisz podać stosunek błędów na xiy; pure TLS to optymalizuje.