Régression polynomiale brute ou orthogonale?

22

Je veux régresser une variable sur x , x 2 , , x 5 . Dois-je le faire en utilisant des polynômes bruts ou orthogonaux? J'ai regardé la question sur le site qui traite de ces derniers, mais je ne comprends pas vraiment quelle est la différence entre les utiliser. yX,X2,,X5

Pourquoi ne puis-je pas simplement faire une régression "normale" pour obtenir les coefficients de y = 5 i = 0 β i x iβjey=je=05βjeXje (avec les valeurs de p et tous les autres trucs sympas) et au lieu de cela, je dois m'inquiéter si j'utilise polynômes bruts ou orthogonaux? Ce choix me semble sortir du cadre de ce que je veux faire.

Dans le livre de statistiques que je lis actuellement (ISLR par Tibshirani et al), ces choses n'étaient pas mentionnées. En fait, ils ont été minimisés d'une certaine manière.
La raison en est, AFAIK, que dans la lm()fonction dans R, utiliser des y ~ poly(x, 2)quantités à utiliser des polynômes orthogonaux et utiliser des y ~ x + I(x^2)quantités à utiliser des polynômes . Mais à la page 116, les auteurs disent que nous utilisons la première option parce que cette dernière est "lourde", ce qui ne laisse aucune indication que ces commandes ont en fait des choses complètement différentes (et ont des sorties différentes en conséquence).
(troisième question) Pourquoi les auteurs de l'ISLR confondraient-ils leurs lecteurs comme ça?

l7ll7
la source
1
@Sycorax Je sais que cela polya quelque chose à voir avec les polynômes orthogonaux et je (x ^ 2) ne le fait pas (bien que je ne connaisse pas les détails) - mais quand même, pourquoi les auteurs de l'ISLR recommanderaient-ils alors une méthode qui ne fonctionne pas ? Cela semble très trompeur si les deux commandes semblent faire la même chose, mais une seule est en fait correcte.
l7ll7
1
@gung J'ai regardé la documentation de polyce problème et y ai déjà passé un certain temps, mais je n'arrive pas à comprendre pourquoi poly (x, 2) et x + I (x ^ 2) font une différence? Pourriez-vous m'éclairer ici dans les commentaires, si la question est hors sujet?
l7ll7
1
@gung J'ai fait une réédition complète de ma question. Ce choix brut / orthogonal me déroute encore plus - auparavant, je pensais que c'était juste une Rtechnicité mineure , que je ne comprenais pas, mais maintenant il semble que ce soit un problème de statistiques à part entière qui m'empêche de coder une régression qui ne devrait pas être difficile à coder.
l7ll7
2
@gung Cela m'a en fait plus troublé que cela n'a aidé. Auparavant, je pensais que je devais simplement utiliser des polynômes orthogonaux, car cela semblait être la bonne façon, mais dans cette réponse, des polynômes bruts sont utilisés. Étonnamment, tout le monde sur le net crie "RTFM", mais il n'y a en fait pas de réponse claire, quand utiliser quoi. (Votre lien ne donne pas non plus de réponse à cela, juste un exemple, lorsque la pol.
Orth
2
À moins que vous ne travailliez dans un domaine physique ou d'ingénierie qui indique que la réponse sera un polynôme quintique, la bonne approche est presque certainement de ne pas faire de régression polynomiale en premier lieu. Investissez vos degrés de liberté dans une spline ou quelque chose qui serait beaucoup plus flexible et stable que l'ajustement polynomial.
whuber

Réponses:

10

Je crois que la réponse est moins sur la stabilité numérique (bien que cela joue un rôle) et plus sur la réduction de la corrélation.

Essentiellement - le problème se résume au fait que lorsque nous régressons contre un groupe de polynômes d'ordre élevé, les covariables contre lesquelles nous régressons deviennent fortement corrélées. Exemple de code ci-dessous:

x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
cor(raw.poly)
cor(orthogonal.poly)

C'est extrêmement important. À mesure que les covariables deviennent plus corrélées, notre capacité à déterminer lesquelles sont importantes (et quelle est la taille de leurs effets) s'érode rapidement. Ceci est généralement appelé le problème de la multicolinéarité. À la limite, si nous avions deux variables qui étaient entièrement corrélées, lorsque nous les régressons par rapport à quelque chose, il est impossible de faire la distinction entre les deux - vous pouvez considérer cela comme une version extrême du problème, mais ce problème affecte nos estimations pour moindre degré de corrélation également. Ainsi, dans un vrai sens - même si l'instabilité numérique n'était pas un problème - la corrélation des polynômes d'ordre supérieur endommage énormément nos routines d'inférence. Cela se manifestera par des erreurs standard plus importantes (et donc des statistiques t plus petites) que vous verriez autrement (voir l'exemple de régression ci-dessous).

y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
summary(raw.mod)
summary(orthogonal.mod)

Si vous exécutez ce code, l'interprétation est un peu difficile car les coefficients changent tous et les choses sont donc difficiles à comparer. En regardant les statistiques T, nous pouvons voir que la capacité à déterminer les coefficients était BEAUCOUP plus grande avec les polynômes orthogonaux. Pour les 3 coefficients pertinents, j'ai obtenu des statistiques t de (560,21,449) pour le modèle orthogonal, et seulement (28, -38,121) pour le modèle polynomial brut. C'est une énorme différence pour un modèle simple avec seulement quelques termes polynomiaux d'ordre relativement faible qui importaient.

Cela ne veut pas dire que cela vient sans frais. Il y a deux coûts principaux à garder à l'esprit. 1) nous perdons une certaine interprétabilité avec les polynômes orthogonaux. Nous pourrions comprendre ce que x**3signifie le coefficient sur , mais interpréter le coefficient sur x**3-3x(le troisième poly hermite - pas nécessairement ce que vous utiliserez) peut être beaucoup plus difficile. Deuxièmement - lorsque nous disons que ces polynômes sont orthogonaux - nous voulons dire qu'ils sont orthogonaux par rapport à une certaine mesure de distance. Choisir une mesure de distance adaptée à votre situation peut être difficile. Cependant, cela dit, je crois que la polyfonction est conçue pour choisir de telle sorte qu'elle soit orthogonale par rapport à la covariance - ce qui est utile pour les régressions linéaires.

user5957401
la source
3
-1. Les erreurs standard les plus importantes que vous voyez sur les coefficients d'ordre inférieur sont un hareng rouge. Les coefficients d'ordre inférieur dans vos deux modèles estiment des choses complètement différentes, donc comparer leurs erreurs standard n'a aucun sens. Le coefficient d'ordre le plus élevé est le seul à estimer la même chose dans les deux modèles, et vous verrez que la statistique t est identique, que les polynômes soient orthogonaux ou non. Vos deux modèles sont statistiquement équivalents en termes de valeurs ajustées, R ^ 2, etc., ils diffèrent principalement uniquement par l'interprétation des coefficients
Jake Westfall
@JakeWestfall, je ne pense pas être d'accord avec vous. Tout d'abord, l'exécution du code produit des valeurs qui sont différentes pour tous les ordres polynomiaux, pas tous sauf un - en substance, il prend le polynôme et fait PCA dessus. Deuxièmement, et plus important encore, les statistiques t sont sensiblement différentes - l'exécution du code dans ma réponse confirmera que - fonctionnellement, nous résolvons le problème de multicolinéarité. Vous avez raison de dire que les valeurs ajustées, R ^ 2, les tests F, etc. ne changent pas. En fait, c'est une raison d'orthogonaliser - cela ne change rien, sauf notre capacité à détecter les termes polynomiaux .
user5957401
1
Re: le premier point, désolé, je voulais parler du t-stat du terme d'ordre le plus élevé, pas de son coefficient. Ce prédicteur est mis à l'échelle + décalé entre les modèles, donc oui, le coef change, mais il teste le même effet de fond, comme le montre t
Jake Westfall
Concernant le deuxième point, la raison pour laquelle "les statistiques t sont sensiblement différentes" pour les termes d'ordre inférieur est, encore une fois, parce qu'elles estiment des choses complètement différentes dans les deux modèles. Considérons l'effet linéaire: raw.modil estime la pente de la courbe à x = 0, orthogonal.modil estime la pente marginale (c'est-à-dire identique à l' lm(y ~ poly(x,1))endroit où les termes d'ordre supérieur sont omis). Il n'y a aucune raison pour que les estimations de ces estimations complètement différentes aient des erreurs-types comparables. On peut facilement construire un contre-exemple où les raw.modstatistiques t sont beaucoup plus élevées
Jake Westfall
@JakeWestfall. Je pense toujours que vous manquez la multicolinéarité. Cependant, nous semblons parler les uns des autres, et il y a peut-être une solution. Vous dites que vous pouvez facilement construire un contre-exemple, veuillez le faire. Je pense que voir le dgp que vous avez en tête clarifierait beaucoup pour moi. Pour le moment, les seules choses que j'ai pu trouver qui pourraient se comporter comme vous le décrivez impliquent une mauvaise spécification du modèle.
user5957401
8

Pourquoi ne puis-je pas simplement faire une régression "normale" pour obtenir les coefficients?

0.40.4000000059604644775390625

L'utilisation d'un polynôme brut causera un problème car nous en aurons un nombre énorme. Voici une petite preuve: nous comparons le nombre de conditions de matrice avec un polynôme brut et orthogonal.

> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13

Vous pouvez également consulter ma réponse ici pour un exemple.

Pourquoi existe-t-il des coefficients importants pour les polynômes d'ordre supérieur

Haitao Du
la source
6
Vous semblez utiliser des flotteurs à simple précision et les citer à une précision quadruple! Comment est-ce arrivé? À l'exception des GPU, presque tous les calculs statistiques utilisent au moins une double précision. Par exemple, dans Rla sortie de print(0.4, digits=20)is 0.40000000000000002.
whuber
6

J'ai l'impression que plusieurs de ces réponses ratent complètement le propos. La réponse d'Haitao aborde les problèmes de calcul avec l'ajustement des polynômes bruts, mais il est clair que OP pose des questions sur les statistiques différences entre les deux approches. Autrement dit, si nous avions un ordinateur parfait qui pourrait représenter toutes les valeurs exactement, pourquoi préférerions-nous une approche à l'autre?

R2XOuiX=0X=0X

data("iris")

#Raw:
fit.raw <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
summary(fit.raw)

#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        1.1034     0.1304   8.464 2.50e-14 ***
#> Petal.Width        1.1527     0.5836   1.975  0.05013 .  
#> I(Petal.Width^2)   1.7100     0.5487   3.116  0.00221 ** 
#> I(Petal.Width^3)  -0.5788     0.1408  -4.110 6.57e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3898 on 146 degrees of freedom
#> Multiple R-squared:  0.9522, Adjusted R-squared:  0.9512 
#> F-statistic: 969.9 on 3 and 146 DF,  p-value: < 2.2e-16

#Orthogonal
fit.orth <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), data = iris)

#Marginal effect of X at X=0 from orthogonal model
library(margins)
summary(margins(fit.orth, variables = "Petal.Width", 
                at = data.frame(Petal.Width = 0)))
#> Warning in check_values(data, at): A 'at' value for 'Petal.Width' is
#> outside observed data range (0.1,2.5)!
#>       factor Petal.Width    AME     SE      z      p  lower  upper
#>  Petal.Width      0.0000 1.1527 0.5836 1.9752 0.0482 0.0089 2.2965

Créé le 2019-10-25 par le package reprex (v0.3.0)

L'effet marginal de Petal.Widthà 0 de l'ajustement orthogonal et son erreur standard sont exactement égaux à ceux de l'ajustement polynomial brut. L'utilisation de polynômes orthogonaux n'améliore pas la précision des estimations de la même quantité entre les deux modèles.

OuiXOuiX

library(jtools)
data("iris")

fit.raw3 <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
fit.raw1 <- lm(Petal.Length ~ Petal.Width, data = iris)

round(summ(fit.raw3, part.corr = T)$coef, 3)
#>                    Est.  S.E. t val.     p partial.r part.r
#> (Intercept)       1.103 0.130  8.464 0.000        NA     NA
#> Petal.Width       1.153 0.584  1.975 0.050     0.161  0.036
#> I(Petal.Width^2)  1.710 0.549  3.116 0.002     0.250  0.056
#> I(Petal.Width^3) -0.579 0.141 -4.110 0.000    -0.322 -0.074

round(summ(fit.raw1, part.corr = T)$coef, 3)
#>              Est.  S.E. t val. p partial.r part.r
#> (Intercept) 1.084 0.073 14.850 0        NA     NA
#> Petal.Width 2.230 0.051 43.387 0     0.963  0.963

fit.orth3 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), 
               data = iris)
fit.orth1 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3)[,1], 
               data = iris)

round(summ(fit.orth3, part.corr = T)$coef, 3)
#>                                Est.  S.E.  t val. p partial.r part.r
#> (Intercept)                   3.758 0.032 118.071 0        NA     NA
#> stats::poly(Petal.Width, 3)1 20.748 0.390  53.225 0     0.975  0.963
#> stats::poly(Petal.Width, 3)2 -3.015 0.390  -7.735 0    -0.539 -0.140
#> stats::poly(Petal.Width, 3)3 -1.602 0.390  -4.110 0    -0.322 -0.074

round(summ(fit.orth1, part.corr = T)$coef, 3)
#>                                    Est.  S.E. t val. p partial.r part.r
#> (Intercept)                       3.758 0.039 96.247 0        NA     NA
#> stats::poly(Petal.Width, 3)[, 1] 20.748 0.478 43.387 0     0.963  0.963

Créé le 2019-10-25 par le package reprex (v0.3.0)

0,0010,0030,0050,9270,9270,0200,0050,927. À partir du modèle polynomial orthogonal mais pas du modèle polynomial brut, nous savons que la majeure partie de la variance expliquée dans le résultat est due au terme linéaire, avec très peu venant du terme carré et encore moins du terme cubique. Les valeurs polynomiales brutes ne racontent pas cette histoire.

Maintenant, que vous souhaitiez cet avantage d'interprétation par rapport à l'avantage d'interprétation d'être réellement capable de comprendre les coefficients du modèle, alors vous devriez utiliser des polynômes orthogonaux. Si vous préférez regarder les coefficients et savoir exactement ce qu'ils signifient (bien que je doute que ce soit généralement le cas), alors vous devriez utiliser les polynômes bruts. Si vous ne vous souciez pas (c.-à-d., Vous voulez seulement contrôler la confusion ou générer des valeurs prédites), alors cela n'a vraiment pas d'importance; les deux formulaires contiennent les mêmes informations concernant ces objectifs. Je dirais également que les polynômes orthogonaux devraient être préférés dans la régularisation (par exemple, le lasso), car la suppression des termes d'ordre supérieur n'affecte pas les coefficients des termes d'ordre inférieur, ce qui n'est pas vrai avec les polynômes bruts,

Noé
la source
1
Excellente contribution. Je ne peux pas répliquer vos résultats marginaux (la fonction margin affiche une erreur sur poly lorsque j'essaie d'exécuter votre premier bloc de code - je ne connais pas le package margin) - mais ils sont exactement ce que j'attends. Comme petite suggestion - vous devez également inclure la sortie de l'analyse de marge sur le modèle brut. Votre argument est miné (légèrement) par le changement des valeurs de p du résumé aux fonctions de marge (changeant nos conclusions non moins!) - qui semble être causé par l'utilisation d'une distribution normale au lieu d'une distribution. Votre point de régularisation est excellent.
user5957401
1
Merci pour les mots gentils. Je pense que vous devez inclure le stats::dans l'appel à poly()en lm()pour marginsle reconnaître ( ce qui est stupide). Je voulais concentrer mon argumentation sur les estimations ponctuelles et les erreurs standard, et je sais qu'il y a beaucoup d'informations étrangères et distrayantes présentées, mais j'espère que le texte illustre mes points.
Noah
Ce n'est pas ça. J'ai copié votre code exactement et vous utilisezstats::poly() . L'erreur dit 'degree' must be less than number of unique points- ce qui ne m'aide pas beaucoup. Néanmoins, margin()sauvegarde des déclarations prouvables donc ce n'est pas important.
user5957401
4

Je corrobore l'excellente réponse de @ user5957401 et j'ajoute des commentaires sur l'interpolation, l'extrapolation et les rapports.

Même dans le domaine des valeurs de paramètres stables, les coefficients / paramètres modélisés par les polynômes orthogonaux auront des erreurs standard sensiblement plus petites que les coefficients / paramètres modélisés par les paramètres bruts. Essentiellement, les polynômes orthogonaux sont un ensemble libre de descripteurs à covariance nulle. C'est PCA gratuitement!

Le seul inconvénient potentiel est de devoir l'expliquer à quelqu'un qui ne comprend pas la vertu des descripteurs à covariance nulle. Les coefficients ne sont pas immédiatement interprétables dans le contexte d'effets de premier ordre (de type vitesse) ou de second ordre (de type accélération). Cela peut être assez accablant dans un environnement professionnel.

Dans une simulation rapide avec un polynôme de degré 5, N = 1000 points, des valeurs de paramètres aléatoires (réduites de dix- pour maintenir la variable de réponse dans les 2 ordres de grandeur), et le bruit non corrélé la moitié de la variation globale de la variable de réponse, la R2Les deux modèles étaient identiques. Tout commeunej R2's. Le pouvoir prédictif est le même. Mais les erreurs standard des valeurs des paramètres pour le modèle orthogonal étaient égales ou des ordres de grandeur inférieurs au modèle brut.

Je serais donc "d'ordre de grandeur" plus confiant en rapportant le modèle orthogonal que le modèle brut. Dans la pratique, j'interpolerais avec l'un ou l'autre modèle, mais j'extrapolerais seulement avec le modèle orthogonal.

Peter Leopold
la source
1

J'aurais juste fait un commentaire pour le mentionner, mais je n'ai pas assez de représentants, je vais donc essayer de développer une réponse. Vous pourriez être intéressé de voir que dans la section de laboratoire 7.8.1 dans "Introduction à l'apprentissage statistique" (James et. Al., 2017, correction de la 8ème impression), ils discutent quelques différences entre l'utilisation ou non de polynômes orthogonaux, qui utilise le raw=TRUEou raw=FALSEdans la poly()fonction. Par exemple, les estimations de coefficient changeront, mais les valeurs ajustées ne:

# using the Wage dataset in the ISLR library
fit1 <- lm(wage ~ poly(age, 4, raw=FALSE), data=Wage)
fit2 <- lm(wage ~ poly(age, 4, raw=TRUE), data=Wage)
print(coef(fit1)) # coefficient estimates differ
print(coef(fit2))
all.equal(predict(fit1), predict(fit2)) #returns TRUE    

Le livre explique également comment, lorsque des polynômes orthogonaux sont utilisés, les valeurs de p obtenues à l'aide du anova()test F imbriqué (pour explorer le degré de polynôme pouvant être justifié) sont les mêmes que celles obtenues lors de l'utilisation du test t standard, produites par summary(fit). Ceci illustre que la statistique F est égale au carré de la statistique t dans certaines situations.

Shoeburg
la source
Les commentaires ne doivent jamais être utilisés comme réponses, quels que soient vos numéros de réputation.
Michael R. Chernick
En ce qui concerne votre dernier point, cela vaut également pour les polynômes non orthogonaux. Le coefficient t-test est égal au F-test comparant un modèle avec le coefficient en lui et un modèle sans pour tous les coefficients en régression (pris un à la fois).
Noah du