廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Missing Not At Random(MNAR):R と Stan で欠測が欠測データに依存する場合のパラメータ推定

測定機器かなにかの都合上、観察対象の値が小さくなると欠損が出やすくなる状況を考えます。

R で以下のようにしてシミュレーションデータを生成しました。

set.seed(1)
N <-200
X<-rnorm(N,2)
X2 <-ifelse(runif(N)<plogis(1-2*X),NA,X)

プロットしてみます。

h1 <-hist(X, col = "#ff00ff40", border = "#ff00ff",main="")
hist(X2, col = "#0000ff40", border = "#0000ff", breaks = h1$breaks, add = TRUE)
par(family="Osaka")
legend("topright",c("完全データ","観測データ"),fill=c("#ff00ff","#0000ff"))

f:id:abrahamcow:20161123102901p:plain

欠損データを除外して平均や分散を求めるとバイアスが生じるであろうことは想像できます。

ではどうするか。

d_i を観測データ i が欠損のとき 1, そうでなければ 0 の値をとる変数として, 欠損のしやすさのモデルとして以下を考えることにしました。

 \Pr(d_i=1)=\rm{logit}^{-1}(\beta_0-\beta_1 z_{i})

ここで z_i は観測されなかった変数も含めてデータを拡大して表したものです。

z_i がそれぞれ独立に正規分布に従うとして, 完全なデータが得られたときの尤度関数は、

\prod_i \phi (z_i|\mu,\sigma) (\rm{logit}^{-1}(\beta_0+\beta_1 z_{i}))^{d_i} (1-\rm{logit}^{-1}(\beta_0+\beta_1 z_{i}))^{1-d_i}

です。φ正規分布の密度関数です。

Stan では完全なデータが得られた状況での尤度がわかっていればサンプリングが可能です。

欠損値の真の値 z_i はすべて未知パラメータとして扱います。

コーディングすると以下のようになりました。

# missing_normal.stan
data{
  int<lower=0> N;
  int<lower=0,upper=N> m;
  real z[N-m];
}
parameters{
  real mu;
  vector[m] z_miss;
  real<lower=0> sigma;
  real<lower=0> beta[2];
}
model{
  for(i in 1:(N-m)){
    target += normal_lpdf(z[i]|mu,sigma) +
      bernoulli_logit_lpmf(0|beta[1]-beta[2]*z[i]);
  }
  for(i in 1:m){
   target += normal_lpdf(z_miss[i]|mu,sigma) + 
     bernoulli_logit_lpmf(1|beta[1]-beta[2]*z_miss[i]); 
  }
  beta ~ exponential(0.1);
}

観察対象の値が小さくなると欠損が出やすくなることは事前に知っていたので、beta には正の値しか取らない指数分布を仮定することにしました。

サンプリングを実行します。

dat <- list(N=N,m=sum(is.na(X2)),z=X2[!is.na(X2)])
fit_normal_N200 <-stan("missing_normal.stan",data = dat,iter=10000)

MCMC が収束しているかはトレースプロットと Rhat 値で確かめました。

traceplot(fit_normal_N200,pars=c("mu","sigma","beta","z_miss[1]","z_miss[2]"))


f:id:abrahamcow:20161123104401p:plain

単純に欠損値を除いて平均と標準偏差を推定した場合(これを簡便法と呼ぶことにします)と、事後期待値を比べると、バイアスが改善していることがわかります。

> round(mean(X2,na.rm = TRUE),2)
[1] 2.2
> round(sd(X2,na.rm = TRUE),2)
[1] 0.92
> get_posterior_mean(fit_normal_N200,pars=c("mu","sigma","beta"))
        mean-chain:1 mean-chain:2 mean-chain:3 mean-chain:4 mean-all chains
mu          1.954692     1.952811     1.955300     1.954167        1.954243
sigma       1.030481     1.033981     1.030133     1.033480        1.032019
beta[1]     2.028536     2.134627     2.024334     2.066202        2.063425
beta[2]     4.411043     4.758496     4.365761     4.402222        4.484380

ただし beta の推定値に関してはまだまだです(真値は 1 と 2)。

サンプルサイズを増やして計算したところ推定値は以下のようになりました。

サンプルサイズ 簡便法 mu 簡便法 sigma Stan mu Stan sigma Stan beta[1] Stan beta[2]
500 2.22 0.92 2.00 1.04 1.66 2.55
1000 2.20 0.92 2.00 1.02 0.95 1.89

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)