【Rによるデータサイエンス】非線形回帰分析　〜　ロジスティック回帰

非線形回帰分析とは。。。
線形回帰分析以外のことらしいｗ。
まぁそこは深く突っ込む気はないので、その程度に捉えておく。

◎ロジスティック回帰
カラーテレビの普及率のデータを打ち込む。

> 年度<-c(1966:1984)
> 普及率<-c(0.003,0.016,0.054,0.139,0.263,0.423,0.611,0.758,0.859,0.903,0.937,0.954,0.978,0.978, 0.982,0.985,0.989,0.988,0.992)

可視化する。

> plot(年度,普及率)

横軸の区間[1966,1984]全体を見ると線形ではない関数に見えるので、非線形回帰分析を行ってみる。

経験的に、普及率や成長率のデータはこのような非線形曲線が多く、ロジスティック関数にあてはまるらしい。「ふーん」と思いながら深く考えずに、とりあえず指示通りやってみる。

ロジスティック関数は次の通り。

$y=\frac{a}{1+be^{cx}}$

ロジスティック関数を用いる回帰分析をロジスティック回帰と呼ぶ。

Rで非線形回帰分析を行うには関数nls*1を利用できる。

nls(formula,start,trace)

引数の説明は次の通り。

引数formula:被説明関数と説明変数との関係式
引数start:formulaにて記述した関係式に含まれる回帰係数（パラメータ）の名前と初期値
引数trace:パラメータの計算過程を返すか否かを指定

注意事項として、次のことがあげられる。

関数nlsでの回帰係数の決定には、線形回帰分析の場合と同様に、目的変数の実測値と予測値とのサを最小とする最小2乗法を使用するが、解がうまく求まらずに計算に失敗する場合がある。この場合、試行錯誤が必要になる。
それでは実際に非線形回帰分析を実施する。
> fm<-nls(普及率~a/(1+b*exp(c*年度)),start=c(a=1,b=1,c=-1),trace=TRUE)
 以下にエラー nlsModel(formula, mf, start, wts) : 
   パラメータの初期値で勾配行列が特異です
失敗した・・・。
失敗の原因は、関係式の中のexpが非常に大きな値になったことらしいが、どうすれば分かるんだろう・・・。
とりあえず教科書通りに、説明変数1966〜1984を1〜19に置き換えて実施する。
> fm<-nls(普及率~a/(1+b*exp(c*1:19)),start=c(a=1,b=1,c=-1),trace=TRUE)
3.905671 :   1  1 -1
2.387674 :   0.9824052  0.4300442 -0.1029666
1.743185 :   0.8872618  0.8264732 -0.2623701
0.7740848 :   0.9841109  2.3123040 -0.2310466
0.5578214 :   0.9271411  7.5149745 -0.5270324
0.09229081 :   1.0001338 17.1134431 -0.4350728
0.06874582 :   0.9606817 40.3886048 -0.6606493
0.01653944 :   0.9826601 75.9510196 -0.7160221
0.003486704 :    0.9806949 110.6878618  -0.7509771
0.001959816 :    0.9804580 123.8500779  -0.7565368
0.00194977 :    0.9806034 123.8048223  -0.7553703
0.001949752 :    0.9806268 123.6621028  -0.7551686
0.001949752 :    0.9806279 123.6609455  -0.7551647
サマリーを確認する。
> summary(fm)

Formula: 普及率 ~ a/(1 + b * exp(c * 1:19))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a   0.98063    0.00384 255.401  < 2e-16 ***
b 123.66094   13.56739   9.115 9.82e-08 ***
c  -0.75516    0.01742 -43.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01104 on 16 degrees of freedom

Number of iterations to convergence: 12 
Achieved convergence tolerance: 6.987e-06
係数はa=0.98063,b=123.66094,c=-0.75516と求められた。
決定係数（Residual standard error: 0.01104）は1に近いとはいえない。
可視化してみる。
予測値は関数fittedで得られる。

悪くない。
実測値を重ねてみる。

さて、今度は関数nlsの引数formulaに指定する関数に、Rが用意している既存の関数を指定して非線形回帰分析を実施してみる。
その場合、関数名を指定すればよく、関数を記述する手間が省ける。
ロジスティック関数SSlogisを試してみる。
>fm1<-nls(普及率~SSlogis(年度,Asym,xmid,scal))
> summary(fm1)

Formula: 普及率 ~ SSlogis(年度, Asym, xmid, scal)

Parameters:
      Estimate Std. Error  t value Pr(>|t|)    
Asym 9.806e-01  3.840e-03   255.40   <2e-16 ***
xmid 1.971e+03  3.533e-02 55791.87   <2e-16 ***
scal 1.324e+00  3.055e-02    43.35   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01104 on 16 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 6.041e-06
先ほど求めた数式を実際に打ち込んだ場合の結果fmと比較する。
> data.frame(fitted(fm),fitted(fm1))
   fitted.fm. fitted.fm1.
1  0.01658921  0.01658922
2  0.03464022  0.03464024
3  0.07088854  0.07088857
4  0.13947542  0.13947545
5  0.25576569  0.25576573
6  0.42053977  0.42053979
データを見ると見た目は凄く似ている。
相関係数を計算してみる。
> cor(fitted(fm),fitted(fm1))
[1] 1
Rの計算精度では相関係数の値は1。同一視できるね。

*1:Nonlinear Least Squares