Comment effectuer une régression orthogonale (total des moindres carrés) via l'ACP?

J'utilise toujours lm()en R pour effectuer une régression linéaire de sur . Cette fonction renvoie un coefficient tel que$y$$x$$\beta$

$$y = \beta x.$$

Aujourd'hui, j'ai appris le nombre total de moindres carrés et cette princomp()fonction (analyse en composantes principales, ACP) peut être utilisée pour l'exécuter. Cela devrait être bon pour moi (plus précis). J'ai fait quelques tests en utilisant princomp(), comme:

r <- princomp( ~ x + y)

Mon problème est: comment interpréter ses résultats? Comment obtenir le coefficient de régression? Par "coefficient", je veux dire le nombre que je dois utiliser pour multiplier la valeur pour donner un nombre proche de .$\beta$$x$$y$

Moindres carrés ordinaires vs moindres carrés totaux

Considérons d'abord le cas le plus simple d'une seule variable prédictive (indépendante) . Par souci de simplicité, laissez x et y être centrés, c'est-à-dire que l'interception est toujours nulle. La différence entre la régression OLS standard et la régression TLS "orthogonale" est clairement montrée sur cette figure (adaptée par moi) de la réponse la plus populaire dans le fil le plus populaire sur PCA:$x$$x$$y$

OLS ajuste l'équation en réduisant au minimum les distances au carré entre les valeurs observées y et les valeurs prédites y . TLS correspond à la même équation en minimisant les distances au carré entre les points ( x , y ) et leur projection sur la ligne. Dans ce cas le plus simple, la ligne TLS est simplement le premier composant principal des données 2D. Pour trouver β , faites l'ACP sur les points ( x , y ) , c'est-à-dire construisez la matrice de covariance 2 × 2 Σ et trouvez son premier vecteur propre v$y=\beta x$$y$$\hat y$$(x,y)$$\beta$$(x,y)$$2\times 2$$\boldsymbol \Sigma$ ; alors β = v y / v x .$\mathbf v = (v_x, v_y)$$\beta = v_y/v_x$

En Matlab:

 v = pca([x y]);    //# x and y are centered column vectors
 beta = v(2,1)/v(1,1);

Dans R:

 v <- prcomp(cbind(x,y))$rotation
 beta <- v[2,1]/v[1,1]

Soit dit en passant, cela produira une pente correcte même si et y n'étaient pas centrés (car les fonctions PCA intégrées effectuent automatiquement le centrage). Pour récupérer l'ordonnée à l'origine, calculez β 0 = ˉ y - β ˉ x .$x$$y$$\beta_0 = \bar y - \beta \bar x$

OLS vs TLS, régression multiple

Étant donné une variable dépendante et de nombreuses variables indépendantes x i (encore une fois, toutes centrées pour plus de simplicité), la régression correspond à une équation y = β 1 x 1 + … + β p x p . OLS est la place en minimisant les erreurs quadratiques entre les valeurs observées de y et les valeurs prédites y . TLS fait l'ajustement en minimisant les distances au carré entre observées ( x , y ) ∈ R p + $y$$x_i$

$$y= \beta_1 x_1 + \ldots + \beta_p x_p.$$ $y$$\hat y$$(\mathbf x, y)\in\mathbb R^{p+1}$ points et les points les plus proches sur le plan de régression / hyperplan.

Notez qu'il n'y a plus de "ligne de régression"! L'équation ci-dessus spécifie un hyperplan : c'est un plan 2D s'il y a deux prédicteurs, un hyperplan 3D s'il y a trois prédicteurs, etc. Donc la solution ci-dessus ne fonctionne pas: on ne peut pas obtenir la solution TLS en prenant uniquement le premier PC (qui est une ligne). Pourtant, la solution peut être facilement obtenue via PCA.

Comme précédemment, l'ACP est effectuée sur points. On obtient p + 1 vecteurs propres dans les colonnes de V . Les p premiers vecteurs propres définissent un hyperplan H à p dimensions dont nous avons besoin; le dernier (nombre p + 1 ) vecteur propre v p + 1 lui est orthogonal. La question est de savoir comment transformer la base de H donnée par les p premiers vecteurs propres en coefficients β .$(\mathbf x, y)$$p+1$$\mathbf V$$p$$p$$\mathcal H$$p+1$$\mathbf v_{p+1}$$\mathcal H$$p$$\boldsymbol \beta$

Observe que si l' on pose pour tout i ≠ k et que x k = 1 , alors Y = ß k , soit le vecteur ( 0 , ... , 1 , ... , β k ) ∈ H se trouve dans l'hyperplan H . Par contre, on sait que v p + 1 = ( v 1 , … , v p + $x_i=0$$i \ne k$$x_k=1$$\hat y=\beta_k$

$$(0,\ldots, 1, \ldots, \beta_k) \in \mathcal H$$ $\mathcal H$ est orthogonal. C'est-à-dire que leur produit scalaire doit être nul: 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}.$$

En Matlab:

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

Dans R:

 v <- prcomp(cbind(X,y))$rotation
 beta <- -v[-ncol(v),ncol(v)] / v[ncol(v),ncol(v)]

Encore une fois, cela produira des pentes correctes même si et y n'étaient pas centrés (car les fonctions PCA intégrées effectuent automatiquement le centrage). Pour récupérer l'ordonnée à l'origine, calculez β 0 = ˉ y - ˉ x β .$x$$y$$\beta_0 = \bar y - \bar {\mathbf x} \boldsymbol \beta$

$x$$(x,y)$$v^{(1)}_y/v^{(1)}_x=-v^{(2)}_x/v^{(2)}_y$

Solution de formulaire fermé pour TLS

$\boldsymbol \beta$

$\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),$$ $$\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.$$

Régression multiple multivariée

La même formule peut être généralisée au cas multivarié, mais même pour définir ce que fait le TLS multivarié, il faudrait une certaine algèbre. Voir Wikipedia sur TLS . La régression OLS multivariée équivaut à un ensemble de régressions OLS univariées pour chaque variable dépendante, mais dans le cas TLS, ce n'est pas le cas.

Basé sur l'implémentation naïve de GNU Octave trouvée ici , quelque chose comme ça pourrait (grain de sel, il est tard) fonctionner.

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

princompexécute l' analyse de la composante principale au lieu de la régression des moindres carrés totaux. Pour autant que je sache, il n'y a pas de fonction R ni de package qui fait TLS; il y a tout au plus une régression de Deming dans MethComp .
Pourtant, veuillez considérer cela comme une suggestion que cela ne vaut probablement pas la peine.