Wie werden Koeffizienten aus einer Polynommodellanpassung interpretiert?

36

Ich versuche, ein Polynom zweiter Ordnung zu erstellen, das zu einigen meiner Daten passt. Angenommen, ich zeichne diese Übereinstimmung mit ggplot():

ggplot(data, aes(foo, bar)) + geom_point() + 
       geom_smooth(method="lm", formula=y~poly(x, 2))

Ich bekomme:

Parabelform mit Konfidenzband auf Streudiagramm

Eine Passung zweiter Ordnung funktioniert also ganz gut. Ich berechne es mit R:

summary(lm(data$bar ~ poly(data$foo, 2)))

Und ich bekomme:

lm(formula = data$bar ~ poly(data$foo, 2))
# ...
# Coefficients:
#                     Estimate Std. Error t value Pr(>|t|)    
# (Intercept)         3.268162   0.008282 394.623   <2e-16 ***
# poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206
# poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***
# ....

Nun würde ich annehmen, dass die Formel für meine Passform lautet:

bar=3.2680.122foo+1.575foo2

Aber das gibt mir nur die falschen Werte. Zum Beispiel, wenn 3 ist, würde ich erwarten, dass ungefähr 3.15 wird. Durch Einfügen in die obige Formel erhalte ich jedoch: Barfoobar

bar=3.2680.1223+1.57532=17.077

Was gibt? Interpretiere ich die Koeffizienten des Modells falsch?

user13907
quelle
2
Diese Frage wird in mehreren Threads beantwortet, die durch Durchsuchen unserer Website nach orthogonalem Polynom
whuber
6
@whuber Wenn ich gewusst hätte, dass das Problem bei "orthogonalen Polynomen" liegt, hätte ich wahrscheinlich eine Antwort gefunden. Aber wenn Sie nicht wissen, wonach Sie suchen sollen, ist es etwas schwierig.
user13907
2
Sie können auch nach Antworten suchen, indem Sie nach poly suchen , was in Ihrem Code an erster Stelle steht. Ich habe solche Informationen aus zwei Gründen in die Kommentare aufgenommen: (1) Die Links können zukünftigen Lesern und Ihnen helfen und (2) sie können Ihnen helfen, unser (etwas eigenwilliges) Suchsystem zu nutzen.
whuber
7
Sie haben eine Frage zu Ihrer Verwendung von gepostet, polyohne ?polyzuerst R einzugeben? Am oberen Rand steht in großen, freundlichen Buchstaben " Orthogonale Polynome berechnen ".
Glen_b
4
@Glen_b Ja, gut, ich habe Typ in ?polyder Syntax zu verstehen. Zugegeben, ich weiß nur wenig über die dahinter stehenden Konzepte. Ich wusste nicht, dass es noch etwas anderes gibt (oder einen so großen Unterschied zwischen "normalen" Polynomen und orthogonalen Polynomen), und die Beispiele, die ich online gesehen habe, wurden alle poly()zum Anpassen verwendet, insbesondere mit ggplot- warum sollte ich dann nicht einfach das und verwenden? verwirrt sein, wenn das Ergebnis "falsch" war? Wohlgemerkt, ich bin kein Mathematiker - ich wende nur das an, was ich von anderen gesehen habe, und versuche es zu verstehen.
user13907

Antworten:

55

Meine ausführliche Antwort ist unten, aber die allgemeine (dh echte) Antwort auf diese Art von Frage lautet: 1) Experimentieren Sie, schrauben Sie herum, schauen Sie sich die Daten an, Sie können den Computer nicht kaputt machen, egal was Sie tun. . . Experiment; oder 2) RTFM .

Hier ist ein RCode, der das in dieser Frage identifizierte Problem mehr oder weniger reproduziert:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/95939/
# 
# It is an exploration of why the result from lm(y_x+I(x^2))
# looks so different from the result from lm(y~poly(x,2))

library(ggplot2)


epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

summary(lm(y~x+I(x^2)))       # Looks right
summary(lm(y ~ poly(x, 2)))   # Looks like garbage

# What happened?
# What do x and x^2 look like:
head(cbind(x,x^2))

#What does poly(x,2) look like:
head(poly(x,2))

Der erste lmgibt die erwartete Antwort zurück:

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.92734    0.15376  25.542  < 2e-16 ***
x           -0.53929    0.11221  -4.806 5.62e-06 ***
I(x^2)       0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Der zweite lmgibt etwas Seltsames zurück:

Call:
lm(formula = y ~ poly(x, 2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.24489    0.02241 144.765  < 2e-16 ***
poly(x, 2)1  0.02853    0.22415   0.127    0.899    
poly(x, 2)2  1.09835    0.22415   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Da lmes in beiden Aufrufen dasselbe ist, müssen die Argumente lmunterschiedlich sein. Schauen wir uns also die Argumente an. Offensichtlich yist das auch so. Es sind die anderen Teile. Schauen wir uns die ersten Beobachtungen zu den Variablen auf der rechten Seite im ersten Aufruf von an lm. Die Rückkehr von head(cbind(x,x^2))sieht so aus:

            x         
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Das ist wie erwartet. Erste Spalte ist xund zweite Spalte ist x^2. Wie wäre es mit dem zweiten Anruf von lm, dem mit Poly? Die Rückkehr von head(poly(x,2))sieht so aus:

              1         2
[1,] -0.1714816 0.2169976
[2,] -0.1680173 0.2038462
[3,] -0.1645531 0.1909632
[4,] -0.1610888 0.1783486
[5,] -0.1576245 0.1660025
[6,] -0.1541602 0.1539247

OK, das ist wirklich anders. Die erste Spalte ist nicht xund die zweite Spalte ist nicht x^2. Also, was poly(x,2)auch immer tut, es kommt nicht zurück xund x^2. Wenn wir wissen wollen, was polypassiert, können wir zunächst die Hilfedatei lesen. Also sagen wir help(poly). Die Beschreibung lautet:

Gibt orthogonale Polynome vom Grad 1 bis zum Grad über der angegebenen Menge von Punkten x zurück oder wertet sie aus. Diese sind alle orthogonal zum konstanten Polynom vom Grad 0. Alternativ können Sie Rohpolynome auswerten.

Entweder wissen Sie, was "orthogonale Polynome" sind, oder Sie wissen es nicht. Wenn Sie dies nicht tun, verwenden Sie Wikipedia oder Bing (natürlich nicht Google, weil Google böse ist - natürlich nicht so schlecht wie Apple, aber immer noch schlecht). Oder Sie entscheiden, dass es Ihnen egal ist, was orthogonale Polynome sind. Möglicherweise bemerken Sie den Ausdruck "Rohpolynome" und etwas weiter unten in der Hilfedatei polyeine Option, rawdie standardmäßig gleich ist FALSE. Diese beiden Überlegungen könnten Sie dazu inspirieren, herauszufinden, head(poly(x, 2, raw=TRUE))welche Ergebnisse Sie erzielen:

            1        2
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Angeregt durch diese Entdeckung (es sieht richtig, jetzt, nicht wahr?), Könnten Sie gehen , um zu versuchen summary(lm(y ~ poly(x, 2, raw=TRUE))) Das gibt:

Call:
lm(formula = y ~ poly(x, 2, raw = TRUE))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              3.92734    0.15376  25.542  < 2e-16 ***
poly(x, 2, raw = TRUE)1 -0.53929    0.11221  -4.806 5.62e-06 ***
poly(x, 2, raw = TRUE)2  0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Es gibt mindestens zwei Ebenen für die obige Antwort. Zuerst habe ich deine Frage beantwortet. Zweitens, und noch viel wichtiger, habe ich Ihnen gezeigt, wie Sie solche Fragen selbst beantworten sollen. Jede einzelne Person, die "weiß, wie man programmiert", hat eine Sequenz wie die über sechzig Millionen Mal durchlaufen. Sogar Leute, die so deprimierend schlecht programmieren wie ich, gehen diese Sequenz die ganze Zeit durch. Es ist normal, dass Code nicht funktioniert. Es ist normal zu missverstehen, was Funktionen tun. Der Umgang damit besteht darin, herumzudrehen, zu experimentieren, die Daten und RTFM zu betrachten. Verlassen Sie den Modus "einem Rezept gedankenlos folgen" und wechseln Sie in den Modus "Detektiv".

Rechnung
quelle
7
Ich denke, das verdient eine +6. Ich werde versuchen, mich in ein paar Tagen zu erinnern, wann das möglich wird. FTR, ich denke, es muss nicht ganz so sarkastisch sein, aber es macht einen guten Job darin, zu zeigen, was orthogonale Polynome sind / wie sie funktionieren, und den Prozess zu zeigen, den Sie verwenden, um solche Dinge herauszufinden.
gung - Wiedereinsetzung von Monica
13
Tolle Antwort, danke. Ich bin zwar ein wenig beleidigt von einem "RTFM" (aber vielleicht bin es nur ich): Das Problem ist, dass in allem, was ich gelesen habe, zumindest im Hinblick auf die lineare Regression in R, manchmal Leute dies tun, andere das tun. Ehrlich gesagt verstehe ich den Wikipedia-Eintrag über orthogonale Polynome nicht. Mir fällt nicht ein, warum man dies für eine Regression verwenden würde, wenn die Koeffizienten, die Sie erhalten, "falsch" sind. Ich bin kein Mathematiker - ich versuche, die Rezepte zu befolgen, weil ich kein gelehrter Koch bin, aber ich muss trotzdem etwas essen.
user13907
12
@ user13907, das bist nicht nur du. Dies ist in der Tat eine gute Antwort, die eine Aufwertung verdient, aber es wäre von Vorteil, einen schöneren Ton zu haben.
Waldir Leoncio
8
Sie müssen nicht wirklich verstehen, was orthogonale Polynome hier sind - Sie müssen nur verstehen, dass sie nicht das sind, was Sie wollen. Warum möchte jemand orthogonale Polynome? Senden Sie cov (poly (x, 2)), um festzustellen, dass die Kovarianz zwischen den beiden Termen im Polynom Null ist (bis zum Rundungsfehler). Dies ist die Schlüsseleigenschaft orthogonaler Polynome - ihre Terme haben untereinander keine Kovarianz. Manchmal ist es praktisch, wenn Ihre RHS-Variablen keine Korrelation aufweisen. Ihre Koeffizienten sind nicht falsch, sie müssen nur unterschiedlich interpretiert werden.
Bill
2
Oh, okay, diese Erklärung in einfachem Englisch macht jetzt Sinn. Vielen Dank.
user13907
5

Es gibt einen interessanten Ansatz zur Interpretation der polynomialen Regression von Stimson et al. (1978) . Es beinhaltet das Umschreiben

Y=β0+β1X+β2X2+u

wie

Y=m+β2(fX)2+u

Dabei ist ist das Minimum oder Maximum (abhängig vom Vorzeichen von ) und ist der Fokuswert. Es transformiert im Grunde die dreidimensionale Kombination von Hängen in eine Parabel in zwei Dimensionen. Ihre Arbeit gibt ein Beispiel aus der Politikwissenschaft. β 2 f = - β 1 / 2 β 2m=β0β12/4β2β2f=β1/2β2

Durden
quelle
4

Wenn Sie nur einen Anstoß in die richtige Richtung wünschen, ohne so viel Urteilsvermögen: poly()Erstellen Sie stattdessen orthogonale (nicht korrelierte) Polynome, bei I()denen die Korrelation zwischen den resultierenden Polynomen vollständig ignoriert wird. Die Korrelation zwischen Prädiktorvariablen kann in linearen Modellen ein Problem sein (siehe hier für weitere Informationen darüber, warum die Korrelation problematisch sein kann). Daher ist es wahrscheinlich (im Allgemeinen) besser, sie poly()anstelle von zu verwenden I(). Warum sehen die Ergebnisse nun so anders aus? Nun, beide poly()und I()x nehmen und es in eine neue x umwandeln (im Fall I(), ist die neue x nur x ^ 1 oder x ^ 2, im Falle poly()werden die neuen xs viel komplizierter (wenn Sie wollen wissen , Woher sie kommen (und wahrscheinlich auch nicht), können Sie loslegenhier oder die oben genannte Wikipedia-Seite oder ein Lehrbuch). Der Punkt ist, dass Sie beim Berechnen (Vorhersagen) von y auf der Grundlage eines bestimmten Satzes von x-Werten die konvertierten x-Werte verwenden müssen, die entweder von poly()oder erzeugt wurden I()(je nachdem, welcher in Ihrem linearen Modell vorhanden war). So:

library(ggplot2)    

set.seed(3)
epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
   geom_smooth(method = "lm", formula = y ~ poly(x, 2))

modI <- lm(y~x+I(x^2)) 
summary(modI) # Looks right
modp <- lm(y ~ poly(x, 2))
summary(modp)  # Looks like garbage

# predict y using modI
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

# predict y using modp
# calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

In diesem Fall geben beide Modelle dieselbe Antwort zurück, was darauf hindeutet, dass die Korrelation zwischen Prädiktorvariablen Ihre Ergebnisse nicht beeinflusst. Wenn die Korrelation ein Problem wäre, würden die beiden Methoden unterschiedliche Werte vorhersagen.

filups21
quelle
1

'poly' führt eine Graham-Schmidt-Orthonormalisierung für die Polynome 1, x, x ^ 2, ..., x ^ deg durch. Diese Funktion macht beispielsweise das Gleiche wie 'poly', ohne natürlich die Attribute 'coef' zurückzugeben.

MyPoly <- 
function(x, deg)
{
    n <- length(x)
    ans <- NULL
    for(k in 1:deg)
    {
        v <- x^k
        cmps <- rep(0, n)
        if(k>0) for(j in 0:(k-1)) cmps <- cmps + c(v%*%ans[,j+1])*ans[,j+1]
        p <- v - cmps
        p <- p/sum(p^2)^0.5
        ans <- cbind(ans, p)
    }
    ans[,-1]
}

Ich bin auf diesem Thread gelandet, weil ich an der funktionalen Form interessiert war. Wie drücken wir das Ergebnis von 'poly' als Ausdruck aus? Kehren Sie einfach die Graham-Schmidt-Prozedur um. Sie werden mit einem Chaos enden!

izmirlig
quelle