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

29

J'utilise toujours lm()en R pour effectuer une régression linéaire de sur . Cette fonction renvoie un coefficient tel queyxβ

y=β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 .βxy

Dail
la source
Un instant les gars, je suis un peu confus. regardez: zoonek2.free.fr/UNIX/48_R/09.html Cela s'appelle PCA (Analyse en Composantes Principales, alias "régression orthogonale" ou "sommes perpendiculaires de carrés" ou "total des moindres carrés") donc je pense que nous parlons à propos de TLS avec princomp () Non?
Dail
Non; ce sont deux choses différentes, voir l'article de wikipedia sur PCA. Le fait qu'il soit utilisé ici est un hack (je ne sais pas exactement, mais je vais le vérifier); c'est pourquoi l'extraction complexe des coefficients.
1
Une question connexe: stats.stackexchange.com/questions/2691/… et un article de blog est référencé par l'une des réponses: cerebralmastication.com/2010/09/…
Jonathan

Réponses:

48

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:xxy

OLS vs TLS

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=βxyy^(x,y)β(x,y)2×2Σ ; alors β = v y / v x .v=(vx,vy)β=vy/vx

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 .xyβ0=y¯β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 + 1yxi

y=β1x1++βpxp.
yy^(x,y)Rp+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 β .(x,y)p+1VppHp+1vp+1Hpβ

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 + 1xi=0ikxk=1y^=βk

(0,,1,,βk)H
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 .
vp+1=(v1,,vp+1)H
vk+βkvp+1=0βk=-vk/vp+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 β .Xyβ0=y¯-X¯β

X(X,y)vy(1)/vX(1)=-vX(2)/vy(2)

Solution de formulaire fermé pour TLS

β

Xyvp+1[Xy]σp+12-vp+1/vp+1=(β-1)

(XXXyyXyy)(β-1)=σp+12(β-1),
βTLS=(XX-σp+12je)-1Xy,
βOLS=(XX)-1Xy.

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.

amibe dit réintégrer Monica
la source
1
Je ne connais pas R, mais je voulais quand même fournir des extraits R pour référence future. Il y a beaucoup de gens compétents ici en R. N'hésitez pas à modifier mes extraits si nécessaire! Merci.
Amoeba dit Reinstate Monica
(0,,1,,βk)
XjeXk=1y=βjXjy=βk1=βk(0,,1,βk)y=βjXj
Amoeba dit Reinstate Monica
Il me semble avoir mal lu cette partie mais maintenant c'est clair. Merci aussi pour la clarification.
JohnK
2
Dans R, vous préférerez peut-être " vecteurs propres (cov (cbind (x, y))) $ " à "prcomp (cbind (x, y)) $ rotation" car le premier est beaucoup plus rapide pour les vecteurs plus grands.
Thomas Browne
9

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)
}
cashoes
la source
4

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.


la source
Je pensais que Deming dans le package MethComp était TLS - quelle est la différence?
mark999
Vous devez lui donner le ratio d'erreurs sur x et y; pure TLS optimise cela.