1．非線形回帰分析とは
線形モデルは被説明変数を説明変数の線形関数で表現するモデルである。本稿では、線形モデル以外のモデルを非線形モデルと呼ぶことにする。
非線形回帰分析には、説明変数や被説明変数に何らかの変換を施し、線形関係に置き換え線形回帰分析を行う方法、非線形の関数を当てはめる方法、加法モデル、樹木モデル、ニューラルネットワーク法など多数の方法が提案されている。
2．非線形回帰モデルと関数nls
非線形回帰分析の方法の最も基本的な方法は、非線形の関数式を用いてデータを当てはめる方法である。
Rには、自由に関数式を指定することができる非線形回帰分析の関数nlsがある。関数nlsの簡単な書き式を次に示す。
nls(formula,data,start,trace)
引数 formulaの書き式は、関数lmやglmと若干異なる。関数lmとglmでは非説明変数と説明変数のみを指定すればよかったが、関数nlsでは被説明関数と説明変数との関係式を具体的に書く必要がある。
例えば、データに関数
を当てはめる非線形回帰分析の場合はfourmulaの書き式は次のように関係式を書く。
y~ a/(1+b*exp(c*x))
このfourmulaの中のa、b、cが求める回帰係数(パラメータ)である。
関数nlsでは、線形回帰分析の場合と同じく、被説明変数の実測値と予測値(あるいは推測値)との差を最小とする最小2乗法で係数(パラメータ)を求めている。
関数nlsでは、パラメータの名前とその初期値を指定する必要がある。初期値の指定は、引数startで設定した方が良い。上記のformulaに対応するstartの記述例を次に示す。
nls(y~ a/(1+b*exp(c*x)),start=c(a=1,b=1,c=-1))
初期値は、データ解析者の勘と経験に頼るのが殆どである。
traceはパラメータの計算の過程を返すか否かを指定する引数である。計算過程の結果を返す場合は、TRUE(あるいはT)を指定する。デフォルトはFALSEになっている。
非線形回帰分析は、線形回帰分析より計算が難しいため、係数パラメータを求める際に、数値解がうまく求められず、プログラム実行が失敗する場合がしばしばある。その際には、モデルを換えて試行錯誤しなければならない。
ここでは一般化線モデルを説明する際に用いた日本におけるカラーテレビの普及率(先月の号を参照)の例を用いて非線形回帰関数nlsの使用方法を説明する。
>年度<-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)
このような普及率、成長率のデータは、通常ロジスティック関数が多く用いられている。ここでは、次に示すロジスティック関数を用いてみよう。
上記に説明した手順で、次のように指定して実行するとエラメッセージが返される。
>fm<-nls(普及率 ~ a/(1+b*exp(c*年度)),start=c(a=1,b=1,c=-1),trace
= TRUE)
これは、関係式の中の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.7740847
: 0.9841109
2.3123040 -0.2310466
0.5578214
: 0.9271411 7.5149748 -0.5270324
0.09229084
: 1.0001338
17.1134434 -0.4350728
0.06874581
: 0.9606817
40.3886048 -0.6606493
0.01653944
: 0.9826601
75.9510234 -0.7160221
0.003486704
: 0.9806949 110.6878664
-0.7509771
0.001959816
: 0.9804580 123.8500804 -0.7565368
0.001949770
: 0.9806034 123.8048180 -0.7553703
0.001949752
: 0.9806268 123.6621023 -0.7551686
0.001949752
: 0.9806279 123.6609507 -0.7551648
このようにパラメータの計算過程が返され、パラメータa、b、cの最終推定値はそれぞれ0.9807489、123.54058867、-0.7549528となっている。関数summaryで回帰の要約が返される。
>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.66095
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
Correlation of Parameter Estimates:
a b
b -0.2950
c 0.3967
-0.9709
予測値は関数lmの場合と同じく関数fittedで返される。
次のコマンドで、実測値と非線形回帰関数nlsによる予測値をグラフで示す。図1から、予測値は実測値に非常によく近似していることが分かる。
>plot(年度,普及率,cex=2)
>lines(年度,普及率)
>lines(年度,fitted(fm),col=2,lty=2,lwd=2)
>legend(1975,0.5,c("実測値","予測値"),
col=1:2,lty=1:2,lwd=2)
図１ テレビの普及率の実測値と予測値
残差は関数residで返され、関数plotで予測値対標準化された残差の散布図が作成される。
Rには頻繁に使用されている関数式が用意されているものもある。関数nlsでは、その関数式をリンクして、非線形回帰分析を行うことができる。用意されている関数をリンクして回帰分析を行う際には、その関数の具体的な関数式を記述する必要がなく、関数のオブジェクトの名前を引数として用いる。例えば、ロジスティック関数のオブジェクトの名前はSSlogis、ワイブール関数はSSweibullとなっている。
カラーテレビの普及率のデータをロジスティック関数SSlogisを用いる場合は、次のような書き式で非線形回帰分析ができる。
>fm1<-nls(普及率~SSlogis(1:19, Asym, xmid, scal ),list( Asym = 1, xmid = 0, scal = 1
),trace = TRUE)
上記の具体的な関係式を記述して求めた予測値とロジスティック関数SSlogisを用いて求めた予測値を次に示す。異なる両方法で得られた予測値は非常によく近似している。
> data.frame(fitted(fm),fitted(fm1))
fitted.fm
fitted.fm1
1
0.01660369 0.01660370
2
0.03466288 0.03466289
＜中略＞
18 0.98059694 0.98059694
19 0.98067744 0.98067745
関数nlsによる非線形回帰分析の説明のため、人工データを用いた例を次に示す。まず次のような人工データを作成し、そのデータの散布図を図2に示す。
>x<-seq(-5, 5,0.1)
>y<-10*x^3+100*rnorm(x,0,1)
>plot(x,y)
図2 3次多項式の人口データ
データの散布図から縦軸yと横軸xの関係は3次多項式の関係であることが分かる。
3次多項式の一般式は
である。この式を用いたnlsは次ぎのように記述する。
>fm3=nls(y~a+b*x+c*x^2+d*x^3,start=c(a=1,b=1,c=1,d=1),trace=T)
20490809
: 1 1 1 1
1089805
: -20.576319 24.692240 1.206357
8.678613
> summary(fm3)
Formula: y ~
a + b * x + c * x^2 + d * x^3
Parameters:
Estimate Std. Error t value
Pr(>|t|)
a
-20.5763 15.8218 -1.301 0.19651
b 24.6922 9.0471 2.729
0.00754 **
c 1.2064 1.3875 0.869
0.38674
d 8.6786 0.5421 16.011 <
2e-16 ***
---
Signif.
codes: 0 `***' 0.001 `**' 0.01 `*'
0.05 `.' 0.1 ` ' 1
Residual
standard error: 106 on 97 degrees of freedom
Correlation
of Parameter Estimates:
a b
c
b 8.778e-09
c -7.454e-01
-3.082e-09
d -2.804e-08
-9.166e-01 2.545e-08
返された結果に基づいた多項式回帰モデルを次に示す。
読者が同じ方法で上記の操作を行っても本稿と同じ結果が得られないのは、間違いではなく、人工データを作成する際に正規乱数を生成して用いたためである。
統計モデルを作成する際には、パラメータの数を少なくするのが良い。そこで、簡潔な3次式
を用いて同じデータについて非線形回帰分析を行ってみよう。
>(fm4=nls(y~a+b*x^3,start=c(a=1,b=1),trace=T))
20715546 : 1 1
1181990 : -10.32228 10.03462
Nonlinear regression model
model: y ~ a + b * x^3
data: parent.frame()
a
b
-10.32228 10.03462
residual
sum-of-squares: 1181990
返された結果のモデルは
となる。この異なる2つの多項式モデルによるAIC値を次に示す。
> AIC(fm3)
[1] 1232.551
> AIC(fm4)[1]
[1] 1236.752
AICの値には大きな差がないことから、両モデルによる当てはめの良さには明らかな差がないことが分かる。
視覚的に考察するため、両モデルによる予測値を図2に加えるコマンドとその結果を図3に示す。
>plot(x,y)
>lines(x,fitted(fm3),lty=1,col=2)
>lines(x,fitted(fm4),lty=2,col=4,lwd=2)
>legend(1,-700,c("一般式","簡潔式"), lty=1:2,col=c(2,4),lwd=2)
図3から分かるように、両モデルの予測値には差がほとんど見られない。
図3 人工データの実値と両モデルの予測値
初期値は大まかなものであるので、多少異なってもプログラムは作動する。初期値の設定に左右されるのは計算量である。
このような多項式による回帰分析は、線形回帰関数lmを次のように記述することで等価の回帰モデルを求めることができる。
lm(y~poly(x,3))
上記の書き式の中のpoly(x,3)は、3次多項式(Polynomial)のリンクに関する記述である。
3．その他の非線形回帰モデル
(1) 平滑化
関数lm、glm、nlsが対応できない、より複雑なデータにおける回帰モデルの説明のため、次の人工データを用いる。
>x1=seq(-10,10,0.1)
>y1=50*sin(x1)+x1^2+10*rnorm(length(x1),0,1)
>plot(x1,y1)
上記のコマンドで生成されたx1を横軸、y1を縦軸とした散布図を図４に示す。データy1を被説明変数、x1を説明変数とした非線形回帰モデル考える場合、散布図から思い付くのは多項式回帰モデルであろう。
次の3種類の多項式回帰モデルのコマンドとその回帰モデルによる予測値を散布図に加えるコマンドを示す。
>lm.p2<-fitted(lm(y1~poly(x1,2)))
>lines(x1,lm.p2,lty=4,col=4,lwd=2)
#2次多項式モデルの予測値
>lmp5<-fitted(lm(y1~poly(x1,5)))
>lines(x1,lmp5,lty=2,col=2,lwd=2)
#5次多項式モデルの予測値
>lmp8<-fitted(lm(y1~poly(x1,8)))
>lines(x1,lmp8)
#8次多項式モデルの予測値
>legend(-2.5,140,c("poly(x1,2)","poly(x1,5)","poly(x1,8)"),col=c(4,2,1),lty=c(4,2,1),lwd=2)
#凡例の表示
図4から分かるように、このデータにおいては、多項式によるモデルの当てはめの精度が良いとは言いがたい。
図4 3種類の多項式回帰モデルの予測値
このように複雑に変化する被説明変数を推定する方法として平滑化回帰の方法が提案されている。平滑化回帰を詳しく説明する紙面がないので、ここではRで提供している関数smooth.spline、ksmooth、supsmu、lowessを用いた平滑化された回帰曲線を図５に示す。
>par(mfrow=c(2,2),oma=c(2,2,2,2),mar=c(2,2,2,2))
#1つの画面に4コマのグラフを配置
>plot(x1,y1,main="関数smooth.splineによる平滑結果")
>lines(smooth.spline(x1,y1),col=2,lwd=2)
#smooth.splineによる平滑化
>plot(x1,y1,main="関数ksmoothによる平滑結果")
>lines(ksmooth(x1,y1),col=4,lwd=2)
#ksmoothによる平滑化
>plot(x1,y1,main="関数supsmuによる平滑結果")
>lines(supsmu(x1,y1),col=3,lwd=2)
#supsmuによる平滑化
>plot(x1,y1,main="関数lowessによ平滑結果")
>lines(lowess(x1,y1,f=0.1),col=5,lwd=2)
#lowessによる平滑化
図5 平滑化の例
パッケージfieldsの中の関数sregを用いると平滑化スプライン(smoothing spline )回帰を行うことができる。
(2) 加法モデル
加法回帰モデル(additive
regression model)は、次のような数式で表されるモデルである。
式の中の
は被説明変数、
は説明変数、
は説明変数の変換を行う関数、
は残差である。
線形回帰モデル、nlsによる非線形回帰モデルは加法モデルの特殊なケースと考えられる。
が全て線形関数の場合は、線形回帰モデルとなる。
一般化線形モデルを加法モデル化したとき、一般化加法モデル(GAM:
generalized additive model)と呼ぶ。
Rのパッケージmgcvに一般化加法モデルの関数gam がある。関数gamでは、それぞれの説明変数を、平滑化関数をリンクして用いることもできる。
データairqualityを用いてその使用法を例示する。ここでは、用いる全ての説明変数について平滑化スプラインをリンクすることにする。平滑化スプラインsmooth.splineの略記号はsである。次にその書き式を示す。加法回帰モデルの要約はsummaryで返される。
> library(mgcv); data(airquality)
>airq.gam<-gam(Ozone~s(Solar.R)+s(Wind)+s(Temp),data=airquality)
>summary(airq.gam)
Family:
gaussian
Link
function: identity
Formula:
Ozone ~
s(Solar.R) + s(Wind) + s(Temp)
Parametric
coefficients:
Estimate std. err. t ratio Pr(>|t|)
(Intercept) 42.099 1.663 25.32 < 2.22e-16
Approximate
significance of smooth terms:
edf
chi.sq
p-value
s(Solar.R) 2.76
13.7
0.0038243
s(Wind) 2.91 50.559 5.3379e-09
s(Temp) 3.833 55.901 3.4155e-09
R-sq.(adj)
= 0.723 Deviance explained = 74.7%
GCV score
= 338.9 Scale est. = 306.83 n = 111
説明変数を平滑化スプライン関数にリンクした加法モデルの自由度調整済みの重相関係数は0.73である。
同じの説明変数を用いた関数lmによる線形回帰モデルの自由度調整済みの重相関係数は0.6059である(先月号を参照)ので、上記の加法モデルが線形回帰モデルより当てはめが良いと言えよう。
(3) その他の方法
上記の説明以外にも多くの回帰方法が提案されている。その全てを詳細に説明する余裕がないので、局所回帰モデル、射影追跡回帰モデル、樹木モデル、ニューラルネットワークモデルについて簡潔に紹介する。
局所回帰モデル(local regression model )は、局所的にデータに回帰関数、あるいは回帰曲面を当てはめる1つの方法である。平滑化関数lowessは局所多項式を用いている。Rには局所回帰関数loessがパッケージstatsの中に含まれている。関数loessは局所多項式による回帰モデルである。
射影追跡回帰モデル(projection pursuit regression model)は、スタンフォード大学のFriedmanらによって1980年代に提案された、射影した説明変数を加法モデルに適用する方法である。RのパッケージMASSに射影追跡回帰関数pprがある。
樹木モデル(tree-based model)は、社会科学、機械学習の分野で提案された決定木(decision trees)モデルである。統計学で決定木に対して注目し始めたのは、1980年代からである。Rには樹木モデルのパッケージtree、rpart、mvpartがある。
ニューラルネットワークモデル(neural networks mode)は人工知能の分野で発展してきたニューロンをモデル化したものであり、多くの方法が提案されているが、回帰モデルに用いるのは階層型ニューラルネットワークである。Rには、階層型ニューラルネットワークのパッケージnnetがある。
これらの回帰モデルに関する参考書として下記の参考文献がある。
参考文献
[1] S-PLUSによるデータ解析、伊藤幹夫・その他 訳、シュプリンガーフェアラーク東京。
[2] Sと統計モデル−データ科学の新しい波−、柴田里程 訳、共立出版株式会社。
[3] 平滑化とノンパラメトリック回帰への招待、竹澤邦夫・大森宏 訳、（財）農林統計協会。
[4] みんなのためのノンパラメトリック回帰（上）（下）、竹澤邦夫 著、吉岡書店。