Rohe oder orthogonale Polynomregression?

Ich möchte eine Variable auf x , x 2 , … , x 5 zurückführen . Sollte ich dazu rohe oder orthogonale Polynome verwenden? Ich habe mir Fragen auf der Website angesehen, die sich mit diesen beschäftigen, aber ich verstehe nicht wirklich, was der Unterschied zwischen der Verwendung ist. $y$$x,x^2,\ldots,x^5$

Warum kann ich nicht einfach eine "normale" Regression durchführen, um die Koeffizienten von y = ∑ 5 i = 0 β i x $\beta_i$$y=\sum_{i=0}^5 \beta_i x^i$ (zusammen mit p-Werten und all den anderen netten Dingen) und muss mir stattdessen Sorgen machen, ob ich sie verwende? rohe oder orthogonale Polynome? Diese Wahl scheint mir außerhalb des Rahmens dessen zu liegen, was ich tun möchte.

In dem Statistikbuch, das ich gerade lese (ISLR von Tibshirani et al.), Wurden diese Dinge nicht erwähnt. Eigentlich wurden sie in gewisser Weise heruntergespielt.
Der Grund ist, AFAIK, dass in der lm()Funktion in R die Verwendung von y ~ poly(x, 2)Beträgen der Verwendung von orthogonalen Polynomen und die Verwendung von y ~ x + I(x^2)Beträgen der Verwendung von rohen Polynomen entspricht. Aber auf S. 116 sagen die Autoren, dass wir die erste Option verwenden, weil die letztere "umständlich" ist, was keinen Hinweis darauf hinterlässt, dass diese Befehle tatsächlich völlig verschiedene Dinge betreffen (und folglich unterschiedliche Ausgaben haben).
(dritte Frage) Warum würden die Autoren von ISLR ihre Leser so verwirren?

Ich glaube, bei der Antwort geht es weniger um numerische Stabilität (obwohl dies eine Rolle spielt) als vielmehr darum, die Korrelation zu verringern.

Im Wesentlichen läuft das Problem darauf hinaus, dass die Kovariaten, gegen die wir regressieren, in hohem Maße korrelieren, wenn wir uns gegen eine Reihe von Polynomen höherer Ordnung zurückbilden. Beispielcode unten:

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

Das ist enorm wichtig. Mit zunehmender Korrelation der Kovariaten schwindet unsere Fähigkeit, zu bestimmen, welche wichtig sind (und wie groß ihre Auswirkungen sind), rapide. Dies wird typischerweise als das Problem der Multikollinearität bezeichnet. Wenn wir zwei Variablen hatten, die vollständig korreliert waren, und wenn wir sie gegen etwas regressieren, ist es unmöglich, zwischen den beiden zu unterscheiden - Sie können sich dies als eine extreme Version des Problems vorstellen, aber dieses Problem wirkt sich auf unsere Schätzungen für aus auch geringere Korrelationsgrade. Im eigentlichen Sinne - auch wenn numerische Instabilität kein Problem war - beschädigt die Korrelation aus Polynomen höherer Ordnung unsere Inferenzroutinen enorm. Dies äußert sich in größeren Standardfehlern (und damit kleineren t-Werten), die Sie sonst sehen würden (siehe Beispiel für eine Regression unten).

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)

Wenn Sie diesen Code ausführen, ist die Interpretation ein wenig schwierig, da sich die Koeffizienten alle ändern und es daher schwierig ist, die Ergebnisse zu vergleichen. Wenn wir uns die T-Statistiken ansehen, können wir sehen, dass die Fähigkeit, die Koeffizienten zu bestimmen, bei orthogonalen Polynomen VIEL größer war. Für die 3 relevanten Koeffizienten erhielt ich t-Werte von (560,21,449) für das Orthogonalmodell und nur (28, -38,121) für das Rohpolynommodell. Dies ist ein großer Unterschied für ein einfaches Modell, bei dem nur wenige Polynomterme niedriger Ordnung von Bedeutung sind.

Das heißt nicht, dass dies ohne Kosten kommt. Es sind zwei Hauptkosten zu berücksichtigen. 1) Wir verlieren etwas Interpretierbarkeit mit orthogonalen Polynomen. Wir verstehen vielleicht, was der Koeffizient auf x**3bedeutet, aber die Interpretation des Koeffizienten auf x**3-3x(das dritte Einsiedler-Poly - nicht unbedingt das, was Sie verwenden werden) kann viel schwieriger sein. Zweitens - wenn wir sagen, dass diese Polynome orthogonal sind - meinen wir, dass sie in Bezug auf ein gewisses Maß an Distanz orthogonal sind. Es kann schwierig sein, ein für Ihre Situation relevantes Entfernungsmaß auszuwählen. Allerdings glaube ich, dass die polyFunktion so gewählt werden soll, dass sie in Bezug auf die Kovarianz orthogonal ist - was für lineare Regressionen nützlich ist.

Warum kann ich nicht einfach eine "normale" Regression durchführen, um die Koeffizienten zu erhalten?

$0.4$$0.4000000059604644775390625$ hier

Die Verwendung eines rohen Polynoms führt zu Problemen, da wir eine große Anzahl haben werden. Hier ist ein kleiner Beweis: Wir vergleichen die Matrixbedingungsnummer mit dem rohen und dem orthogonalen Polynom.

> 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

Sie können meine Antwort auch hier für ein Beispiel überprüfen.

Warum gibt es große Koeffizienten für Polynome höherer Ordnung?

Ich habe das Gefühl, dass einige dieser Antworten völlig verfehlen. Die Antwort von Haitao befasst sich mit den Rechenproblemen beim Anpassen von Rohpolynomen, aber es ist klar, dass OP nach dem fragt statistischen Unterschieden zwischen den beiden Ansätzen . Wenn wir also einen perfekten Computer hätten, der alle Werte genau darstellen könnte, warum würden wir dann einen Ansatz dem anderen vorziehen?

$R^2$$X$$Y$$X=0$$X=0$erhalten Sie genau die gleiche Steigung und den gleichen Standardfehler, obwohl sich der Koeffizient und der Standardfehler des Terms erster Ordnung in der orthogonalen Polynomregression vollständig von seinem Wert in der rohen Polynomregression unterscheiden. Das heißt, wenn versucht wird, aus beiden Regressionen die gleichen Größen zu erhalten (dh Größen, die auf die gleiche Weise interpretiert werden können), sind die Schätzungen und Standardfehler identisch. Die Verwendung von orthogonalen Polynomen bedeutet nicht, dass Sie auf magische Weise mehr Gewissheit über die Steigung von haben$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

^{Erstellt am 25.10.2019 durch das Paket reprex (v0.3.0)}

Der marginale Effekt von Petal.Widthat 0 aus der orthogonalen Anpassung und sein Standardfehler sind genau gleich denen aus der rohen Polynomanpassung. Die Verwendung von orthogonalen Polynomen verbessert nicht die Genauigkeit von Schätzungen derselben Größe zwischen den beiden Modellen.

$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

^{Erstellt am 25.10.2019 durch das Paket reprex (v0.3.0)}

$0.001$$0.003$$0.005$$0.927$$0.927$$0.020$$0.005$$0.927$. Aus dem orthogonalen Polynommodell, aber nicht dem rohen Polynommodell, wissen wir, dass der größte Teil der im Ergebnis erklärten Varianz auf dem linearen Term beruht, wobei nur sehr wenig vom quadratischen Term und noch weniger vom kubischen Term kommt. Die rohen Polynomwerte erzählen diese Geschichte nicht.

Wenn Sie nun diesen Interpretationsvorteil gegenüber dem Interpetationsvorteil der tatsächlichen Fähigkeit, die Koeffizienten des Modells zu verstehen, wünschen, sollten Sie orthogonale Polynome verwenden. Wenn Sie lieber die Koeffizienten betrachten und genau wissen möchten, was sie bedeuten (obwohl ich bezweifle, dass dies normalerweise der Fall ist), sollten Sie die rohen Polynome verwenden. Wenn es Ihnen egal ist (dh Sie möchten nur die Verwechslung kontrollieren oder vorhergesagte Werte generieren), ist dies wirklich egal. Beide Formulare enthalten dieselben Informationen zu diesen Zielen. Ich würde auch argumentieren, dass orthogonale Polynome bei der Regularisierung bevorzugt werden sollten (z. B. Lasso), da das Entfernen von Termen höherer Ordnung die Koeffizienten der Termen niedrigerer Ordnung nicht beeinflusst, was bei rohen Polynomen nicht zutrifft.

Ich bestätige die ausgezeichnete Antwort von @ user5957401 und füge Kommentare zu Interpolation, Extrapolation und Berichterstellung hinzu.

Selbst im Bereich stabiler Parameterwerte weisen die durch die orthogonalen Polynome modellierten Koeffizienten / Parameter wesentlich kleinere Standardfehler auf als die durch die Rohparameter modellierten Koeffizienten / Parameter. Im Wesentlichen sind die orthogonalen Polynome eine freie Menge von Null-Kovarianz-Deskriptoren. Das ist PCA kostenlos!

Der einzige mögliche Nachteil besteht darin, dies jemandem erklären zu müssen, der die Tugend der Null-Kovarianz-Deskriptoren nicht versteht. Die Koeffizienten sind im Kontext von Effekten erster Ordnung (geschwindigkeitsabhängig) oder zweiter Ordnung (beschleunigungsabhängig) nicht sofort interpretierbar. Dies kann in einem Geschäftsumfeld ziemlich schädlich sein.

$10^{-d}$$\sim$$R^2$$adj ~R^2$ ‚s. Die Vorhersagekraft ist dieselbe. Die Standardfehler der Parameterwerte für das orthogonale Modell waren jedoch gleich oder um Größenordnungen niedriger als das Rohmodell.

Ich würde also "Größenordnungen" sicherer sein, das orthogonale Modell zu beschreiben als das rohe. In der Praxis würde ich mit beiden Modellen interpolieren , aber ich würde nur mit dem orthogonalen extrapolieren .

Ich hätte das nur kommentiert, aber ich habe nicht genug Repräsentanten, also werde ich versuchen, eine Antwort zu finden. Es könnte Sie interessieren zu sehen, dass in Laborabschnitt 7.8.1 in "Einführung in das statistische Lernen" (James et al., 2017, korrigierter achter Ausdruck) einige Unterschiede zwischen der Verwendung von orthogonalen Polynomen und der Verwendung von nicht-orthogonalen Polynomen erörtert werden raw=TRUEoder raw=FALSEin der poly()Funktion. Beispielsweise ändern sich die Koeffizientenschätzungen, die angepassten Werte jedoch nicht:

# 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    

In dem Buch wird auch erläutert, wie bei Verwendung von orthogonalen Polynomen die mit dem anova()verschachtelten F-Test (um zu untersuchen, inwieweit ein Polynom gerechtfertigt sein kann) erhaltenen p-Werte mit denen des Standard-t-Tests von identisch sind summary(fit). Dies zeigt, dass die F-Statistik in bestimmten Situationen dem Quadrat der t-Statistik entspricht.