Régression polynomiale brute ou orthogonale?

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. $y$$x,x^2,\ldots,x^5$

Pourquoi ne puis-je pas simplement faire une régression "normale" pour obtenir les coefficients de y = ∑ 5 i = 0 β i x $\beta_i$$y=\sum_{i=0}^5 \beta_i x^i$ (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?

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.

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

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

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?

$R^2$$X$$Y$$X=0$$X=0$$X$

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.

$Y$$X$$Y$$X$

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.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.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,

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 $10^{-d}$ pour maintenir la variable de réponse dans les 2 ordres de grandeur), et le bruit non corrélé $\sim$ la moitié de la variation globale de la variable de réponse, la $R^2$Les deux modèles étaient identiques. Tout comme$adj ~R^2$'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.

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.