Rで行えるプロビット回帰に関するまとめ~最尤法、ベイズ~

今回はRでプロビット回帰を行う際に使用できる分析方法と、その注意点を上げておきます。

プロビット回帰とは
最尤法
ベイズ
注意点とコメント
関連記事

プロビット回帰とは

予測変数が0か1しかないとき、例えば倒産する(1)、倒産しない(0)や、購買する(1)、購買しない(0)等のケースでは、普通線形回帰分析は利用しません。まあ別に利用することも出来るけど、予測値に1以上の値や0未満の値が出て来てどう考えたら良いかわからなくなるからです。
そこで、ロジスティック回帰やプロビット回帰といったような手法を利用します。

最尤法

モデル

プロビットモデルは次のようなモデルです。
まず、1が出る確率 $P(Y_i=1|X)$ を次のように定義します。

$P(Y_i=1|X) =\int_{-\infty}^{\alpha+\beta X} \frac{1}{\sqrt{2 \pi}} exp(-\frac{z^2}{2}) dz$

要は正規分布の分布関数です。分布関数であれば値は $[0,1]$ の範囲に収まるので、これを利用します。

後は $P(Y_i|X)$ の $i=1,2,\cdots,n$ に関する同時分布を考えて、それを尤度とします。
解は解析的には得られないので、数値最適化を利用することになります。

Rによるコード

Rでは、基本パッケージにプロビット回帰を行う関数が用意されており、簡単に行うことが出来ます。

glm(y~X1+X2+X3,family=binomial(probit))

シミュレーション

それでは、サンプルサイズを500から100000まで動かして、係数比を見てやることにします。

n_max=100000
nseq=seq(500,n_max,500)
glmbeta1_seq=c()
glmbeta2_seq=c()
glmbeta3_seq=c()
glmbeta4_seq=c()
beta0=-0.5
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=beta0+beta1*X1+beta2*X2+beta3*X3
probabl=pnorm(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]~X1[1:i]+X2[1:i]+X3[1:i],family=binomial(probit))
	glmbeta1_seq[which(nseq==i)]=coef(res)[1]
	glmbeta2_seq[which(nseq==i)]=coef(res)[2]
	glmbeta3_seq[which(nseq==i)]=coef(res)[3]
	glmbeta4_seq[which(nseq==i)]=coef(res)[4]
}

f:id:bananarian:20190117113303p:plain
赤線が真の値です。
結構近い値が出力されています。

ただ、ロジスティック回帰もそうですが、サンプルサイズがそこそこ大きくならないと使えない感じがありますね。
特にプロビット回帰だと正規分布の裾がロジットよりも薄いので尚更です。

ベイズ

当然ベイズでも出来ます。Stanを使えばそれ用のコードがあるので結構簡単です。

model="
data{
	int<lower=1> N;
	int<lower=1> K;
	matrix[N,(K-1)] 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];
		for(k in 2:K){
			mu[n] = mu[n]+beta[k]*X[n,(k-1)];
		}
	}
}

model{
	Y ~bernoulli(Phi(mu));
}
"


data=list(N=length(y),K=4,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: c368d88db0144dcae18220ad714858ba.
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%
beta[1]    -2.18    0.02 0.25    -2.64    -2.36    -2.19    -2.01    -1.68
beta[2]     0.56    0.00 0.01     0.54     0.55     0.56     0.56     0.57
beta[3]     0.93    0.00 0.01     0.90     0.92     0.93     0.94     0.95
beta[4]     0.92    0.00 0.01     0.88     0.91     0.92     0.92     0.94
lp__    -3218.43    0.09 1.43 -3222.04 -3219.18 -3218.06 -3217.44 -3216.77
        n_eff Rhat
beta[1]   260    1
beta[2]   151    1
beta[3]   154    1
beta[4]   151    1
lp__      260    1

Samples were drawn using NUTS(diag_e) at Thu Jan 17 14:59:04 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).

比をとってやると、真の値の比とほぼ同じ値をとっていることがわかります。

> 0.3/0.5
[1] 0.6
> 0.56/0.93
[1] 0.6021505
> 0.56/0.92
[1] 0.6086957

注意点とコメント

・プロビット回帰はロジット回帰よりも分布の裾が軽いので、0か1がまれにしか出ないようなデータには当てはまりが悪くなります。
・多変量に拡張しない限りは、正直ロジットとプロビットの差は裾の部分にしかないので、単なる二項選択ならロジットを使うのがラクチンではある。
詳しくはこちら
人間の選択をモデル化するプロビット・ロジットモデルの違いと経済学的解釈法 - バナナでもわかる話