Partielle Regression der kleinsten Quadrate in R: Warum entspricht PLS für standardisierte Daten nicht der Maximierung der Korrelation?

12

Ich bin sehr neu in Partial Least Squares (PLS) und versuche, die Ausgabe der R-Funktion plsr()im plsPaket zu verstehen . Lassen Sie uns Daten simulieren und den PLS ausführen:

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)

Ich hatte erwartet, dass die folgenden Nummern a und 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

werden berechnet, um zu maximieren

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

aber das ist nicht genau der Fall:

> 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

Handelt es sich um einen numerischen Fehler, oder verstehe ich die Natur von und falsch?ab ?

Ich würde auch gerne wissen, was diese Koeffizienten sind:

> p$coef
, , 1 comps

           yy
xx1 0.6672848
xx2 0.6368604 

EDIT : Jetzt sehe ich was p$coefist:

> 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 

Ich denke also, ich habe Recht mit der Natur von undab .

EDIT: In Anbetracht der Kommentare von @chl denke ich, dass meine Frage nicht klar genug ist, also lass mich mehr Details angeben. In meinem Beispiel gibt es einen Vektor von Antworten und eine zweispaltige Matrix von Prädiktoren, und ich verwende die normalisierte Version von und die normalisierte Version von (zentriert und dividiert durch Standardabweichungen). Die Definition der ersten PLS-Komponente ist wobei und gewählt sind, dass ein Maximalwert des inneren Produkts erhalten wird. Daher entspricht dies der Maximierung der Korrelation zwischenX ~ Y Y ~ X X t 1 t 1 = a ~ X 1 + b ~ X 2 a b t 1 , ~ Yt 1 YYXY~YX~Xt1t1=aX~1+bX~2abt1,Y~ .t1 und , nicht wahr?Y

Stéphane Laurent
quelle
2
Die PLS-Regression maximiert die Faktorwerte (die als Produkt von Rohdaten mit Kovarianz des Ladevektors (der Ladevektoren) berechnet werden , nicht die Korrelation (wie in der kanonischen Korrelationsanalyse). plsIn diesem JSS-Dokument finden Sie einen guten Überblick über das Paket und die PLS-Regression .
Chl
1
Da alle Vektoren zentriert und normalisiert sind, ist Kovarianz eine Korrelation, nicht wahr? Tut mir leid, aber das JSS-Papier ist für Anfänger zu technisch.
Stéphane Laurent
Im Allgemeinen gibt es einen asymmetrischen Deflationsprozess (resultierend aus der Regression der linearen Kombination eines Blocks auf die lineare Kombination des anderen), der die Dinge etwas kompliziert. Ich lieferte in dieser Antwort ein schematisches Bild . Hervé Abdi gab einen allgemeinen Überblick über die PLS-Regression, und auch Wegelins Umfrage zu PLS-Methoden (Partial Least Squares) ist sehr nützlich. An dieser Stelle sollte ich wahrscheinlich alle diese Kommentare in eine Antwort umwandeln ...
chl
In meinem Beispiel gibt es einen Vektor von Antworten und eine zweispaltige Matrix X von Prädiktoren, und ich verwende die normalisierte Version ˜ Y von Y und die normalisierte Version ˜ X von X (zentriert und dividiert durch Standardabweichungen). Meine Definition der ersten PLS-Komponente t 1 ist t 1 = a X 1 + b X 2, wobei a und b so gewählt sind, dass ein Maximalwert des Skalarprodukts t 1 erhalten wird.YXY~YX~Xt1t1=aX~1+bX~2ab. Ist es nicht die gute Definition? t1,Y~
Stéphane Laurent
Entschuldigung, @Stéphane, meine obigen Kommentare haben nicht berücksichtigt, dass Sie nur eine Komponente angefordert haben (daher spielt die Deflation hier keine entscheidende Rolle). Es scheint jedoch, dass Ihre Optimierungsfunktion keine Einheitsnorm-Gewichtungsvektoren auferlegt, so dass am Ende . (Übrigens werden Sie mehr Informationen darüber erhalten, was diese "Koeffizienten" sind, aber Sie haben das anscheinend bereits selbst entdeckt.)a2+b21?coef.mvr
chl

Antworten:

17

Die PLS-Regression beruht auf iterativen Algorithmen (z. B. NIPALS, SIMPLS). Ihre Beschreibung der Hauptideen ist richtig: Wir suchen einen (PLS1, eine Antwortvariable / mehrere Prädiktoren) oder zwei (PLS2, mit verschiedenen Modi, mehreren Antwortvariablen / mehreren Prädiktoren) Vektoren von Gewichten, (und v ) Zum Beispiel, um eine lineare Kombination der ursprünglichen Variablen zu bilden, so dass die Kovarianz zwischen Xu und Y (Yv, für PLS2) maximal ist. Konzentrieren wir uns darauf, das erste Gewichtungspaar zu extrahieren, das der ersten Komponente zugeordnet ist. Formal lautet das Kriterium zur Optimierung max cov ( X u , Y v )uv In Ihrem Fall Y ist univariate, so dass es zu maximierenbeträgt COV ( X u , y ) Var ( X u ) 1 / 2 × cor ( X u , y ) × Var ( y ) 1 / 2 ,

maxcov(Xu,Yv).(1)
Y Da Var ( y ) hängt nicht von u , müssen wir maximieren Var ( X u ) 1 / 2 × cor ( X u , y ) . Überlegen wir uns, wo die Daten individuell standardisiert sind (ich habe anfangs den Fehler gemacht, Ihre Linearkombination anstelle von x 1 und x 2 separat zuskalieren!), So dass Var ( x 1 )
cov(Xu,y)Var(Xu)1/2×cor(Xu,y)×Var(y)1/2,st.u=1.
Var(y)uVar(Xu)1/2×cor(Xu,y)X=[x_1;x_2]x1x2 ; Jedoch Var ( X u ) 1 undabhängig von u . Zusammenfassend lässt sich sagen, dass dieMaximierung der Korrelation zwischen der latenten Komponente und der Antwortvariablen nicht zu den gleichen Ergebnissen führtVar(x1)=Var(x2)=1Var(Xu)1u .

Ich sollte Arthur Tenenhaus danken der mich in die richtige Richtung wies.

pls. regressionpls.pcru , aber der ChemometrieVignette bietet auch eine gute Diskussion (S. 26-29). Von besonderer Bedeutung ist auch die Tatsache, dass die meisten PLS-Routinen (zumindest die, die ich in R kenne) davon ausgehen, dass Sie nicht standardisierte Variablen bereitstellen, da die Zentrierung und / oder Skalierung intern erfolgt (dies ist besonders wichtig, wenn Sie beispielsweise eine Kreuzvalidierung durchführen ).

uu=1u

u=XyXy.

Mit einer kleinen Simulation kann es wie folgt erhalten werden:

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

Sie können die obigen Ergebnisse ( u=[0.5792043;0.8151824]insbesondere) mit den Ergebnissen von R-Paketen vergleichen. Wenn Sie beispielsweise NIPALS aus dem Chemometrics- Paket verwenden (eine andere mir bekannte Implementierung ist im mixOmics- Paket verfügbar ), erhalten Sie Folgendes :

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

Ähnliche Ergebnisse würden mit plsrund seinem Standard-Kernel-PLS-Algorithmus erzielt :

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

u Länge 1 hat.

Vorausgesetzt, Sie ändern Ihre Funktion, um sie auf eine lesbare zu optimieren

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

und udanach normalisieren (u <- u/sqrt(crossprod(u)) ), du solltest näher an der obigen Lösung sein.

maxuXYv,
uXY
svd(crossprod(X, y))$u

Im allgemeineren Fall (PLS2) kann man zusammenfassend sagen, dass die ersten kanonischen PLS-Vektoren die beste Approximation der Kovarianzmatrix von X und Y in beide Richtungen sind.

Verweise

  1. Tenenhaus, M (1999). L'approche PLS . Revue de Statistique Appliquée , 47 (2), 5-40.
  2. ter Braak, CJF und de Jong, S. (1993). Die objektive Funktion der partiellen Regression der kleinsten Quadrate . Journal of Chemometrics , 12, 41–54.
  3. Abdi, H (2010). Partielle Regression der kleinsten Quadrate und Projektion auf die latente Strukturregression (PLS-Regression) . Wiley Interdisciplinary Reviews: Computational Statistics , 2, 97-106.
  4. Boulesteix, AL und Strimmer, K (2007). Partial Least Squares: ein vielseitiges Tool zur Analyse hochdimensionaler Genomdaten . Briefings in Bioinformatics , 8 (1), 32-44.
chl
quelle
Danke chl. Ich werde Ihre Antwort nach Möglichkeit lesen (und mit Sicherheit zustimmen und auf das Häkchen klicken!)
Stéphane Laurent
Ich habe gerade Ihre Antwort gelesen - herzlichen Glückwunsch und vielen Dank.
Stéphane Laurent