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 = \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$

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

Rozważmy najpierw najprostszy przypadek tylko jednej zmiennej predykcyjnej (niezależnej) . Dla uproszczenia, niech x i y 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 dopasowuje się do równania przez minimalizowanie kwadratów odległości pomiędzy obserwowanymi wartościami ý i przewidywanych wartości y . TLS pasuje do tego samego równania, minimalizując kwadratowe odległości między punktami ( x , y ) 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 ( x , y ) , tj. Skonstruuj macierz kowariancji 2 × 2 i znajdź pierwszy wektor własny$y=\beta x$$y$$\hat y$$(x,y)$$\beta$$(x,y)$$2\times 2$v = ( v x , v y ) β = v y / v $\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 .y β 0 = ˉ y - β ˉ $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.x i y = β 1 x 1 + … + β p x p . y y ( x , y ) ∈ R s + $y$$x_i$

$$y= \beta_1 x_1 + \ldots + \beta_p x_p.$$ $y$$\hat y$$(\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.p + 1 V p p H p + 1 v p + 1 H p $(\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:I ≠ k x k = 1 y = β K ( 0 , ... , 1 , ... , β k ) ∈ H H V P + 1 = ( V, 1 , ... , v p + 1$x_i=0$$i \ne k$$x_k=1$$\hat y=\beta_k$

$$(0,\ldots, 1, \ldots, \beta_k) \in \mathcal H$$ $\mathcal H$v k + β k v p + 1 = 0 ⇒ β k = - v k / v p + 1 $$\mathbf v_{p+1}=(v_1, \ldots, v_{p+1}) \:\bot\: \mathcal H$$ $$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 .y β 0 = ˉ y - ˉ x$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 , y ), V ( 1 ) Y / V ( 1 ) x = - V ( 2 ) x / v ( 2 )$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$

$$\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 $$\boldsymbol \beta_\mathrm{TLS} = (\mathbf X^\top \mathbf X - \sigma^2_{p+1}\mathbf I)^{-1} \mathbf X^\top \mathbf 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.

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)
}

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.