R函數`poly`真的做了什麼？

我已閱讀手冊頁?poly（我承認我沒有完全理解），並閱讀Introduction to Statistical Learning中的函數說明。R函數`poly`真的做了什麼？

我目前的理解是撥打poly(horsepower, 2)應該等同於編寫horsepower + I(horsepower^2)。然而，這似乎是由下面的代碼的輸出相矛盾：

library(ISLR) 

summary(lm(mpg~poly(horsepower,2), data=Auto))$coef 

    #      Estimate Std. Error t value  Pr(>|t|) 
    #(Intercept)   23.44592 0.2209163 106.13030 2.752212e-289 
    #poly(horsepower, 2)1 -120.13774 4.3739206 -27.46683 4.169400e-93 
    #poly(horsepower, 2)2 44.08953 4.3739206 10.08009 2.196340e-21 

summary(lm(mpg~horsepower+I(horsepower^2), data=Auto))$coef 

    #     Estimate Std. Error t value  Pr(>|t|) 
    #(Intercept)  56.900099702 1.8004268063 31.60367 1.740911e-109 
    #horsepower  -0.466189630 0.0311246171 -14.97816 2.289429e-40 
    #I(horsepower^2) 0.0.0001220759 10.08009 2.196340e-21

我的問題是，爲什麼不輸出匹配，什麼是真正做poly？

來源

2013-10-20 merlin2011

看看在這個問題的答案：http://mathoverflow.net/questions/38864/visualizing-orthogonal-polynomials – Kieran

我只是看着（不完全理解），但我在這種情況下，d仍然想知道'poly（馬力，2）'產生的封閉形式公式究竟是什麼？ – merlin2011

嘗試'聚（馬力，度= 2，原始=真）';你傳遞2作爲錯誤的參數，'raw'默認爲FALSE。 – baptiste

當在統計模型中引入多項式項時，通常的動機是確定響應是否爲「彎曲」以及在添加該項時曲率是否「顯着」。投擲+I(x^2)項的結果是根據它們的位置，小的偏差可能被擬合過程「放大」，並且當它們僅僅是數據範圍的一端或另一端的波動時被誤解爲由於曲率項。這導致不適當地分配「重要性」的聲明。

如果你只是拋出一個與I(x^2)的平方項，必要性也將與x高度相關，至少在x > 0的域中。改爲使用：poly(x,2)會創建一個「曲線」變量集，其中線性項與x不具有高度相關性，曲率在數據範圍內大致相同。（如果您想了解統計理論，請搜索「正交多項式」。）只需鍵入poly(1:10, 2)並查看兩列。

poly(1:10, 2) 
       1   2 
[1,] -0.49543369 0.52223297 
[2,] -0.38533732 0.17407766 
[3,] -0.27524094 -0.08703883 
[4,] -0.16514456 -0.26111648 
[5,] -0.05504819 -0.34815531 
[6,] 0.05504819 -0.34815531 
[7,] 0.16514456 -0.26111648 
[8,] 0.27524094 -0.08703883 
[9,] 0.38533732 0.17407766 
[10,] 0.49543369 0.52223297 
attr(,"degree") 
[1] 1 2 
attr(,"coefs") 
attr(,"coefs")$alpha 
[1] 5.5 5.5 

attr(,"coefs")$norm2 
[1] 1.0 10.0 82.5 528.0 

attr(,"class") 
[1] "poly" "matrix"

的「二次」一詞的中心在5.5和線性項已被向下移動，因此在相同的x點是0（在模型中的隱式(Intercept)術語在被依賴用於換檔一切恢復當時預測被要求）。

來源

2013-10-20 23:59:58

這是我最喜歡的主題之一，我很高興地看到@ G.Grothendieck的回答，因爲我很佩服他的知識深度。他的答案優於我的，因爲他使用提供的數據集作爲問題中的「挑戰」。 –

要得到普通的多項式，就像在問題中使用raw = TRUE一樣。不幸的是，迴歸中的普通多項式存在一個不理想的方面。如果我們擬合二次曲線，比方說，然後是三次曲線，那麼三次曲線的低階係數全都不同於二次曲線的曲線，即56.900099702，-0.466189630,0.0，二次曲線對比6.068478e + 01，-5.688501e-01， 2.079011e-03用下面的立方體改裝後。

library(ISLR) 
fm2raw <- lm(mpg ~ poly(horsepower, 2, raw = TRUE), Auto) 
cbind(coef(fm2raw)) 
##           [,1] 
## (Intercept)      56.900099702 
## poly(horsepower, 2, raw = TRUE)1 -0.466189630 
## poly(horsepower, 2, raw = TRUE)2 0.0

fm3raw <- lm(mpg ~ poly(horsepower, 3, raw = TRUE), Auto) 
cbind(coef(fm3raw)) 
##           [,1] 
## (Intercept)      6.068478e+01 
## poly(horsepower, 3, raw = TRUE)1 -5.688501e-01 
## poly(horsepower, 3, raw = TRUE)2 2.079011e-03 
## poly(horsepower, 3, raw = TRUE)3 -2.146626e-06

我們真正喜歡的是添加在這樣一種方式，即採用二次適合低階係數保持不變，與三次擬合改裝後的三次項。爲此，採用poly(horsepower, 2, raw = TRUE)列的線性組合，並對poly(horsepower, 3, raw = TRUE)進行相同的操作，以使二次擬閤中的列彼此正交，並且類似於三次擬合。這足以保證當我們添加更高階係數時，低階係數不會改變。注意前三係數現在是如何在下面的兩套相同的（而高於它們的區別）：

fm2 <- lm(mpg ~ poly(horsepower, 2), Auto) 
cbind(coef(fm2)) 
##       [,1] 
## (Intercept)   23.44592 
## poly(horsepower, 2)1 -120.13774 
## poly(horsepower, 2)2 44.08953 
fm3 <- lm(mpg ~ poly(horsepower, 3), Auto) 
cbind(coef(fm3)) 
##        [,1] 
## (Intercept)   23.445918 
## poly(horsepower, 3)1 -120.137744 
## poly(horsepower, 3)2 44.089528 
## poly(horsepower, 3)3 -3.948849

重要的，因爲poly(horsepwer, 2)的列（的poly(horsepower, 2, raw = TRUE)兩種二次模型columnns只是線性組合正交和原始）表示相同的模型（即它們給出相同的預測）並且僅在參數化方面不同。例如，擬合值是相同的：

all.equal(fitted(fm2), fitted(fm2raw)) 
## [1] TRUE

我們也可以驗證多項式確實有正交列這也是正交的攔截：

nr <- nrow(Auto) 
e <- rep(1, nr)/sqrt(nr) # constant vector of unit length 
p <- cbind(e, poly(Auto$horsepower, 2)) 
round(crossprod(p), 2) 
##  1 2 
## 1 0 0 
## 1 0 1 0 
## 2 0 0 1

來源

2013-10-21 00:35:20

你能解釋爲什麼'e'被計算除以'nr'的平方根？ – HelloWorld

所以它的長度'sqrt（crossprod（e，e））'是1. –

以防萬一，[**以下問題**]（https://stats.stackexchange.com/questions/331607/is-多層次迴歸 - 適合於簡單 - 前後控制 - 設計）可以使用庫'lme4'來回答？ – rnorouzian

一個快速的回答是，poly一個向量是x，它基本上等於矩陣的QR分解，其列的功率爲x（居中之後）。例如：

> x<-rnorm(50) 
> x0<-sapply(1:5,function(z) x^z) 
> x0<-apply(x0,2,function(z) z-mean(z)) 
> x0<-qr.Q(qr(x0)) 
> cor(x0,poly(x,5)) 
       1    2    3    4    5 
[1,] -1.000000e+00 -1.113975e-16 -3.666033e-17 7.605615e-17 -1.395624e-17 
[2,] -3.812474e-17 1.000000e+00 1.173755e-16 -1.262333e-17 -3.988085e-17 
[3,] -7.543077e-17 -7.778452e-17 1.000000e+00 3.104693e-16 -8.472204e-17 
[4,] 1.722929e-17 -1.952572e-16 1.013803e-16 -1.000000e+00 -1.611815e-16 
[5,] -5.973583e-17 -1.623762e-18 9.163891e-17 -3.037121e-16 1.000000e+00

來源

2013-10-21 00:52:16 mrip

這可能比我接受的答案「數學基礎更深」。 –

R函數`poly`真的做了什麼？

回答

相關問題