<<目次へ

ロジスティック回帰分析

　Rでロジスティック回帰分析を行うには最初から実装されているnls()という関数を使用します。この場合,パラメータを推定するための初期値はユーザが試行錯誤しながら考えなければなりません。そのような手間を省くための方法を紹介します。

　まず手始めに群馬大学の青木先生のwebサイトに置いてあるシンプレックス法のソースコードを拝借してください。このソースを貼り付けるなり,読み込むなりしておくことを前提として以降の説明をします。

　ロジスティック曲線を描く関数は以下のように定義されます。このときα,βが推定するパラメータです。なお、次式において分子は1になっているが、これをγとして、γの値を推定することもできる（比率データのように上限値が1と定まっている場合は、γとするより1とするべきである）。

手順1　以下のソースコードを貼り付けてください（読み込んでも良いですが）。

ln.dat<-function(dat,gamma){
        n<-length(dat[,1])
        Y<-c()
        x.data<-1:n
        for(i in 1:n){
                Y[i]<-log(gamma/dat[i,2])
        }
        result<-lm(Y~x.data)
        kaiki<-summary(result)
        return(coefficients(kaiki))
}

K.logistic<-function(dat,gamma,alpha,beta){
        x<-dat[,1]
        y<-dat[,2]
        my.model<-function(x, p){
                return(p[1]/(1+exp(p[2]+p[3]*x)))
        }
        simplex(my.model,c(gamma,alpha,beta),x,y,plot.flag=TRUE)
}

手順2　ln.dat()を使ってαとβの初期値を見つける

　この関数の引数であるdatはデータ行列,gammaは推定するパラメータγの初期値です。出力される切片（Intercept）をαの初期値として,また回帰係数（x.data）をβの初期値と考えてよいでしょう。

> my.dat　#1列目はx(説明変数),2列目はy(=目的変数)
      [,1]       [,2]
[1,] 1.691 0.06779661
[2,] 1.724 0.16666667
[3,] 1.755 0.30645161
[4,] 1.784 0.55357143
[5,] 1.811 0.82539683
[6,] 1.837 0.89830508
[7,] 1.861 0.96774194
[8,] 1.884 1.00000000
> ln.dat(my.dat,1)　#γ=1と仮定する
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  2.4982329  0.3193529  7.822796 0.0002303444
x.data      -0.3721354  0.0632413 -5.884373 0.0010679569

手順3　K.logistic()を使ってパラメータを推定する

　この関数の引数であるdatはデータ行列（上で用いたのと同じもの）,gammaはγの初期値,alphaはαの初期値,betaはβの初期値です。図は自動的に描かれます。

> K.logistic(my.dat, 1, 2.4982329, -0.3721354) #γ=1,α=2.4982329,β=-0.3721354
$converge
[1] TRUE

$parameters　#推定されたパラメータ
[1]   1.027154  61.098532 -34.381595

$residuals
[1] 11.57737

　以上のものは完全なものではありません。場合によってはこの方法でもうまくいかないことがあります。