多項式回帰入門。線形回帰に飽きたらない人へ

『回帰分析』と言えば線形回帰を指すのが一般的です。
線形回帰とは得られたデータの列のど真ん中に直線をエイヤと引くことでデータのモデルを作ること、文字通り『線』を使う分析技術です。

図1:線形回帰のテスト

線形回帰は汎用性も高くモデルの解釈もしやすいとあって非常によく使われます。
しかしデータの間に直線関係=比例関係があると考えるのは最も単純な仮定です。
我々は現実はもっと複雑だとは知っているものの、あまりに複雑で手に終えそうもないので線形回帰を行っているというのが本音ではないでしょうか。
そうならば本当はもっと複雑なモデルを使えば分析の精度は向上し、より現実に近づけるのかもしれません。

【多項式回帰】

線形回帰を思い出します。
線形回帰はデータに対して次の予測モデルを当てはめて分析します。


\displaystyle y=\beta_0+\beta_1x

yが予測したい目的変数、xが予測に使える説明変数です。そして分析の基本方針は学習データ\{(x_1,y_1),(x_2,y_2),\ldots\}を用いて\beta_0,\beta_1がどのような値になるかを検討することです。

これをより複雑なモデルにしようとしたとき、いちばん自然なアイデアは


\displaystyle y=\beta_0+\beta_1x+\beta_2x^2\ \cdots (1)

と一次式だったものを二次式にしてしまうことです。
全く単純な手続きですが、たったこれだけで曲線を使った回帰ができるようになります。
パラメータも2つ(\beta_0,\beta_1)から3つ(\beta_0,\beta_1,\beta_2)となり、モデルの柔軟性も上がっていることが分かります。

ということは、これをさらに


\displaystyle y=\sum_{i=0}^{n}\beta_ix^i

とより複雑にすることも可能なはずです。
これが(n次)多項式回帰と呼ばれるモデルです。
nとして自由に自然数を選ぶことができて、複雑なモデルが必要になれば次数nを大きくすればよい、という非常に理解しやすいシステムになっています。

多項式回帰のシンプルさはx^iを独立した変数と見ると重回帰のように見えるということにもあります。
このため重回帰で使えたテクニック(係数の検定など)をそのまま使うことができますし、線形回帰と同じように最小二乗法でパラメータ推定ができます。

モデル式はシンプル、モデルの複雑さを変数一つで自由に調整でき、分析は簡単そして高速、それが多項式回帰です。
これ1つだけ覚えていれば単純な分析から複雑な分析まで自由になりそうです。
しかし本当にそんなうまい話あるのでしょうか。

【デモ用データ】

実際に分析してみましょう。
分析にはRにデフォルトで入っているデータcarsを使います。
carsは車の速度と停止距離との関係を調査したデータでサンプル数は50件です。
なお数字の単位はmph(マイル/時間)とft(フィート)です。

#R
head(cars)
#  speed dist
#1     4    2
#2     4   10
#3     7    4
#4     7   22
#5     8   16
#6     9   10
library(ggplot2)
gp <- ggplot(cars, aes(x = speed, y = dist)) + geom_point(size = 3) + xlim(0,25)
gp
図2:車の速度と停止距離との関係の調査データ

まず何も考えず線形回帰してみます。
速度がゼロのときは停止距離もゼロになると考えられるので切片がゼロの回帰をすることにします。
切片をゼロにするためにはモデル式に『-1』を追加します。

#R
lm1 <- lm(dist ~ speed - 1, cars)
summary(lm1)
#
#Call:
#lm(formula = dist ~ speed - 1, data = cars)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-26.183 -12.637  -5.455   4.590  50.181 
#
#Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
#speed   2.9091     0.1414   20.58   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 16.26 on 49 degrees of freedom
#Multiple R-squared:  0.8963,    Adjusted R-squared:  0.8942 
#F-statistic: 423.5 on 1 and 49 DF,  p-value: < 2.2e-16
gp + geom_smooth(method = lm, formula = y ~ x - 1, se = FALSE, fullrange = TRUE)
図3:carsを線形回帰

speedの係数は2.91と推定され、その推定標準誤差は\pm0.14です。
つまり車の速度が1マイル増えるごとに停止距離は2.91フィート増えるということです。
P値が十分小さい(***とアスタリスクが3つもついている)ので回帰モデルは有意です。
また誤差自乗平均(Residual Standard Error)は16.26とdistの平均値42.98に対して37%に抑えられています。
R-squaredは0.896で十分大きいでしょう。

と、線形回帰で割と満足のいく精度が出てしまいましたが、これを多項式の次数nを増やしながら分析を進化させていきます。

【二次多項式による回帰分析】

#R
lm2 <- lm(dist ~ poly(speed, degree = 2, raw = TRUE) - 1, cars)
summary(lm2)
#Call:
#lm(formula = dist ~ poly(speed, degree = 2, raw = TRUE) - 1, data = cars)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-28.836  -9.071  -3.152   4.570  44.986 
#
#Coefficients:
#                                     Estimate Std. Error t value Pr(>|t|)   
#poly(speed, degree = 2, raw = TRUE)1  1.23903    0.55997   2.213  0.03171 * 
#poly(speed, degree = 2, raw = TRUE)2  0.09014    0.02939   3.067  0.00355 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 15.02 on 48 degrees of freedom
#Multiple R-squared:  0.9133,    Adjusted R-squared:  0.9097 
#F-statistic: 252.8 on 2 and 48 DF,  p-value: < 2.2e-16
gp + geom_smooth(method = lm, formula = y ~ poly(x, degree = 2, raw = TRUE) - 1, se = FALSE, fullrange = TRUE)
図4:carsの二次多項式回帰

ここでR関数polyは多項式を生成する関数で、degreeが多項式の次数を指定するオプションです。
rawは生成する多項式を直交にするかどうかを指定するオプションで、デフォルトはTRUEになっています。FALSEにすることによって式(1)と同じ形で回帰できます。
多項式の直交性について説明することは本筋から離れるので別記事で説明しますが、基本的には多項式の次数(ここでは2)よりもデータサンプル数(ここでは50)が十分大きければどちらを使ってもあまり問題にならない、と思います。
唯一注意するべき点は、モデル式に-1を追加するためにはraw = TRUEとしている必要があるということです。

回帰結果に話を戻します。
まずグラフの見た目から明らかに回帰が改善しています。
またそれぞれの係数のp値が十分低く説明力は有意であることが確認できます。
誤差が15.02に小さくなっています。

【3次以上の回帰分析】

2次の回帰分析を使えるようになれば3次以上も簡単です。

#R
lm3 <- lm(dist ~ poly(speed, degree = 3, raw = TRUE) - 1, cars)
lm5 <- lm(dist ~ poly(speed, degree = 5, raw = TRUE) - 1, cars)
lm9 <- lm(dist ~ poly(speed, degree = 9, raw = TRUE) - 1, cars)
summary(lm3)
#
#Call:
#lm(formula = dist ~ poly(speed, degree = 3, raw = TRUE) - 1, data = cars)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-27.741  -8.755  -4.049   5.435  45.345 
#
#Coefficients:
#                                      Estimate Std. Error t value Pr(>|t|)
#poly(speed, degree = 3, raw = TRUE)1  2.299945   1.802672   1.276    0.208
#poly(speed, degree = 3, raw = TRUE)2 -0.038399   0.209557  -0.183    0.855
#poly(speed, degree = 3, raw = TRUE)3  0.003638   0.005871   0.620    0.539

#Residual standard error: 15.12 on 47 degrees of freedom
#Multiple R-squared:  0.914,     Adjusted R-squared:  0.9085 
#F-statistic: 166.5 on 3 and 47 DF,  p-value: < 2.2e-16
summary(lm5)
#
#Call:
#lm(formula = dist ~ poly(speed, degree = 5, raw = TRUE) - 1, data = cars)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-23.372  -8.174  -2.402   7.287  42.324 
#
#Coefficients:
#                                       Estimate Std. Error t value Pr(>|t|)
#poly(speed, degree = 5, raw = TRUE)1  4.2222943  9.0943265   0.464    0.645
#poly(speed, degree = 5, raw = TRUE)2 -1.2139330  2.6927510  -0.451    0.654
#poly(speed, degree = 5, raw = TRUE)3  0.1777469  0.2844085   0.625    0.535
#poly(speed, degree = 5, raw = TRUE)4 -0.0094557  0.0126653  -0.747    0.459
#poly(speed, degree = 5, raw = TRUE)5  0.0001711  0.0002015   0.849    0.400
#
#Residual standard error: 15.1 on 45 degrees of freedom
#Multiple R-squared:  0.9178,    Adjusted R-squared:  0.9087 
#F-statistic: 100.5 on 5 and 45 DF,  p-value: < 2.2e-16
summary(lm9)

#Call:
#lm(formula = dist ~ poly(speed, degree = 9, raw = TRUE) - 1, data = cars)
#
#Residuals:
#    Min      1Q  Median      3Q     Max 
#-21.829 -10.157  -1.091   5.990  43.232 
#
#Coefficients:
#                                       Estimate Std. Error t value Pr(>|t|)
#poly(speed, degree = 9, raw = TRUE)1 -2.958e+02  3.133e+02  -0.944    0.351
#poly(speed, degree = 9, raw = TRUE)2  2.220e+02  2.275e+02   0.976    0.335
#poly(speed, degree = 9, raw = TRUE)3 -6.742e+01  6.703e+01  -1.006    0.320
#poly(speed, degree = 9, raw = TRUE)4  1.102e+01  1.063e+01   1.037    0.306
#poly(speed, degree = 9, raw = TRUE)5 -1.070e+00  1.001e+00  -1.069    0.291
#poly(speed, degree = 9, raw = TRUE)6  6.368e-02  5.778e-02   1.102    0.277
#poly(speed, degree = 9, raw = TRUE)7 -2.279e-03  2.008e-03  -1.135    0.263
#poly(speed, degree = 9, raw = TRUE)8  4.505e-05  3.857e-05   1.168    0.250
#poly(speed, degree = 9, raw = TRUE)9 -3.778e-07  3.147e-07  -1.200    0.237
#
#Residual standard error: 15.23 on 41 degrees of freedom
#Multiple R-squared:  0.9239,    Adjusted R-squared:  0.9072 
#F-statistic:  55.3 on 9 and 41 DF,  p-value: < 2.2e-16
gp + geom_smooth(method = lm, formula = y ~ poly(x, degree = 3, raw = TRUE) - 1, se = FALSE, fullrange = TRUE)
gp + geom_smooth(method = lm, formula = y ~ poly(x, degree = 5, raw = TRUE) - 1, se = FALSE, fullrange = TRUE)
gp + geom_smooth(method = lm, formula = y ~ poly(x, degree = 9, raw = TRUE) - 1, se = FALSE, fullrange = TRUE)

図5:carsの3次多項式回帰

図6:carsの5次多項式回帰

図7:carsの9次多項式回帰

グラフを見る限り、9次では明らかに何か問題が発生しているようです。
パラメータ推定結果にも*(アスタリスク)が無くなってしまっています。
何が起こっているのでしょうか?
モデルが高度になるほどおかしな結果が出るとは納得がいかない、これでは勉強した甲斐が無いではないか!

これは機械学習で過学習(overfitting)と呼ばれる現象です。
モデルが複雑になると、あまりに柔軟に当てはめできるので、無理矢理に式を合わせようと機械が暴走しはじめるのです。
特にデータの端の部分でひどい。
端部分はデータが少ないので柔軟なモデルがめちゃくちゃグニャグニャしてしまえば誤差を思いっきり減らせるからです。
このことより柔軟なモデルを使った分析はデータが少なければ事態がより深刻になっていくと推定されます。

我々は複雑なモデルを手に入れても、それを使えば無条件に分析のレベルが上がるという話にはならないようです。
むしろ図7を見る限りレベルが上がるどころか悪化してしまう。
機械学習の世界ではno-free-lunch theorem(タダ飯はおまへん定理)ということが言われます。
複雑なモデルは複雑なりの弱点が有り、単純なモデルの方が優れた結果を出すこともありえるのです。

【次数の決定】

さて、過学習によりnを大きくしすぎると問題が起こることが分かりました。
しかし次数を増やすとモデルが柔軟になってより高度な回帰ができるのも事実です。
ということはnが十分小さいときはnを大きくすると性能は向上するが、性能はどこかで頭打ちになり、ある時点からnを大きくしすぎると性能は悪化していくと予想できます。

図8:多項式の次数と予測誤差?との関係図。次数1から2にかけては誤差が下がるがそれ以降は誤差はむしろ増える。さらにその増え方は加速する。

我々が知りたいのはnがいくつのときに最適になるかです。
どうやって最適なnを見つければよいでしょうか?

まず誤差を使うことは好ましくありません。
一般に次数を増やせば誤差は小さくなるからです。
これは二次多項式回帰の二次の係数をゼロにすれば線形回帰とみなせることから分かります。
重回帰分析で説明変数を増やせば増やすほど誤差が減っていくのと同じです。

シンプルな方法は係数のp値が有意かどうかを調べることです。
複雑なモデル(次数の多いモデル)の方が説明力が高いため、説明力の高いモデルですべてのp値が有意であればそちらを採用するべきでしょう。

もうひとつ伝統的な方法は分散分析を用いる方法です。
理論は省略しますが、Rでは次のように確認します。

#R
anova(lm1, lm2)
#Analysis of Variance Table
#
#Model 1: dist ~ speed - 1
#Model 2: dist ~ poly(speed, degree = 2, raw = TRUE) - 1
#  Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
#1     49 12954                                
#2     48 10831  1    2122.7 9.4069 0.003546 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

次数2のP値が0.354%ですので、次数2の係数を採用することで有意に説明力が改善するという結論です。

分散分析を用いてモデルの説明力を確認する場合は回帰がネストしていることが条件です。
上の例でいえば、二次多項式回帰は線形回帰に二次の項を追加したモデルであり、線形回帰を含んでいるため分散分析を用いることができます。

調子に乗って全部一気に分散分析で検定します。

#R
anova(lm1, lm2, lm3, lm5, lm9)
#Analysis of Variance Table
#
#Model 1: dist ~ speed - 1
#Model 2: dist ~ poly(speed, degree = 2, raw = TRUE) - 1
#Model 3: dist ~ poly(speed, degree = 3, raw = TRUE) - 1
#Model 4: dist ~ poly(speed, degree = 5, raw = TRUE) - 1
#Model 5: dist ~ poly(speed, degree = 9, raw = TRUE) - 1
#  Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
#1     49 12954                                
#2     48 10831  1   2122.66 9.1542 0.004272 **
#3     47 10743  1     87.75 0.3784 0.541845   
#4     45 10263  2    480.05 1.0351 0.364269   
#5     41  9507  4    756.35 0.8155 0.522711   
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

変数が追加されたことにより有意に説明力が向上しているときはアスタリスクが表示されます。
Model2=lm2=二次多項式までは説明力は増えるけれどもそれ以上はそうではありません、つまり二次多項式が採用すべきモデルだということです。

また有名な手法の1つにCross Validationというものもあります。
図8はCross Validationでそれぞれの次数の予測誤差を算出しました。
この手法にはいくつかの危険が潜んでいることも知られてはいますが、とりあえずで予測誤差を評価するのには役立ちます。
Cross Validationについてはいつか別記事で触れたいと思います。

【データの背景を見よ】

分散分析によって我々はめでたく最適なモデルは二次多項式回帰であることを確認することができました。
しかし実は二次多項式回帰を使うべき最大の理由は別のところにあるのです。

そもそもデータは自動車の速度と停止距離のデータでした。
物理学によれば自動車が持つ運動エネルギーは速度の二乗に比例します。
また停止距離は空走距離と制動距離の和で表すことができて、空走距離は車の速度に比例し、制動距離は車のエネルギーに比例します。
すると停止距離は速度の二次式であると考えられます。


停止距離 〜 速度^2

車の重量やブレーキの能力などの影響でここのデータのパラメータが異なりデータに誤差は発生しますが、二次式であることは変わりありません。

データを生み出す物理法則が二次式なのであればモデル式も二次式であるべきです。
そして二次回帰モデルのパラメータ\beta_0,\beta_1,\beta_2を観測データを使って推定することが統計分析の仕事になります。

データ分析をする際はまずデータの素性について考えることが分析の成果を向上させます。
データの素性を深く知れば統計テクニックに頼らずに統計モデルを作れます。
そして分散分析などの統計テクニックを使って求めた統計的に望ましいモデルよりも、原理に基づいて導き出したモデルの方が良いモデルです。
なぜならデータを使ってデータを生み出した原理を推測しようというのが統計学の技術なのであり、データを使わずに原理についての情報が得られるならばそれに越したことはないからです。
そこには有意水準など確率論が入ってくる余地はありません。

関連記事(非線形回帰)

関連記事(他)

参考文献

多項式回帰入門。線形回帰に飽きたらない人へ” への1件のコメント

  1. ピンバック:多項式回帰分析と重回帰分析のどちらを使えばいいか判断基準はありますか【メンターが回答】 | TechAcademyマガジン

コメントを残す