Częściowa regresja najmniejszych kwadratów w R: dlaczego PLS na znormalizowanych danych nie jest równoważny maksymalizacji korelacji?

Jestem bardzo nowy w częściowych najmniejszych kwadratach (PLS) i staram się zrozumieć wynik funkcji R plsr()w plspakiecie. Symulujmy dane i uruchom PLS:

library(pls)
n <- 50
x1 <- rnorm(n); xx1 <- scale(x1) 
x2 <- rnorm(n); xx2 <- scale(x2)
y <- x1 + x2 + rnorm(n,0,0.1); yy <- scale(y)
p <- plsr(yy ~ xx1+xx2, ncomp=1)

Spodziewałem się, że następujące liczby i $a$ $b$

> ( w <- loading.weights(p) )

Loadings:
    Comp 1
xx1 0.723 
xx2 0.690 

               Comp 1
SS loadings       1.0
Proportion Var    0.5
> a <- w["xx1",]
> b <- w["xx2",]
> a^2+b^2
[1] 1

są obliczane w celu maksymalizacji

> cor(y, a*xx1+b*xx2)
          [,1]
[1,] 0.9981291

ale nie jest tak dokładnie:

> f <- function(ab){
+ a <- ab[1]; b <- ab[2]
+ cor(y, a*xx1+b*xx2)
+ }
> optim(c(0.7,0.6), f, control=list(fnscale=-1))
$par
[1] 0.7128259 0.6672870

$value
[1] 0.9981618

Czy to błąd numeryczny, czy też źle rozumiem naturę i ? $a$ $b$

Chciałbym również wiedzieć, jakie są te współczynniki:

> p$coef
, , 1 comps

           yy
xx1 0.6672848
xx2 0.6368604

EDYCJA : Teraz widzę, co p$coefto jest:

> x <- a*xx1+b*xx2
> coef(lm(yy~0+x))
        x 
0.9224208 
> coef(lm(yy~0+x))*a
        x 
0.6672848 
> coef(lm(yy~0+x))*b
        x 
0.6368604

Więc myślę, że mam rację co do natury i . $a$ $b$

EDYCJA: W świetle komentarzy udzielonych przez @chl uważam, że moje pytanie nie jest wystarczająco jasne, dlatego podaję więcej szczegółów. W moim przykładzie jest wektor odpowiedzi i dwukolumnowa macierz predyktorów, a ja używam znormalizowanej wersji z i znormalizowanej wersji z (wyśrodkowana i podzielona przez odchylenia standardowe). Definicja pierwszych PLS składnika jest z i wybiera się w celu uzyskania wartości maksymalnej produktu wewnętrznej . $Y$ $X$ $\tilde Y$ $Y$ $\tilde X$ $X$ $t_1$ $t_1 = a \tilde X_1 + b \tilde X_2$ $a$ $b$ $\langle t_1, \tilde Y \rangle$ Zatem jest to równoważne z maksymalizacją korelacji między i , prawda? $t_1$ $Y$

r regression partial-least-squares

— Stéphane Laurent
źródło

Współczynnik regresji PLS maksymalizuje wyniki (co oblicza się jako iloczyn surowych danych z załadunkiem wektorem (ami)) kowariancji , nie korelacji (jak to ma miejsce w kanoniczna Correlation Analysis). W plstym dokumencie JSS znajduje się dobry przegląd pakietu i regresji PLS .

— chl

Ponieważ wszystkie wektory są wyśrodkowane i znormalizowane, kowariancja jest korelacją, prawda? Przepraszamy, ale artykuł JSS jest zbyt techniczny dla początkującego.

— Stéphane Laurent,

Ogólnie rzecz biorąc, istnieje asymetryczny proces deflacji (wynikający z regresji kombinacji liniowej jednego bloku na kombinację liniową drugiego), który nieco komplikuje sytuację. W tej odpowiedzi przedstawiłem schematyczny obraz . Hervé Abdi przedstawił ogólny przegląd regresji PLS, a metoda Wegelin's Survey of Partial Least Squares (PLS) jest również bardzo przydatna. W tym miejscu prawdopodobnie powinienem przekonwertować wszystkie te komentarze na odpowiedź ...

— chl

W moim przykładzie jest wektor odpowiedzi i dwukolumnowa macierz predyktorów, a ja używam znormalizowanej wersji z i znormalizowanej wersji z (wyśrodkowana i podzielona przez odchylenia standardowe). Moja definicja pierwszego komponentu PLS to z i wybranymi w celu uzyskania maksymalnej wartości iloczynu skalarnego . Czy to nie jest dobra definicja?

Y

$Y$

X

$X$

\tilde{Y}

$\tilde Y$

Y

$Y$

\tilde{X}

$\tilde X$

X

$X$

t_{1}

$t_1$

t_{1} = a {\tilde{X}}_{1} + b {\tilde{X}}_{2}

$t_1 = a \tilde X_1 + b \tilde X_2$

a

$a$

b

$b$

⟨ t_{1}, \tilde{Y} ⟩

$\langle t_1, \tilde Y \rangle$

— Stéphane Laurent,

Przepraszamy, @ Stéphane, ponieważ moje komentarze powyżej nie uwzględniały faktu, że poprosiłeś tylko o jeden komponent (więc deflacja nie odgrywa tutaj kluczowej roli). Wydaje się jednak, że twoja funkcja optymalizacji nie narzuca wektorów masy jednostkowej, tak że ostatecznie . (btw, poda więcej informacji na temat tych „współczynników”, ale najwyraźniej sam to odkryłeś).

a^{2} + b^{2} \neq 1

$a^2+b^2\neq 1$ ?coef.mvr

— chl

Regresja PLS opiera się na algorytmach iteracyjnych (np. NIPALS, SIMPLS). Twój opis głównych pomysłów jest poprawny: szukamy jednego (PLS1, jednej zmiennej odpowiedzi / wielu predyktorów) lub dwóch (PLS2, z różnymi trybami, wielu zmiennych odpowiedzi / wielu predyktorów) wektora (-ów) wag, (i ) powiedzmy, aby utworzyć kombinację liniową pierwotnej zmiennej (zmiennych), tak że kowariancja między Xu i Y (Yv, dla PLS2) jest maksymalna. Skupmy się na wyodrębnieniu pierwszej pary wag powiązanych z pierwszym składnikiem. Formalnie kryterium optymalizacji czyta W twoim przypadku jest jednoznaczne, więc sprowadza się do maksymalizacji $u$ $v$

max cov (X u, Y v) . (1)

$\max\text{cov}(Xu, Yv).\qquad (1)$

Y

$Y$

cov (X u, y) \equiv Var (X u)^{1 / 2} \times cor (X u, y) \times Var (y)^{1 / 2}, s t . ‖ u ‖ = 1.

$\text{cov}(Xu, y)\equiv \text{Var}(Xu)^{1/2}\times\text{cor}(Xu, y)\times\text{Var}(y)^{1/2},\quad st. \|u\|=1.$ Ponieważ nie zależy od , musimy zmaksymalizować . Zastanówmy się , gdzie dane są indywidualnie standaryzowane (początkowo popełniłem błąd skalowania kombinacji liniowej zamiast osobno i !), Tak że ; jednak i zależy od . Podsumowując, maksymalizacja korelacji między składnikiem utajonym a zmienną odpowiedzi nie da takich samych wyników

Var (y)

$\text{Var}(y)$

u

$u$

Var (X u)^{1 / 2} \times cor (X u, y)

$\text{Var}(Xu)^{1/2}\times\text{cor}(Xu, y)$ X=[x_1;x_2]

x_{1}

$x_1$

x_{2}

$x_2$

Var (x_{1}) = Var (x_{2}) = 1

$\text{Var}(x_1)=\text{Var}(x_2)=1$

Var (X u) \neq 1

$\text{Var}(Xu)\neq 1$

u

$u$ .

Powinienem podziękować Arthurowi Tenenhausowi, który wskazał mi właściwy kierunek.

Używanie wektorów wagi jednostkowej nie jest ograniczające, a niektóre pakiety ( pls. regressionw wersji plsgenomics , oparte na kodzie z wcześniejszego pakietu Wehrensa pls.pcr) zwracają niestandardowe wektory wagi (ale z ukrytymi składnikami wciąż o normie 1), jeśli są wymagane. Ale większość pakietów PLS zwróci znormalizowany , w tym ten, którego użyłeś, zwłaszcza te implementujące algorytm SIMPLS lub NIPALS; Znalazłem dobry przegląd obu podejść w prezentacji Barry'ego M. Wise'a, właściwości regresji częściowych najmniejszych kwadratów (PLS) i różnic między algorytmami , ale chemometrii $u$ winieta oferuje również dobrą dyskusję (str. 26–29). Szczególnie ważny jest również fakt, że większość procedur PLS (przynajmniej ta, którą znam w R) zakłada, że udostępniasz niestandardowe zmienne, ponieważ centrowanie i / lub skalowanie jest obsługiwane wewnętrznie (jest to szczególnie ważne na przykład podczas sprawdzania poprawności krzyżowej ).

Biorąc pod uwagę ograniczenie , wektor jest $u'u=1$ $u$

u = \frac{X^{'} y}{‖ X^{'} y ‖} .

$u=\frac{X'y}{\|X'y\|}.$

Za pomocą małej symulacji można uzyskać w następujący sposób:

set.seed(101)
X <- replicate(2, rnorm(100))
y <- 0.6*X[,1] + 0.7*X[,2] + rnorm(100)
X <- apply(X, 2, scale)
y <- scale(y)

# NIPALS (PLS1)
u <- crossprod(X, y)
u <- u/drop(sqrt(crossprod(u)))         # X weights
t  <- X%*%u
p <- crossprod(X, t)/drop(crossprod(t)) # X loadings

Możesz porównać powyższe wyniki ( u=[0.5792043;0.8151824]w szczególności) z tym, co dadzą pakiety R. Np. Używając NIPALS z pakietu chemometrii (inna implementacja, o której wiem, że jest dostępna w pakiecie mixOmics ), otrzymalibyśmy:

library(chemometrics)
pls1_nipals(X, y, 1)$W  # X weights [0.5792043;0.8151824]
pls1_nipals(X, y, 1)$P  # X loadings

Podobne wyniki można uzyskać za pomocą plsrdomyślnego algorytmu PLS jądra:

> library(pls)
> as.numeric(loading.weights(plsr(y ~ X, ncomp=1)))
[1] 0.5792043 0.8151824

We wszystkich przypadkach możemy sprawdzić, czy ma długość 1. $u$

Pod warunkiem, że zmienisz funkcję, aby zoptymalizować na czytającą

f <- function(u) cov(y, X%*%(u/sqrt(crossprod(u))))

a unastępnie normalizować ( u <- u/sqrt(crossprod(u))), powinieneś być bliżej powyższego rozwiązania.

Sidenote : Jako kryterium (1) jest równa znajduje się w lewym pojedynczej wektora z SVD odpowiadający największej wartości własnej:

max u^{'} X^{'} Y v,

$\max u'X'Yv,$

u

$u$

X^{'} Y

$X'Y$

svd(crossprod(X, y))$u

W bardziej ogólnym przypadku (PLS2) sposobem na podsumowanie powyższego jest stwierdzenie, że pierwsze wektory kanoniczne PLS są najlepszym przybliżeniem macierzy kowariancji X i Y w obu kierunkach.

Bibliografia

Tenenhaus, M (1999). L'approche PLS . Revue de Statistique Appliquée , 47 (2), 5-40.
ter Braak, CJF i de Jong, S (1993). Funkcja celu częściowej regresji metodą najmniejszych kwadratów . Journal of Chemometrics , 12, 41–54.
Abdi, H (2010). Częściowa regresja najmniejszych kwadratów i rzut na ukrytą regresję struktury (regresja PLS) . Wiley Interdisciplinary Reviews: Statystyka obliczeniowa , 2, 97-106.
Boulesteix, AL i Strimmer, K (2007). Częściowe najmniejsze kwadraty: wszechstronne narzędzie do analizy wielowymiarowych danych genomowych . Briefings in Bioinformatics , 8 (1), 32-44.

— chl
źródło

Dzięki chl. Przeczytam twoją odpowiedź, kiedy tylko będzie to możliwe (i na pewno głosuję i kliknę znacznik wyboru!)

— Stéphane Laurent,

Właśnie przeczytałem twoją odpowiedź - gratuluję i bardzo dziękuję.

— Stéphane Laurent,