【Rによるデータサイエンス】非線形回帰分析 〜 多項式回帰

前回「【Rによるデータサイエンス】非線形回帰分析　〜　ロジスティック回帰」では非線形回帰分析の方法の1つとしてロジスティック回帰について勉強した。

今回は2つ目の方法「多項式回帰」について勉強する。

次のようにして作成した人工データを使用する。

> x<-seq(-5,5,0.1) #-5から5までを当分割した長さ0.1のベクトルの作成
> y<-10*x^3+100*rnorm(x,0,1)
> plot(x,y)

さて、期待する多項式モデルは
である。
これを回帰モデルとして非線形回帰分析を実施する。

> fm3<-nls(y~a+b*x+c*x^2+d*x^3,start=c(a=1,b=1,c=1,d=1),trace=T)
21050827 :  1 1 1 1
1094019 :  -3.2969538  6.6961094  0.1130614  9.8168854
> summary(fm3)

Formula: y ~ a + b * x + c * x^2 + d * x^3

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  -3.2970    15.8523  -0.208    0.836    
b   6.6961     9.0645   0.739    0.462    
c   0.1131     1.3902   0.081    0.935    
d   9.8169     0.5431  18.076   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 106.2 on 97 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.069e-07

この結果を見ると、係数aのp値は0.836、係数cのp値は0.935であり、係数b,dのp値と比較して大きい。

それなら係数aとcの部分は無視して、

という多項式モデルを使用して非線形回帰分析を行ってみる。

> fm4<-nls(y~b*x+d*x^3,start=c(b=1,d=1),trace=T)
21032713 :  1 1
1094644 :  6.696107 9.816886
> summary(fm4)

Formula: y ~ b * x + d * x^3

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b   6.6961     8.9751   0.746    0.457    
d   9.8169     0.5377  18.256   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 105.2 on 99 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 2.788e-08

若干だが決定係数とp値の値が小さくなった。
じゃあ、回帰係数bの項も削除してみる。

このモデルで非線形回帰分析を実施する。

> fm5<-nls(y~d*x^3,start=c(d=1),trace=T)
21274942 :  1
1100799 :  10.18461
> summary(fm5)

Formula: y ~ d * x^3

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
d  10.1846     0.2145   47.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 104.9 on 100 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.007e-09

最後に、それぞれのAICを計算して変数の選択基準を調べる。

> AIC(fm3)
[1] 1234.941
> AIC(fm4)
[1] 1230.998
> AIC(fm5)
[1] 1229.565

ここまでの結果を整理すると次のようになる。

-	回帰係数4つの多項式モデル	回帰係数を2つにした簡易多項式モデル	回帰係数を1つにした簡易多項式モデル
決定係数	106.2	105.2	104.9
AIC	1234.941	1230.998	1229.565

可視化してみる。

> plot(x,y)
> lines(x,fitted(fm3),lty=1,col=2)
> lines(x,fitted(fm4),lty=2,col=4)
> lines(x,fitted(fm5),lty=3,col=6)
> legend(locator(1),c("多項式","係数2つの簡易版","係数1つの簡易版"),lty=1:3,col=c(2,4,6),lwd=2,cex=1.5)

小さくて何も見えねーなｗ。