Rohe oder orthogonale Polynomregression?

22

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. yx,x2,,x5

Warum kann ich nicht einfach eine "normale" Regression durchführen, um die Koeffizienten von y = 5 i = 0 β i x i zu erhaltenβichy=ich=05βichxich (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?

l7ll7
quelle
1
@Sycorax Ich weiß, dass polydas etwas mit orthogonalen Polynomen zu tun hat und ich (x ^ 2) nicht (obwohl ich die Details nicht kenne) - aber warum würden die Autoren von ISLR dann eine Methode empfehlen, die nicht funktioniert? ? Es scheint sehr irreführend, wenn beide Befehle dasselbe tun, aber nur einer ist tatsächlich in Ordnung.
l7ll7
1
@gung Ich habe die Dokumentation von angeschaut polyund habe bereits eine Weile mit diesem Problem verbracht, aber ich kann nicht herausfinden, warum poly (x, 2) und x + I (x ^ 2) einen Unterschied machen. Könnten Sie mich bitte hier in den Kommentaren aufklären, wenn die Frage offtopisch ist?
l7ll7
1
@gung Ich habe meine Frage komplett überarbeitet. Diese Wahl von roh / orthogonal verwirrt mich noch mehr - früher dachte ich, dies sei nur eine geringfügige Rtechnische Tatsache, die ich nicht verstanden habe, aber jetzt scheint es ein ausgewachsenes Statistikproblem zu sein, das mich daran hindert, eine Regression zu codieren, die nicht sein sollte so schwer zu codieren.
17. Juli,
2
@gung Das hat mich eigentlich mehr verwirrt als geholfen. Früher dachte ich, ich sollte nur orthogonale Polynome verwenden, da dies der richtige Weg zu sein schien, aber in dieser Antwort werden rohe Polynome verwendet. Überraschenderweise schreit jeder im Netz "RTFM", aber es gibt tatsächlich keine klare Antwort, wann man was benutzt. (Ihr Link gibt auch keine Antwort darauf, nur ein Beispiel, wenn orth. Pol.
Schief
2
Sofern Sie nicht in einem physischen oder technischen Bereich arbeiten, in dem die Antwort als Quint-Polynom angegeben ist, ist es mit ziemlicher Sicherheit der richtige Ansatz, überhaupt keine polynomielle Regression durchzuführen. Investieren Sie Ihre Freiheitsgrade in einen Spline oder in etwas, das viel flexibler und stabiler ist als die Polynompassung.
Whuber

Antworten:

10

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.

user5957401
quelle
3
-1. Die größeren Standardfehler, die Sie auf den Koeffizienten niedrigerer Ordnung sehen, sind ein roter Hering. Die Koeffizienten niedrigerer Ordnung in Ihren beiden Modellen schätzen völlig verschiedene Dinge, daher macht der Vergleich ihrer Standardfehler keinen Sinn. Der Koeffizient höchster Ordnung ist der einzige, der in beiden Modellen dasselbe schätzt, und Sie werden feststellen, dass die t-Statistik identisch ist, unabhängig davon, ob die Polynome orthogonal sind oder nicht. Ihre beiden Modelle sind statistisch äquivalent in Bezug auf angepasste Werte, R ^ 2 usw., sie unterscheiden sich hauptsächlich in der Interpretation der Koeffizienten
Jake Westfall
@JakeWestfall, ich glaube nicht, dass ich dir zustimme. Zuallererst erzeugt das Ausführen des Codes Werte, die für alle Polynomreihenfolgen unterschiedlich sind, nicht für alle, sondern nur für eine - im Wesentlichen nimmt es das Polynom und führt PCA darauf aus. Zweitens, und was noch wichtiger ist, die t-Statistiken unterscheiden sich erheblich - wenn ich den Code in meiner Antwort anführe, wird dies bestätigt -, lösen wir das Multikollinearitätsproblem funktionell. Sie haben Recht, dass sich angepasste Werte, R ^ 2, F-Tests usw. nicht ändern. Das ist in der Tat ein Grund zur Orthogonalisierung - es ändert nichts außer unserer Fähigkeit, Polynomterme zu erkennen .
user5957401
1
Betreff: Der erste Punkt, sorry, ich wollte mich auf den t-stat des Terms höchster Ordnung beziehen, nicht auf seinen Koeffizienten. Dieser Prädiktor ist zwischen den Modellen skaliert + verschoben, also ändert sich die Coef, aber er testet den gleichen inhaltlichen Effekt, wie von t
Jake Westfall,
Betreff: Der zweite Punkt, der Grund dafür, dass "die t-Statistiken erheblich voneinander abweichen", ist wiederum, dass sie völlig unterschiedliche Dinge in den beiden Modellen schätzen. Betrachten Sie den linearen Effekt: raw.modEr schätzt die Steigung der Kurve bei x = 0, orthogonal.modschätzt die marginale Steigung (dh identisch mit lm(y ~ poly(x,1))dem Weglassen von Termen höherer Ordnung). Es gibt keinen Grund, warum Schätzungen dieser völlig unterschiedlichen Schätzer vergleichbare Standardfehler aufweisen sollten. Man kann leicht ein Gegenbeispiel konstruieren, in dem raw.modes viel höhere t-Werte gibt
Jake Westfall,
@JakeWestfall. Ich denke immer noch, Sie vermissen die Multikollinearität. Wir scheinen jedoch aneinander vorbeizusprechen, und es gibt vielleicht eine Lösung. Sie sagen, Sie können leicht ein Gegenbeispiel konstruieren, bitte. Ich denke, wenn ich das dgp sehe, an das du denkst, würde das viel für mich klären. Im Moment sind die einzigen Dinge, die mir eingefallen sind und die sich möglicherweise so verhalten, wie Sie es beschreiben, schwere Modellfehlspezifikationen.
User5957401
8

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

0,40,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?

Haitao Du
quelle
6
Sie verwenden anscheinend Schwimmer mit einfacher Genauigkeit und zitieren sie mit vierfacher Genauigkeit! Wie ist das passiert? Mit Ausnahme von GPUs wird für fast alle statistischen Berechnungen mindestens die doppelte Genauigkeit verwendet. ZB in Rder Ausgabe von print(0.4, digits=20)is 0.40000000000000002.
Whuber
6

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?

R2XY.X=0X=0erhalten 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 habenX

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.XY.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,0010,0030,0050,9270,9270,0200,0050,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.

Noah
quelle
1
Hervorragender Beitrag. Ich kann Ihre marginalen Ergebnisse nicht replizieren (die margin-Funktion gibt einen Fehler über poly aus, wenn ich versuche, Ihren ersten Codeblock auszuführen - ich kenne das margin-Paket nicht), aber sie sind genau das, was ich erwarte. Als kleiner Vorschlag - Sie sollten die Ausgabe der Margin-Analyse auch auf das Rohmodell anwenden. Ihr Argument wird (geringfügig) durch die Änderung der p-Werte von der Zusammenfassung zu den Randfunktionen (Änderung unserer Schlussfolgerungen nicht weniger!) Untergraben - was anscheinend durch die Verwendung einer normalen statt einer at-Verteilung verursacht wird. Ihr Regularisierungspunkt ist ausgezeichnet.
user5957401
1
Vielen Dank für die netten Worte. Ich glaube , Sie die umfassen müssen stats::im Aufruf poly()in lm()für marginssie zu erkennen (was dumm ist). Ich wollte mein Argument auf die Punktschätzungen und Standardfehler konzentrieren, und ich weiß, dass eine Menge irrelevanter und ablenkender Informationen präsentiert werden, aber ich hoffe, dass der Text meine Punkte veranschaulicht.
Noah,
Das ist es nicht. Ich habe Ihren Code genau kopiert und Sie verwenden stats::poly(). Der Fehler sagt 'degree' must be less than number of unique points- was mir nicht viel hilft. Trotzdem margin()ist es nicht wichtig, nachweisbare Aussagen zu sichern.
user5957401
4

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-dR2eindj R2 ‚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 .

Peter Leopold
quelle
1

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.

Schuheburg
quelle
Kommentare sollten niemals als Antworten verwendet werden, unabhängig von Ihrer Rufnummer.
Michael R. Chernick
In Bezug auf Ihren letzten Punkt gilt dies auch für nicht-orthogonale Polynome. Der Koeffiziententest ist gleich dem F-Test, der ein Modell mit dem darin enthaltenen Koeffizienten und ein Modell ohne Regressionskoeffizienten vergleicht (jeweils einzeln).
Noah,