時間ができたので一般化線形モデル(ロジスティック回帰)のコードを書いてみた
ロジスティック回帰
ロジスティック回帰モデルといえば、一般化線形モデルの代表格で、統計を少しでも勉強している人は知っている・聞いたことのある手法だと思います。今回は、そのロジスティック回帰を実際にRでプログラムするということをやってみます。ちなみに、Rはもともとglmという関数が入っていて、そこでfamily="binomial"とすれば、それで終わりなんですけど、まぁつくってみる・手を動かすのも大事なのでやってみました。
注)逸脱度の計算と、尤度計算はまだしていませんがちゃんとしないと...
一般化線形モデルの枠組みについて
一般化線形モデルの枠組みについては、Dobson(2002)がわかりやすいです。日本語も出ていて、「一般化線形モデル入門」でググってもらえばすぐにでます。
プログラムの内容
データはDobson(2002)のものをサンプルとして使用します。xは密室内に充満するガスの濃度、nはガスを注入した密室内にいる虫の数、yはそのうちガスの注入で死んだ虫の数です。このとき、死んだ虫の数yは2項分布に従うと考えることができて、これをガスの濃度で説明するモデルを考えます。これは、一般化線形モデルに置けるロジスティック回帰モデルです。ロジスティック回帰モデルというと、目的変数が「0か1」と考える方が多くいらっしゃいますが、大きな枠組みとしては、目的変数の分布に2項分布を仮定することです。つまり、正確な記述としては、
というものになります。ロジスティック回帰の名前の由来というのは、以下の式の左辺を「ロジット変換」と呼ぶからです。それから、この式の左辺を「リンク関数」と呼びます。そして、この形は「正準リンク」と呼ばれていて、このリンク関数を用いると嬉しい性質がいくつかありますが、今回は言及を避けましょう。ただ、皆さんが用いているロジスティック回帰モデルには、色々と背景があるんだなということを知っておいていただければ。
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で使う方法を用いる方がいいのかねぇ。。。