Bagaimana menafsirkan koefisien dari model polinomial?

Saya mencoba membuat polinomial orde kedua untuk beberapa data yang saya miliki. Katakanlah saya merencanakan ini dengan ggplot():

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

Saya mendapat:

Jadi, urutan kedua cocok dengan cukup baik. Saya menghitungnya dengan R:

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

Dan saya mendapatkan:

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

Sekarang, saya akan menganggap formula untuk fit saya adalah:

Tapi itu hanya memberi saya nilai-nilai yang salah. Misalnya, dengan menjadi 3, saya harapkan menjadi sekitar 3,15. Namun, memasukkan ke rumus di atas saya dapatkan: $\text{foo}$$\text{bar}$

Apa yang menyebabkannya? Apakah saya salah menafsirkan koefisien model?

Jawaban terinci saya ada di bawah, tetapi jawaban umum (yaitu nyata) untuk pertanyaan seperti ini adalah: 1) bereksperimen, melihat-lihat, melihat data, Anda tidak dapat merusak komputer apa pun yang Anda lakukan, jadi. . . percobaan; atau 2) RTFM .

Berikut adalah beberapa Rkode yang mereplikasi masalah yang diidentifikasi dalam pertanyaan ini, kurang lebih:

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

Yang pertama lmmengembalikan jawaban yang diharapkan:

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

Yang kedua lmmengembalikan sesuatu yang aneh:

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

Karena lmsama dalam dua panggilan, itu harus menjadi argumen lmyang berbeda. Jadi, mari kita lihat argumennya. Jelas, ysama saja. Itu bagian lain. Mari kita lihat beberapa pengamatan pertama pada variabel sisi kanan pada panggilan pertama lm. Pengembalian head(cbind(x,x^2))terlihat seperti:

            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

Ini seperti yang diharapkan. Kolom pertama adalah xdan kolom kedua adalah x^2. Bagaimana dengan panggilan kedua lm, yang dengan poli? Pengembalian head(poly(x,2))terlihat seperti:

              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

Oke, itu sangat berbeda. Kolom pertama tidak x, dan kolom kedua tidak x^2. Jadi, apa pun yang poly(x,2)terjadi, itu tidak kembali xdan x^2. Jika kita ingin tahu apa yang polyterjadi, kita bisa mulai dengan membaca file bantuannya. Demikian kami katakan help(poly). Deskripsi mengatakan:

Mengembalikan atau mengevaluasi polinomial ortogonal dari derajat 1 ke derajat selama set poin yang ditentukan x. Ini semua ortogonal ke polinomial konstan derajat 0. Atau, mengevaluasi polinomial mentah.

Sekarang, apakah Anda tahu apa itu "polinomial ortogonal" atau tidak. Jika tidak, gunakan Wikipedia atau Bing (bukan Google, tentu saja, karena Google jahat --- tidak seburuk Apple, secara alami, tetapi masih buruk). Atau, Anda mungkin memutuskan bahwa Anda tidak peduli apa polinomial ortogonal itu. Anda mungkin memperhatikan frasa "polinomial mentah" dan Anda mungkin melihat sedikit lebih jauh di dalam file bantuan yang polymemiliki opsi rawyang, secara default, sama dengan FALSE. Dua pertimbangan tersebut dapat menginspirasi Anda untuk mencoba head(poly(x, 2, raw=TRUE))yang kembali:

            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

Gembira dengan penemuan ini (kelihatannya benar, sekarang, ya?), Anda mungkin mencoba summary(lm(y ~ poly(x, 2, raw=TRUE))) ini kembali:

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

Setidaknya ada dua level untuk jawaban di atas. Pertama, saya jawab pertanyaan Anda. Kedua, dan yang jauh lebih penting, saya menggambarkan bagaimana Anda seharusnya menjawab pertanyaan seperti ini sendiri. Setiap orang yang "tahu cara memprogram" telah melalui urutan seperti yang di atas enam puluh juta kali. Bahkan orang-orang yang sangat buruk dalam pemrograman seperti saya melalui urutan ini sepanjang waktu. Itu normal untuk kode tidak bekerja. Adalah normal untuk salah memahami fungsi apa yang dilakukan. Cara untuk mengatasinya adalah dengan bermain-main, bereksperimen, melihat data, dan RTFM. Keluarlah dari mode "tanpa resep" dan masuk ke mode "detektif".

Ada pendekatan yang menarik untuk interpretasi regresi polinomial oleh Stimson et al. (1978) . Ini melibatkan penulisan ulang

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

sebagai

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

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

Jika Anda hanya ingin dorongan ke arah yang benar tanpa penilaian terlalu banyak: poly()menciptakan polinomial ortogonal (tidak berkorelasi), yang bertentangan I(), yang sepenuhnya mengabaikan korelasi antara polinomial yang dihasilkan. Korelasi antara variabel prediktor dapat menjadi masalah dalam model linier (lihat di sini untuk informasi lebih lanjut tentang mengapa korelasi dapat menjadi masalah), jadi mungkin lebih baik (secara umum) untuk digunakan poly()daripada menggunakan I(). Sekarang, mengapa hasilnya terlihat sangat berbeda? Baik, keduanya poly()dan I()ambil x dan ubah menjadi x baru (dalam kasus I(), x baru hanya x ^ 1 atau x ^ 2, dalam kasus poly(), x baru jauh lebih rumit (jika Anda ingin tahu dari mana mereka berasal (dan Anda mungkin tidak), Anda bisa memulaidi sini atau halaman Wikipedia yang disebutkan di atas atau buku teks). Intinya adalah, ketika Anda menghitung (memprediksi) y berdasarkan sekumpulan nilai x tertentu, Anda perlu menggunakan nilai x yang dikonversi yang dihasilkan oleh salah satu poly()atau I()(tergantung mana yang ada dalam model linier Anda). Begitu:

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]

Dalam hal ini, kedua model mengembalikan jawaban yang sama, yang menunjukkan bahwa korelasi antara variabel prediktor tidak mempengaruhi hasil Anda. Jika korelasi merupakan masalah, kedua metode akan memprediksi nilai yang berbeda.

'poli' melakukan Graham-Schmidt orto-normalisasi pada polinomial 1, x, x ^ 2, ..., x ^ deg Misalnya fungsi ini melakukan hal yang sama seperti 'poli' tanpa mengembalikan atribut 'coef' tentunya.

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

Saya mendarat di utas ini karena saya tertarik pada bentuk fungsional. Jadi bagaimana kita mengekspresikan hasil 'poli' sebagai ekspresi? Balikkan prosedur Graham-Schmidt. Anda akan berakhir berantakan!