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+b*exp(c*x))こんな感じです。a、b、cがパラメタですので、これを最小二乗法で求めることになります。traceはパラメタの計算過程を返すかどうか決めるものです。結果を返すときは、TRUEを指定します。以上を引数に書き込むとこんな感じです。

 

nls(y~ a/(1+b*exp(c*x)),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回はお終いです。