読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第36回 非線形回帰分析

多変量解析

第36回は非線形回帰分析について書きます。

これまでは線形回帰分析について書いてきました。線形回帰分析とは、目的変数を説明変数の線形関係で表すことでした。一方、非線形回帰分析では、目的変数を説明変数の非線形関係で表します。

例えば、標準シグモイド関数はよく使われます。y=1/(1+e^-x)で定義されます。シグモイド関数を一般化したロジスティック関数、y=a/(1+be^cx)を使う回帰分析を、ロジスティック回帰分析と呼びます。 

x <- seq(-10,10,1)
y <- 1/(1+exp(-x))
plot(x,y)
f:id:teruaki-sugiura:20150729233926p:plain

次に、適当になんかの製品の普及率を表すようなデータを作成します。ロジャースの普及理論に合うようなデータですね。

nonlinear <- read.csv("nonlinear.csv",header=TRUE)
head(nonlinear)
  year       penetration
1 1990       0.001
2 1991       0.012
3 1992       0.034
attach(nonlinear)
plot(year,penetration)

f:id:teruaki-sugiura:20150729233937p:plain この2つの図を見てみると、明らかに曲線関係が似ていることがわかります。当たり前ですが、シグモイド関数にあうように普及率のデータを作りましたので・・・では、そんなこと気にせず、非線形回帰分析を行なっていきます。使う関数は以下です。 nls(formula,data,start,trace)   引数 formulaには、「関係式を具体的に書く」必要があります。例えば、y~ a/(1+bexp(cx))こんな感じです。a、b、cがパラメタですので、これを最小二乗法で求めることになります。traceはパラメタの計算過程を返すかどうか決めるものです。結果を返すときは、TRUEを指定します。以上を引数に書き込むとこんな感じです。 nls(y~ a/(1+bexp(cx)),start=c(a=1,b=1,c=-1),trace=TRUE)

nls<- nls(penetration~ a/(1+b*exp(c*1:26)),start=c(a=1,b=1,c=-1))
#yearを使うと数値が大きすぎて計算できないので26年分を1:26で代替する
Nonlinear regression model
model: penetration ~ a/(1 + b * exp(c * 1:26))
data: parent.frame()
       a         b         c
  0.9628   135.8741  -0.4771
residual sum-of-squares: 0.008387
Number of iterations to convergence: 14
Achieved convergence tolerance: 1.044e-06
 
plot(year,penetration)
lines(year,penetration)
lines(year,fitted(nls(penetration~ a/(1+b*exp(c*1:26)),start=c(a=1,b=1,c=1))),col=2,lty=2,lwd=2)
legend(1997,0.2,c("observed value","predicted value"),col=1:2,lty=1:5,lwd=1)

f:id:teruaki-sugiura:20150729233949p:plain

関数nlsでは、その関数式をリンクして、非線形回帰分析を行うことができる。用意されている関数(例えば、SSlogis)をリンクして回帰分析を行うなどができる。

(ここらへんは一般化線形モデルの知識が必要と思われます。近々、一般化線形モデルについて書く予定です。)

SSlogis<-nls(penetration~SSlogis(year, Asym, xmid, scal ),list( Asym = 1, xmid = 1, scal = 1 ))
 
SSlogis
Nonlinear regression model
model: penetration ~ SSlogis(year, Asym, xmid, scal)
data: list(Asym = 1, xmid = 1, scal = 1)
     Asym      xmid      scal
   0.9628 1999.2948    2.0960
residual sum-of-squares: 0.008387
Number of iterations to convergence: 0
Achieved convergence tolerance: 4.062e-06
  

関係式を書き表して求めた予測値とSSlogisで求めた予測値を比較する。よく近似していることがわかります。つまり、どちらでやってもほぼ同じ結果となります。

data.frame(fitted(nls),fitted(SSlogis))
   fitted.nls.    fitted.SSlogis.
1   0.01128479      0.01128477
2   0.01805500      0.01805497
3   0.02876417      0.02876413
4   0.04551934      0.04551929

以上で第36回はお終いです。

広告を非表示にする