Stan/RStan を利用したロジスティック回帰 / 二項分布モデルのパラメーター推定

ロジスティック回帰

従属変数が 2 つの値をしか取らない 2 値データは一般にロジスティック回帰で解析を行う場合が多い。2 値データが観測できる実験として、マウスにある薬剤を投与する実験において、投与量 X に応じてマウスの死亡率がどのように変化するかを調べるなどが挙げられる。

ここで次に示すサンプルデータを利用して、ロジスティック回帰を行う例を示す。サンプルデータは、ある薬剤を 10 種類の濃度に調整し、それぞれの濃度 x の薬剤を 10 匹のマウスに投与し、その薬剤の影響で死亡したマウスの数を y とする。

x <- c(0.00, 0.10, 0.21, 0.29, 0.40, 0.49, 0.60, 0.71, 0.80, 0.89, 1.00)
y <- c(0, 0, 1, 1, 3, 6, 6, 7, 8, 10, 10)

GLM によるロジスティック回帰

R の glm 関数によってロジスティック回帰を行う場合は、以下のようにする。

m <- glm(cbind(y, 10 - y) ~ x, family = binomial("logit"))
m
## Call:  glm(formula = cbind(y, 10 - y) ~ x, family = binomial("logit"))
## 
## Coefficients:
## (Intercept)            x
##      -4.114        7.771
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:	    77.8
## Residual Deviance: 5.544 	AIC: 26.55

plot(x, y / 10, xlim = range(x), ylim = c(0, 1),
     xlab = "dose [mg]", ylab = "crisis rate")
par(new = T)
curve(exp(-4.114 + 7.771 * x) / (1 + exp(-4.114 + 7.771 * x)),
      xlim = range(x), ylim = c(0, 1), xlab = "", ylab = "")
R glm 関数によるロジスティック回帰の結果

Stan によるロジスティック回帰

Stan を利用したロジスティック回帰モデルのパラメーター推定では、観測データを二項分布のパラメーターの形に変換する必要がある。Stan コードでは、data ブロックには、入力値であるサンプル数 N、濃度 x および死亡数 y を定義する。次に、parameters ブロックで、ロジスティック回帰モデルに組み込むためのパラメーター beta_0 および beta_1 を定義する。パラメーター beta_0 および beta_1 は、そのまま二項分布のパラメーターとして利用できないので、transformed parameters で変数変換を行う。これらの準備が完了してから、最後に model ブロックで、「確率 p で成功する試行を N 回トライしたあとに、y 回成功した」という二項分布モデルを構築する。

なお、変数変換について、ロジスティック回帰では、確率 p とパラメーター β0, β1 の関係は次のようになっているから

\[ logit(p_{i}) = \log\left(\frac{p_{i}}{1-p_{i}}\right) = \beta_{0} + \beta_{1}x\]

確率 p は次のようにして、パラメーター β0, β1 から計算できる。

\[ p_{i} = logit^{-1}(\beta_{0} + \beta_{1}x) \]

このことに注意して、transformed parameters ブロックは次のように書ける。

data {
    int N;
    real x[N];
    int y[N];
}
parameters {
    real beta_0;
    real beta_1;
}
transformed parameters {
    real p[N];
    for (n in 1:N)
        p[n] = inv_logit(beta_0 + beta_1 * x[n]);
}
model {
    y ~ binomial(N, p);
}

R で次のようにして Stan コードを呼び出して、ベイズ推定を行う。GLM と同じ結果にならなかったものの、やや近い結果が得られたことがわかる。

library(rstan)
d <- list(x = x, y = y, N = length(x))
fit <- stan(file = 'logit.stan', data = d)
fit
## Inference for Stan model: logit.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## beta_0  -3.87    0.02 0.67  -5.28  -4.33  -3.85  -3.40  -2.64   796    1
## beta_1   6.66    0.04 1.11   4.58   5.89   6.63   7.39   8.90   768    1
## p[1]     0.02    0.00 0.02   0.01   0.01   0.02   0.03   0.07   765    1
## p[2]     0.04    0.00 0.02   0.01   0.03   0.04   0.06   0.10   794    1
## ...
## p[11]    0.93    0.00 0.03   0.85   0.92   0.94   0.96   0.98   999    1
## lp__   -51.73    0.03 0.97 -54.37 -52.12 -51.44 -51.03 -50.76  1388    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan 10 12:35:07 2018.
## 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).
R glm 関数によるロジスティック回帰の結果