2 値データの回帰モデルとして一般的に使われている

ロジスティック回帰

ロジット関数をリンク関数として構築した回帰モデルをロジット回帰あるいはロジスティック回帰などとよばれている。一般に特徴の有無、実験動物の生死などのような 2 つの値しか取り得ない 2 値データを解析する際に適用される。

例えば新しく発見された抗生物質 A が黄色ブドウ球菌の増殖を抑制する効果を持つことが期待されている。抗生物質 A の濃度によって効果がどのように変化するかをモデル化したい場合にロジスティック回帰が用いられる。

実験としてまず抗生物質 A を x mM を培地に含ませ、次にその培地に黄色ブドウ球菌の培養液を 10 箇所垂らす。1 日後に発生したコロニーの数を数えることによって、抗生物質 A の効果を推定できる。例えば、コロニーが 10 個確認されれば、その濃度における抗生物質 A は黄色ブドウ球菌の増殖をまったく抑制できなかったことを示す。一方で、コロニーが確認されていなければ、抗生物質 A が高い抑制効果を発揮したことを示す。

抗生物質 A の濃度 x が同じならばどのコロニーも同じ確率で増殖が抑制される。増殖が抑制される確率を p とおく。抗生物質の濃度が x すなわち抑制確率が p のとき、培地に垂らしたスポット数 m のうち、y 個のコロニーの増殖が抑制された時の確率分布は以下のように書ける。(上の例では、m = 10 である)

\[ P(Y = y) = \begin{pmatrix}m \\ y \end{pmatrix}p^{y}(1-p)^{m-y} \]

次に、抗生物質 A の濃度を変化させて x = x1, x2, ..., xn としたとき、その同時確率は以下のように表せる。

\[ P(Y_{1}=y_{i}, \cdots, Y_{n}=y_{n}) = \prod_{i=1}^{n}\begin{pmatrix}m_{i} \\ y_{i} \end{pmatrix}p_{i}^{y_{i}}(1-p_{i})^{m_{i}-y_{i}} \]

抗生物質 A の濃度を変化させた時、抑制されるコロニーの数よりも抑制確率のほうが重要と考えられる。そこで各濃度 xi における抑制確率を πi とすると、E(Yi)/mi = πi により、モデルは以下のように構築できる。

\[E[\mathbf{Y}] = \mathbf{\pi} \] \[g(\mathbf{\pi}) = \mathbf{X}\mathbf{\beta}\]

2 つ目の式について、抗生物質の濃度が変化することによって 1 よりも大きくなると考えられる。そこで、リンク関数 g の逆関数は、その右辺を 0 から 1 の範囲に収める関数でなければならない。そのような関数にはロジット関数やプロビット関数がある。ロジット関数を利用すると以下のように書ける。

\[g(\mathbf{\pi}) = \log\left(\frac{\mathbf{\pi}}{1-\mathbf{\pi}}\right) = \mathbf{X}\mathbf{\beta} = \begin{pmatrix}1 & x\end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2}\end{pmatrix}\]

ここまでで、誤差構造は二項分布、リンク関数はロジット関数と利用することを確認した。

次に、サンプルデータとして次のデータを利用する。x の濃度を 10 段階に設定し、それぞれの段階において培地に垂らしたスポット数を m、増殖が抑制されたコロニー数 y とする。

x <- 1:10
m <- c(10, 12, 10, 9, 12, 11, 10, 9, 9, 10)
y <- c(0,   0,  1, 1,  3,  6,  6, 7, 8, 10)

glm 関数を利用してモデルを構築する。上記のように、スポット数 m が一定で、増殖が抑制されたコロニー数 y がわかれば、抑制されていないコロニー数も計算できる。そのため、誤差構造が二項分布のとき、増殖が抑制されたコロニー数とそうでない数を同時に応答変数として与えることができる。

logit.model <- glm(cbind(y, m - y) ~ x, family = binomial("logit"))

summary 関数を利用してモデルの要約情報を示す。

summary(logit.model)
## Call:
## glm(formula = cbind(y, m - y) ~ x, family = binomial("logit"))
## 
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -0.7730  -0.3857  -0.2131   0.3931   0.8620
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -5.4191     1.0346  -5.238 1.63e-07 ***
## x             0.8693     0.1634   5.320 1.04e-07 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 67.5001  on 9  degrees of freedom
## Residual deviance:  2.6302  on 8  degrees of freedom
## AIC: 23.025
## 
## Number of Fisher Scoring iterations: 4

Coefficients の項を見ると Estimate Intercept すなわち β1 が -5.4191 であり、Estimate xx すなわち β2 が 0.8693 であることがわかる。

各濃度 x における抑制確率 π はモデル式を変更することにより以下のように求まる。

\[ \begin{eqnarray} && \log\left(\frac{\pi}{1-\pi}\right) = \beta{1} + \beta_{2}x \\ &\Leftrightarrow& \pi = \frac{\exp (\beta_{1} + \beta_{2}x)}{1 + \exp (\beta_{1} + \beta_{2}x)} \end{eqnarray} \]

この式を利用すれば実際のデータの散布図上にモデル化された曲線を描くことができる。

plot(x, y / m, xlim = range(x), ylim = c(0, 1),
     xlab = "dose [mg]", ylab = "crisis rate")
par(new = T)
curve(exp(-5.4191 + 0.8693 * x) / (1 + exp(-5.4191 + 0.8693 * x)),
      xlim = range(x), ylim = c(0, 1), xlab = "", ylab = "")

確率に変換してから解析する場合

上では観測したデータをそのまま利用した。しかし、解析する前に、予め増殖が抑制される確率を計算してから解析することは不可能でもない。ただし、観測値を確率に変換するという作業は情報のロスが生じうえ、従属変数が二項分布に従わなくなるため、解析結果に悪影響を及ぼす。

x <- 1:10
m <- c(10, 12, 10, 9, 12, 11, 10, 9, 9, 10)
y <- c(0,   0,  1, 1,  3,  6,  6, 7, 8, 10)
r <- y / m

tmp.model <- glm(r ~ x, family = binomial("logit"))
summary(tmp.model)