Régression des moindres carrés partiels dans R: pourquoi le PLS sur des données standardisées n'est-il pas équivalent à maximiser la corrélation?

Je suis très nouveau dans les moindres carrés partiels (PLS) et j'essaie de comprendre la sortie de la fonction R plsr()dans le plspackage. Simulons les données et exécutons le 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)

Je m'attendais à ce que les nombres suivants $a$ et $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

sont calculés afin de maximiser

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

mais ce n'est pas exactement le cas:

> 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

Est-ce une erreur numérique ou est-ce que je comprends mal la nature de et ? $a$ $b$

J'aimerais aussi savoir quels sont ces coefficients:

> p$coef
, , 1 comps

           yy
xx1 0.6672848
xx2 0.6368604

EDIT : Maintenant, je vois ce qui p$coefest:

> 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

Je pense donc que j'ai raison sur la nature de et . $a$ $b$

EDIT: Au vu des commentaires de @chl, je pense que ma question n'est pas assez claire, alors laissez-moi vous donner plus de détails. Dans mon exemple, il y a un vecteur de réponses et une matrice deux colonnes de prédicteurs et j'utilise la version normalisée de et la version normalisée de (centrée et divisée par les écarts-types). La définition du premier composant PLS est avec et choisis pour avoir une valeur maximale du produit intérieur . $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$ Par conséquent, cela revient à maximiser la corrélation entre et , n'est-ce pas? $t_1$ $Y$

r regression partial-least-squares Stéphane Laurent
la source

La régression PLS maximise les scores des facteurs (qui sont calculées en tant que produit de données brutes avec le vecteur de charges (s)) covariance , pas de corrélation (comme cela se fait dans l' analyse de corrélation canonique). Il y a un bon aperçu du plspackage et de la régression PLS dans ce document JSS .

chl

Puisque tous les vecteurs sont centrés et normalisés, la covariance est une corrélation, n'est-ce pas? Désolé mais le papier JSS est trop technique pour un débutant.

Stéphane Laurent

En général, il y a un processus de déflation asymétrique (résultant de la régression de la combinaison linéaire d'un bloc sur la combinaison linéaire de l'autre) qui complique un peu les choses. J'ai fourni une image schématique dans cette réponse . Hervé Abdi a donné un aperçu général de la régression PLS, et les méthodes d' enquête de Wegelin sur les moindres carrés partiels (PLS) sont également très utiles. À ce stade, je devrais probablement convertir tous ces commentaires en une réponse ...

chl

Dans mon exemple, il y a un vecteur

Y

$Y$ de réponses et une matrice

deux colonnes

X

$X$ de prédicteurs et j'utilise la version normalisée

\tilde{Y}

$\tilde Y$ de

Y

$Y$ et la version normalisée

\tilde{X}

$\tilde X$ de

X

$X$ (centrée et divisée par les écarts-types). Ma définition du premier composant PLS

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$

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

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

$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.$ Puisque ne dépend pas de , nous devons maximiser . Considérons , où les données sont standardisées individuellement (j'ai fait l'erreur de mettre à l'échelle votre combinaison linéaire au lieu de et séparément!), De sorte que ; cependant, et dépend de . En conclusion, maximiser la corrélation entre la composante latente et la variable de réponse ne donnera pas les mêmes résultats

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

Je dois remercier Arthur Tenenhaus qui m'a pointé dans la bonne direction.

L'utilisation de vecteurs de poids unitaire n'est pas restrictive et certains packages ( pls. regressiondans plsgenomics , basés sur le code du package précédent de Wehrens) renverront des pls.pcrvecteurs de poids non standardisés (mais avec des composants latents toujours de la norme 1), si demandé. Mais la plupart des packages PLS renverront normalisé , y compris celui que vous avez utilisé, notamment ceux implémentant l'algorithme SIMPLS ou NIPALS; J'ai trouvé un bon aperçu des deux approches dans la présentation de Barry M. Wise, Propriétés de la régression des moindres carrés partiels (PLS), et des différences entre les algorithmes , mais la chimiométrie $u$ la vignette offre également une bonne discussion (pp. 26-29). Le fait que la plupart des routines PLS (au moins celle que je connais dans R) supposent que vous fournissez des variables non standardisées car le centrage et / ou la mise à l'échelle sont gérés en interne est particulièrement important (cela est particulièrement important lors de la validation croisée, par exemple). ).

Étant donné la contrainte , le vecteur est $u'u=1$ $u$

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

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

En utilisant une petite simulation, elle peut être obtenue comme suit:

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

Vous pouvez comparer les résultats ci-dessus ( u=[0.5792043;0.8151824]en particulier) avec ce que les packages R donneraient. Par exemple, en utilisant NIPALS du chimiométrie package ( une autre mise en œuvre que je sais est disponible dans le mixOmics package), on obtiendrait:

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

Des résultats similaires seraient obtenus avec plsrson algorithme PLS de noyau par défaut:

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

Dans tous les cas, on peut vérifier que est de longueur 1. $u$

À condition de changer votre fonction pour optimiser celle qui lit

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

et normalisez uensuite ( u <- u/sqrt(crossprod(u))), vous devriez être plus proche de la solution ci-dessus.

Sidenote : Comme le critère (1) est équivalent à peut être trouvé comme le vecteur singulier gauche de la SVD de correspondant à la plus grande valeur propre:

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

$\max u'X'Yv,$

u

$u$

X^{'} Y

$X'Y$

svd(crossprod(X, y))$u

Dans le cas plus général (PLS2), une façon de résumer ce qui précède est de dire que les premiers vecteurs canoniques PLS sont la meilleure approximation de la matrice de covariance de X et Y dans les deux directions.

Les références

Tenenhaus, M (1999). L'approche PLS . Revue de Statistique Appliquée , 47 (2), 5-40.
ter Braak, CJF et de Jong, S (1993). La fonction objective de la régression des moindres carrés partiels . Journal of Chemometrics , 12, 41–54.
Abdi, H (2010). Régression des moindres carrés partiels et projection sur la régression de la structure latente (régression PLS) . Wiley Interdisciplinary Reviews: Computational Statistics , 2, 97-106.
Boulesteix, AL et Strimmer, K (2007). Moindres carrés partiels: un outil polyvalent pour l'analyse de données génomiques de haute dimension . Briefings in Bioinformatics , 8 (1), 32-44.

chl
la source

Merci chl. Je vais lire votre réponse chaque fois que possible (et sûrement voter et cliquer sur la coche!)

Stéphane Laurent

Je viens de lire votre réponse - félicitations et merci beaucoup.

Stéphane Laurent

Régression des moindres carrés partiels dans R: pourquoi le PLS sur des données standardisées n'est-il pas équivalent à maximiser la corrélation?

Réponses:

Les références