59

私はマニュアルページを読み?poly（完全には理解していなかったことを認めます）、本Introduction to Statistical Learningの関数の説明も読みました。

私の現在の理解では、への呼び出しは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.001230536 0.0001220759  10.08009  2.196340e-21

私の質問は、出力が一致しないのはなぜですか、poly実際には何をしているのでしょうか?

4

3 に答える 3

69

生の多項式

質問のように通常の多項式を取得するには、 を使用しますraw = TRUE。残念ながら、回帰における通常の多項式には望ましくない側面があります。たとえば、2 次を当てはめ、次に 3 次を当てはめた場合、3 次の低次係数はすべて 2 次とは異なります。 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.001230536

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)、2 次近似の列が互いに直交するようにします。3 次近似についても同様です。これは、高次の係数を追加しても低次の係数が変わらないことを保証するのに十分です。最初の 3 つの係数が、下の 2 つのセットで同じになっていることに注意してください (上では異なります)。つまり、どちらの場合も、3 つの下位係数は 23.44592、-120.13774、および 44.08953 です。

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

同じ予測

重要なことに、の列は2 つの 2 次モデル (直交モデルと生モデル) の列のpoly(horsepwer, 2)線形結合にすぎないpoly(horsepower, 2, raw = TRUE)ため、同じモデル (つまり、同じ予測を与える) を表し、パラメーター化のみが異なります。たとえば、近似値は同じです。

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))
zapsmall(crossprod(p))
##   e 1 2
## e 1 0 0
## 1 0 1 0
## 2 0 0 1

残差二乗和

直交多項式のもう 1 つの優れた特性はpoly、列が単位長さを持ち、相互に直交する (切片列にも直交する) 行列を生成するという事実により、3 次項の追加による残差平方和の減少が単純になることです。モデル行列の 3 次列への応答ベクトルの射影の長さの 2 乗。

# these three give the same result

# 1. squared length of projection of y, i.e. Auto$mpg, on cubic term column
crossprod(model.matrix(fm3)[, 4], Auto$mpg)^2
##         [,1]
## [1,] 15.5934

# 2. difference in sums of squares
deviance(fm2) - deviance(fm3) 
## [1] 15.5934

# 3. difference in sums of squares from anova
anova(fm2, fm3) 
## 
## Analysis of Variance Table
## 
## Model 1: mpg ~ poly(horsepower, 2)
## Model 2: mpg ~ poly(horsepower, 3)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    389 7442.0                           
## 2    388 7426.4  1    15.593 0.8147 0.3673  <-- note Sum of Sq value
于 2013-10-21T00:35:20.013 に答える
36

統計モデルに多項式項を導入する場合、通常の動機は、応答が「曲がっている」かどうか、およびその項が追加されたときに曲率が「有意」であるかどうかを判断することです+I(x^2)。それらの位置に応じてフィッティングプロセスによって拡大され、曲率項が原因であると誤解され、データ範囲の一端または他端での変動に過ぎませんでした。これにより、「重要性」の宣言が不適切に割り当てられます。

の 2 乗項を入れるとI(x^2)、必然的に、少なくともx > 0. 代わりに使用:poly(x,2)線形項が x とあまり相関しておらず、曲率がデータの範囲全体でほぼ同じである変数の「曲がった」セットを作成します。(統計理論について詳しく知りたい場合は、「直交多項式」で検索してください。) 入力poly(1:10, 2)して 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-20T23:59:58.383 に答える
24

簡単な答えは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-21T00:52:16.527 に答える