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

36

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ę:

wykres parabolicznego dopasowania z pasmem pewności na wykresie rozrzutu

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:

bar=3.2680.122foo+1.575foo2

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 foofoobar

bar=3.2680.1223+1.57532=17.077

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

użytkownik13907
źródło
2
Na to pytanie można odpowiedzieć w kilku wątkach, które można znaleźć, przeszukując naszą stronę w poszukiwaniu ortogonalnego wielomianu
whuber
6
@ whuber Gdybym wiedział, że problem dotyczy „wielomianów ortogonalnych”, prawdopodobnie znalazłbym odpowiedź. Ale jeśli nie wiesz, czego szukać, jest to trochę trudne.
user13907
2
Możesz również znaleźć odpowiedzi, wyszukując na poli , który pojawia się w widocznym miejscu w twoim kodzie. Umieszczam takie informacje w komentarzach z dwóch powodów: (1) linki mogą pomóc przyszłym czytelnikom, a także tobie i (2) mogą pomóc ci pokazać, jak wykorzystać nasz (nieco dziwny) system wyszukiwania.
whuber
7
Zadałeś pytanie związane z używaniem polybez wpisywania ?polyR najpierw? Na górze napis „ Oblicz ortogonalne wielomiany ” dużymi przyjaznymi literami.
Glen_b
4
@Glen_b Tak, ja zrobiłem wpisać ?polyzrozumieć składnię. Trzeba przyznać, że mam niewielką wiedzę na temat koncepcji. Nie wiedziałem, że istnieje coś jeszcze (lub tak duża różnica między „normalnymi” wielomianami i ortogonalnymi wielomianami), a przykładami, które widziałem online, wszystkie używane poly()do dopasowania, szczególnie z ggplot- więc dlaczego nie miałbym tego użyć i być zdezorientowanym, jeśli wynik był „zły”? Pamiętaj, że nie mam umiejętności matematycznych - po prostu stosuję to, co widzieli inni, i staram się to zrozumieć.
user13907

Odpowiedzi:

55

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”.

Rachunek
źródło
7
Myślę, że to zasługuje na +6. Spróbuję przypomnieć sobie za kilka dni, kiedy stanie się to możliwe. FTR, myślę, że nie musi to być takie sarkastyczne, ale dobrze sprawdza się, pokazując, jakie są wielomiany ortogonalne / jak one działają, i pokazując proces, którego używasz do rozszyfrowania takich rzeczy.
gung - Przywróć Monikę
13
Świetna odpowiedź, dziękuję. Chociaż jestem trochę obrażony przez „RTFM” (ale może to tylko ja): Problemem jest to, że we wszystkim, co przeczytałem, przynajmniej w odniesieniu do regresji liniowej w R, ludzie czasem to robią, inni to robią. Szczerze mówiąc, nie rozumiem wpisu Wikipedii na temat wielomianów ortogonalnych. Nie przyszło mi do głowy, dlaczego można to wykorzystać do regresji, jeśli otrzymane współczynniki są „złe”. Nie jestem matematykiem - staram się przestrzegać przepisów, bo nie jestem uczonym kucharzem, ale mimo wszystko muszę coś zjeść.
user13907
12
@ user13907, to nie tylko ty. To rzeczywiście dobra odpowiedź, która zasługuje na głosowanie, ale skorzystałaby z ładniejszego tonu.
Waldir Leoncio
8
Naprawdę nie musisz rozumieć, jakie są wielomiany ortogonalne --- musisz po prostu zrozumieć, że nie są one tym, czego chcesz. Dlaczego ktoś może chcieć wielomianów ortogonalnych? Prześlij cov (poli (x, 2)), aby stwierdzić, że kowariancja między tymi dwoma członami wielomianu wynosi zero (do błędu zaokrąglenia). Jest to kluczowa właściwość wielomianów ortogonalnych - ich warunki mają zerową kowariancję względem siebie. Czasami wygodne jest, aby zmienne RHS miały zerową korelację ze sobą. Ich współczynniki nie są błędne, tak naprawdę należy je interpretować inaczej.
Bill
2
Och, okej, to wyjaśnienie zwykłym angielskim ma teraz sens. Dziękuję Ci.
user13907
5

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

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

tak jak

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

m=β0β12/4β2β2f=β1/2β2

Durden
źródło
2
+1 Aby zapoznać się z powiązanymi analizami, zobacz stats.stackexchange.com/questions/28730 i stats.stackexchange.com/questions/157629 .
whuber
4

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.

filups21
źródło
1

„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!

izmirlig
źródło