時間ができたので一般化線形モデル（ロジスティック回帰）のコードを書いてみた

ロジスティック回帰

ロジスティック回帰モデルといえば、一般化線形モデルの代表格で、統計を少しでも勉強している人は知っている・聞いたことのある手法だと思います。今回は、そのロジスティック回帰を実際にRでプログラムするということをやってみます。ちなみに、Rはもともとglmという関数が入っていて、そこでfamily="binomial"とすれば、それで終わりなんですけど、まぁつくってみる・手を動かすのも大事なのでやってみました。
注）逸脱度の計算と、尤度計算はまだしていませんがちゃんとしないと...

一般化線形モデルの枠組みについて

一般化線形モデルの枠組みについては、Dobson(2002)がわかりやすいです。日本語も出ていて、「一般化線形モデル入門」でググってもらえばすぐにでます。

プログラムの内容

データはDobson(2002)のものをサンプルとして使用します。xは密室内に充満するガスの濃度、nはガスを注入した密室内にいる虫の数、yはそのうちガスの注入で死んだ虫の数です。このとき、死んだ虫の数yは2項分布に従うと考えることができて、これをガスの濃度で説明するモデルを考えます。これは、一般化線形モデルに置けるロジスティック回帰モデルです。ロジスティック回帰モデルというと、目的変数が「0か1」と考える方が多くいらっしゃいますが、大きな枠組みとしては、目的変数の分布に2項分布を仮定することです。つまり、正確な記述としては、

$Y_{i} \sim binomial(n_{i}, p_{i})\\ E(Y_{i}) = \mu_{i} = n_{i}p_{i}\\ \log(\frac{p_{i}}{1-p_{i}}) = x_{i}^{T}\beta = \eta_{i}$

というものになります。ロジスティック回帰の名前の由来というのは、以下の式の左辺を「ロジット変換」と呼ぶからです。それから、この式の左辺を「リンク関数」と呼びます。そして、この形は「正準リンク」と呼ばれていて、このリンク関数を用いると嬉しい性質がいくつかありますが、今回は言及を避けましょう。ただ、皆さんが用いているロジスティック回帰モデルには、色々と背景があるんだなということを知っておいていただければ。

$\log(\frac{p_{i}}{1-p_{i}}) = x_{i}^{T}\beta = \eta_{i}$

Rで書いた関数

#sample data
x=c(1.69,1.72,1.75,1.78,1.81,1.83,1.86,1.88)
n=c(59,60,62,56,63,59,62,60)
y=c(6,13,18,28,50,53,58,59)


#GLM for binomial

p_calc=function(x,b){
	exp=exp(x%*%b)
	pp=exp/(1+exp)
	return(pp)
}

glm_bin=function(x,n,y,eps=1e-08,fig=T){
	#data frame
	#if(is.vector(x)==TRUE){
		x=cbind(rep(1,length(x)),x) #intercept
	#}else{
	#	x=cbind(rep(1,nrow(x)),x) #intercept
	#}
	#start value
	beta=rep(0,ncol(x)) #regression coef
	pp=p_calc(x,beta) #probability
	mu=n*pp #expectation of y
	W = diag(as.numeric(n*pp*(1-pp))) #inverse of var(y)
	z = x%*%beta + solve(W)%*%(y-mu)
	#calc
	for(i in 1:25){
		stop = i #save i
		bbsave=beta
		beta = solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%z
		pp=p_calc(x,beta)
		mu=n*pp
		W = diag(as.numeric(n*pp*(1-pp)))
		z = x%*%beta + solve(W)%*%(y-mu)
		sqdf=sum(t(bbsave-beta)%*%(bbsave-beta))
		if(sqdf < eps){
			break;
		}
	}
	dev=2*(sum(y*log(y/(mu))+(n-y)*log(n-y)-(n-y)*log(n-mu)))
lp=x%*%beta
	#names(beta)=c("intercept","X1")
	ans=list(iter=i,beta=beta,prob=pp,linear.predictor=lp,y.fitted=mu,deviance=dev)
	return(ans)
}

#plot
#library
library(ggplot2)

#dataframe 
dat = data.frame(gas=x,rate=y/n)
gas=seq(min(x)-0.1,max(x)+0.1,0.001)
sam=cbind(rep(1,length(gas)),gas)
pline=sam%*%rr$beta
pline=exp(pline)/(1+exp(pline))
dat2 = data.frame(gas,pline)

#plot
(fig=ggplot(data=dat,aes(x=gas,y=rate))+geom_point()+geom_line(data=dat2,aes(x=gas,y=pline)))

結果の確認

#自作のプログラムの結果
$iter
[1] 6

$beta
       [,1]
  -60.81594
x  34.40013

$prob
          [,1]
[1,] 0.0641806
[2,] 0.1614176
[3,] 0.3507569
[4,] 0.6025951
[5,] 0.8097357
[6,] 0.8943813
[7,] 0.9596238
[8,] 0.9792919

$linear.predictor
           [,1]
[1,] -2.6797215
[2,] -1.6477177
[3,] -0.6157139
[4,]  0.4162900
[5,]  1.4482938
[6,]  2.1362964
[7,]  3.1683003
[8,]  3.8563028

$y.fitted
          [,1]
[1,]  3.786655
[2,]  9.685059
[3,] 21.746927
[4,] 33.745328
[5,] 51.013350
[6,] 52.768495
[7,] 59.496674
[8,] 58.757512

$deviance
[1] 1551.011

#Rに元々入っているものの結果
> rr2=glm(cbind(y,n-y)~x,data=dd,family="binomial",epsilon=1e-08)
> rr2

Call:  glm(formula = cbind(y, n - y) ~ x, family = "binomial", data = dd, 
    epsilon = 1e-08)

Coefficients:
(Intercept)            x  
     -60.82        34.40  

Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
Null Deviance:	    284.2 
Residual Deviance: 9.699 	AIC: 39.9

Coefficientが一致しているので問題ない。
ただ、iterationは通常のglm関数の方が少ないので、なんかもう少し関数考えた方がいいかもしれない．．．wood(2006)のGAMで使う方法を用いる方がいいのかねぇ。。。