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

12

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 iab

> ( 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 ? bab

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

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 .X ˜ Y Y ˜ X X t 1 t 1 = a ˜ X 1 + b ˜ X 2YXY~YX~Xt1t1=aX~1+bX~2b t 1 , ~ YT 1 Tabt1,Y~Zatem jest to równoważne z maksymalizacją korelacji między i , prawda?t1Y

Stéphane Laurent
źródło
2
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
1
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? X ~ Y T ~ X X T 1 T 1 = ~ X 1 + b ~ X 2 b t 1 , ~ YYXY~YX~Xt1t1=aX~1+bX~2abt1,Y~
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ś).a2+b21?coef.mvr
chl

Odpowiedzi:

17

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 uv

maxcov(Xu,Yv).(1)
Y
cov(Xu,y)Var(Xu)1/2×cor(Xu,y)×Var(y)1/2,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ówVar(y)uVar(Xu)1/2×cor(Xu,y)X=[x_1;x_2]x1x2Var(x1)=Var(x2)=1Var(Xu)1u.

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 chemometriiuwinieta 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 jestuu=1u

u=XyXy.

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:

maxuXYv,
uXY
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

  1. Tenenhaus, M (1999). L'approche PLS . Revue de Statistique Appliquée , 47 (2), 5-40.
  2. ter Braak, CJF i de Jong, S (1993). Funkcja celu częściowej regresji metodą najmniejszych kwadratów . Journal of Chemometrics , 12, 41–54.
  3. 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.
  4. 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,