Jak interpretować współczynniki z dopasowania modelu wielomianowego?

Próbuję utworzyć wielomian dopasowania drugiego rzędu do niektórych danych, które mam. Powiedzmy, że knuję to dopasowanie z ggplot():

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

Dostaję:

Tak więc dopasowanie drugiego rzędu działa całkiem dobrze. Obliczam to za pomocą R:

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

I dostaję:

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 ***
# ....

Zakładam, że wzór na moje dopasowanie to:

Ale to po prostu daje mi złe wartości. Na przykład, gdy ma wartość 3, spodziewałbym się, że stanie się czymś w pobliżu 3.15. Jednak wstawiając do powyższej formuły otrzymuję: bar $\text{foo}$$\text{bar}$

Co daje? Czy niewłaściwie interpretuję współczynniki modelu?

Moja szczegółowa odpowiedź znajduje się poniżej, ale ogólna (tj. Prawdziwa) odpowiedź na tego rodzaju pytanie brzmi: 1) eksperymentuj, przekręć się, spójrz na dane, nie możesz uszkodzić komputera bez względu na to, co robisz, więc. . . eksperyment; lub 2) RTFM .

Oto Rkod, który mniej więcej odzwierciedla problem zidentyfikowany w tym pytaniu:

# 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))

Pierwszy lmzwraca oczekiwaną odpowiedź:

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

Drugi lmzwraca coś dziwnego:

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

Ponieważ lmjest tak samo w dwóch wywołaniach, muszą to być lmróżne argumenty . Spójrzmy więc na argumenty. Oczywiście yjest tak samo. To inne części. Spójrzmy na kilka pierwszych obserwacji zmiennych po prawej stronie w pierwszym wywołaniu lm. Powrót head(cbind(x,x^2))wygląda następująco:

            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

To jest zgodne z oczekiwaniami. Pierwsza kolumna to xdruga kolumna x^2. Co powiesz na drugie połączenie z lmtym z poli? Powrót head(poly(x,2))wygląda następująco:

              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, to naprawdę coś innego. Pierwsza kolumna nie jest x, a druga kolumna nie x^2. Cokolwiek więc poly(x,2)zrobi, nie powróci xi x^2. Jeśli chcemy wiedzieć, co polyrobi, możemy zacząć od przeczytania pliku pomocy. Tak mówimy help(poly). Opis mówi:

Zwraca lub ocenia ortogonalne wielomiany stopnia 1 do stopnia względem określonego zestawu punktów x. Wszystkie są prostopadłe do stałego wielomianu stopnia 0. Ewentualnie oceń surowe wielomiany.

Teraz albo wiesz, co to są „wielomiany ortogonalne”, albo nie. Jeśli nie, skorzystaj z Wikipedii lub Binga (oczywiście nie Google, ponieważ Google jest zły - nie tak zły jak Apple, oczywiście, ale nadal zły). Lub możesz zdecydować, że nie obchodzi Cię, jakie są wielomiany ortogonalne. Możesz zauważyć wyrażenie „surowe wielomiany” i możesz zauważyć nieco dalej w pliku pomocy polyz opcją, rawktóra domyślnie jest równa FALSE. Te dwa względy mogą zainspirować Cię do wypróbowania, head(poly(x, 2, raw=TRUE))które zwroty:

            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

Podekscytowany tym odkryciem (teraz wygląda na to, prawda?), Możesz spróbować. summary(lm(y ~ poly(x, 2, raw=TRUE))) Zwraca:

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

Powyższa odpowiedź ma co najmniej dwa poziomy. Najpierw odpowiedziałem na twoje pytanie. Po drugie i, co ważniejsze, zilustrowałem, w jaki sposób powinieneś sam odpowiadać na takie pytania. Każda osoba, która „umie programować”, przechodziła przez sekwencję taką jak ta ponad sześćdziesiąt milionów razy. Nawet ludzie tak przygnębiająco słabi w programowaniu jak ja cały czas przechodzę przez tę sekwencję. To normalne, że kod nie działa. Nieporozumienie, jakie funkcje pełnią, jest normalne. Sposobem poradzenia sobie z tym jest przekręcenie, eksperymentowanie, przeglądanie danych i RTFM. Wyjdź z trybu „bezmyślnego przestrzegania przepisu” i przejdź do trybu „detektywistycznego”.

Interesujące podejście do interpretacji regresji wielomianowej opracowali Stimson i in. (1978) . Wymaga przepisania

$Y = \beta_{0} + \beta_{1} X + \beta_{2} X^{2} + u$

tak jak

$Y = m + \beta_{2} \left( f - X \right)^{2} + u$

$m = \beta_{0} - \left. \beta_{1}^{2} \right/ 4 \beta_{2}$$\beta_{2}$$f = \left. -\beta_{1} \right/ 2 \beta_{2}$

Jeśli chcesz po prostu popchnąć we właściwym kierunku bez dość osądu: poly()tworzy wielomianów ortogonalnych (nie skorelowanych), w przeciwieństwie do I(), który całkowicie ignoruje korelację między wynikowymi wielomiany. Korelacja pomiędzy predyktorami może być problem w modelach liniowych (patrz tutaj , aby uzyskać więcej informacji na temat dlaczego korelacja może być problematyczne), więc to chyba lepiej (w ogóle), aby używać poly()zamiast I(). Dlaczego wyniki wyglądają tak inaczej? Cóż, zarówno poly()i I()podjąć X i przekształcić go w nowym X (w przypadku I(), nowa x jest tylko x ^ 1 lub x ^ 2, w przypadku poly(), nowy X są dużo bardziej skomplikowane (jeśli chcesz wiedzieć skąd pochodzą (a prawdopodobnie nie), możesz zacząćtutaj lub wyżej wymieniona strona Wikipedii lub podręcznik). Chodzi o to, że gdy obliczasz (przewidujesz) y na podstawie określonego zestawu wartości x, musisz użyć przekonwertowanych wartości x wytworzonych przez jeden z nich ( poly()lub w I()zależności od tego, która z nich była w twoim modelu liniowym). Więc:

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]

W takim przypadku oba modele zwracają tę samą odpowiedź, co sugeruje, że korelacja między zmiennymi predykcyjnymi nie wpływa na wyniki. Gdyby korelacja była problemem, dwie metody przewidywałyby różne wartości.

„poly” wykonuje orto-normalizację Grahama-Schmidta na wielomianach 1, x, x ^ 2, ..., x ^ deg Na przykład ta funkcja działa tak samo jak „poli” bez zwracania oczywiście atrybutów „coef”.

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]
}

Wylądowałem na tym wątku, ponieważ interesowała mnie forma funkcjonalna. Jak zatem wyrazić wynik „poli” jako wyrażenia? Po prostu odwróć procedurę Grahama-Schmidta. Skończy się bałagan!