今回はロジスティック回帰について、Rで出来ることを数式込みでまとめておきます。
目次
スポンサーリンク
ロジスティック回帰とは
予測変数が0か1しかないとき、例えば倒産する(1)、倒産しない(0)や、購買する(1)、購買しない(0)等のケースでは、普通線形回帰分析は利用しません。まあ別に利用することも出来るけど、予測値に1以上の値や0未満の値が出て来てどう考えたら良いかわからなくなるからです。
そこで、ロジスティック回帰やプロビット回帰といったような手法を利用します。
最尤法
モデル
ロジスティック回帰は次のようなモデルです。
まず、1が出る確率を次のように定義します。
要は、1が出るか0が出るかという確率の構造の中にXが関係していると考えます。
とデータを取ってきたとするとその尤度
は次のように得られるはずです。
この尤度を最大化するパラメータの値を数値最適化を利用して探索的に(又は解析的に)求めます。
Rによるコード
Rではglm関数を利用してロジスティック回帰を簡単に行うことが出来ます。
glm(y~X1+X2+X3,family=binomial(logit))
familyにlogitであることを指定してやればOKです。
又は、optim関数を使って探索的に解を探すことも可能です。
fr=function(beta,X1,X2,X3,Y){ sum(Y*(beta[1]*X1+beta[2]*X2+beta[3]*X3)-log(1+exp(beta[1]*X1+beta[2]*X2+beta[3]*X3))) } optim(c(0,0,0),fr,X1=X1,X2=X2,X3=X3,Y=y,method="BFGS",hessian=TRUE,control=list(fnscale=-1))
Penalized Likelihood (Firth)
数値的に 0 か 1 である確率が生じました
実際にglmを利用してロジスティック回帰をおこなっていると、時々次のようなエラーが出てきます。
数値的に 0 か 1 である確率が生じました
これ、何かというと、完全分離または準完全分離という状況が起こっている場合に出てくるエラーで、簡単に言ってしまうと、「現状のデータセットをモデルで完全に(又は境界を除いて完全に)予測出来てしまうようなケース」を指します。
この状況になると何が問題かというと、最尤推定を行った時に、値が一意にならず、探索する際の収束結果が安定しません。
つまり、推定結果がロバストではなくなります。
これは、ロジスティック回帰における小標本問題といって、まあもっとたくさん標本を取ってくれば解決するわけですが、これ以上標本を用意することはコスト的にも時間的にも厳しいという場合も多々あります。この場合にもいくつか改善策が提案されており、
・Exact logistic regression
・Penalized Likelihood logistic regression
などがあったりします。
※Exact logistic regressionを行うRパッケージは現状(2019年1月14日現在)無い。
一応elrmパッケージというものがあったっぽいけど、公式からremoveされているので、使わない方が良さげな気がする。
Penalized Likelihoodの方は、パッケージがあって、簡単に利用できるので、こっちを説明します。
この辺の詳しいpdfや論文は結構ネット上にあるので、リンクを貼っておきます。
https://www.jstor.org/stable/2336390
http://www.epistat.m.u-tokyo.ac.jp/admin/wp-content/uploads/2014/06/20140115tajima.pdf
http://www.applstat.gr.jp/jjas/V36N23_2.html
http://ir.c.chuo-u.ac.jp/repository/search/binary/p/4466/s/2278/
モデル
で、Penalized Likelihood logistic regressionですが、先ほどの最尤法で利用した尤度関数に対してフィッシャー情報行列を使って罰則を加えてやるモデルです。
フィッシャー情報行列を用いて
こうすることの利点は、まず完全分離においても安定した最尤推定値を出せるようになること、更に理論的にはExact logistic regressionに対する近似値を求めてることに等しく、ベイズの枠組みではジェッフェリーズ無情報事前分布を用いた推定と同値になります。
Rによるコード
Rでは、logistfパッケージに入っています。
library(logistf)
使い方はglmと同じでfirth=TRUEと指定してやればOKです。
logistf(y~0+X1+X2+X3,firth=TRUE)
両者のシミュレーションによる比較
データセット
とりあえず、サンプルサイズを1000から1000000まで1000刻みで動かして、その推定結果を見てみようと思う。
n_max=1000000 nseq=seq(1000,n_max,1000) glmbeta1_seq=c() glmbeta2_seq=c() glmbeta3_seq=c() firthbeta1_seq=c() firthbeta2_seq=c() firthbeta3_seq=c() beta1=0.3 beta2=0.5 beta3=0.5 X1=sample((-50):50,n_max,replace=TRUE) X2=sample(50:70,n_max,replace=TRUE) X3=sample((-70):(-50),n_max,replace=TRUE) mu=beta1*X1+beta2*X2+beta3*X3 probabl=exp(mu)/(1+exp(mu)) y=c() for(p in probabl){ if(p>=runif(1)){ y[which(probabl==p)]=1 } else{ y[which(probabl==p)]=0 } }
大標本での結果
とりあえず推定。
for(i in nseq){ res=glm(y[1:i]~0+X1[1:i]+X2[1:i]+X3[1:i],family=binomial(logit)) glmbeta1_seq[which(nseq==i)]=coef(res)[1] glmbeta2_seq[which(nseq==i)]=coef(res)[2] glmbeta3_seq[which(nseq==i)]=coef(res)[3] res2=logistf(y[1:i]~0+X1[1:i]+X2[1:i]+X3[1:i],firth=TRUE) firthbeta1_seq[which(nseq==i)]=coef(res2)[1] firthbeta2_seq[which(nseq==i)]=coef(res2)[2] firthbeta3_seq[which(nseq==i)]=coef(res2)[3] }
結果
ほとんど同じ動きをしているっぽい。
真の値と微妙にズレているけど、係数比で見ると真の値に近くなっていることが確認出来ます。
値が安定してくるのがサンプルサイズ20000とか30000なのも気になりますね。実際のデータセットはかなり大規模なものでないとかなり誤差が大きくなりそう。
小標本での結果
Firthの手法が有用なのは小標本のケースなので、最大サンプルサイズ200でも試してみました。
結果はこんな感じ。
Firthの方法は安定しているけど、普通の方法ではなんか変なとこに飛んでいってますね。
ベイズ
Stanコード(普通の無情報事前分布)
library(rstan) model=" data{ int<lower=1> N; int<lower=1> K; matrix[N,K] X; int<lower=0,upper=1> Y[N]; } parameters{ vector[K] beta; } transformed parameters{ vector[N] mu; for(n in 1:N){ mu[n]=beta[1]*X[n,1]; for(k in 2:K){ mu[n] = mu[n]+beta[k]*X[n,k]; } } } model{ Y ~bernoulli_logit(mu); } "
無情報事前分布を与えた際の実行と結果
データセット
サンプルサイズは10000で、さっきと同じようにデータを作ってみました。
n=10000 X1=sample((-50):50,n,replace=TRUE) X2=sample(50:70,n,replace=TRUE) X3=sample((-70):(-50),n,replace=TRUE) beta0=0 beta1=0.3 beta2=0.5 beta3=0.5 sigma=sqrt(3) mu=beta1*X1+beta2*X2+beta3*X3 probabl=exp(mu)/(1+exp(mu)) y=c() for(p in probabl){ if(p>=runif(1)){ y[which(probabl==p)]=1 } else{ y[which(probabl==p)]=0 } } data=list(N=length(y),K=3,X=cbind(X1,X2,X3),Y=y)
実行
fit=stan(model_code=model,data=data,seed=1234,pars="beta",chains=1)
結果
> fit Inference for Stan model: ab047d734a690c818ea61b55250e593c. 1 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=1000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] 0.22 0.00 0.01 0.21 0.22 0.22 0.22 0.23 285 1.01 beta[2] 0.37 0.00 0.01 0.35 0.37 0.37 0.38 0.40 233 1.01 beta[3] 0.37 0.00 0.01 0.35 0.37 0.37 0.38 0.40 230 1.01 lp__ -1468.07 0.07 1.09 -1470.91 -1468.59 -1467.81 -1467.25 -1466.76 265 1.00 Samples were drawn using NUTS(diag_e) at Mon Jan 14 14:09:11 2019. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
とりあえず、係数の大小関係は合っているので、オッズ比で考えれば問題なさそう。
小標本の時はpenalized likelihoodの理論的背景を元にジェッフェリーズ事前分布を与えた方が良いのかな？という感じですかね。
注意点
・ロジスティック回帰の係数を解釈する際は、パラメータの比でみる必要がある。つまり、例えばを見た場合、
の変化によって生じる選択確率への影響と、
の変化によって生じる選択確率への影響の比を見ていることになる。
・オッズ比による解釈が一般的。
・もし、係数を単独で解釈したいのであれば、全ての説明変数のスケールを揃えるために標準化してからパラメータを推定するか、次の線形モデルの形で解釈する必要がある。
コメント
・実際のモデリングにおいて、背後の生成過程はわからないので、確実に何かしらのノイズが発生するはず。具体的には次のような感じ。
切片項を用意しても、このノイズには対応できないので、その辺どうすればいいのかが気になるけど、結局オッズ比で解釈することになるし、係数値の正確さは大小関係さえあっていればそんなに気にする必要はない？うーん。場合に寄りけり？でも、係数値を使わなければいけない場面では、その辺りに対応できるよう同時分布を考えた方が良いのかもしれないとも思いました。
