3.4 ロジスティック回帰モデル

いま，サンプルサイズ$n$のデータが与えられているとする． $\boldsymbol y = (y_1, \ldots, y_n)^{\top}$を目的変数，その他の変数をまとめた変数ベクトル$\boldsymbol x_i = (1, x_i^{(1)}, x_i^{(2)}, \ldots, x_i^{(m)})^\top$，$i=1,2,\ldots,n$とする．つまり，目的変数は$m+1$個と考える．

$$\begin{align} \boldsymbol y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \in \{0, 1\}^n, \hspace{3mm} X = \begin{bmatrix} 1 & x_1^{(1)} & x_1^{(2)} & \cdots & x_1^{(m)} \\ 1 & x_2^{(1)} & x_2^{(2)} & \cdots & x_2^{(m)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n^{(1)} & x_n^{(2)} & \cdots & x_n^{(m)} \\ \end{bmatrix} \in R^{n \times (m+1)} \end{align}$$

目的変数$y_i$は$0,1$どちらかの値のみを取る．すなわち$y_i \in \{0,1\}$である． この時ロジスティック回帰モデルでは次のようなモデルを仮定する．すなわち$y_i=1$となるような確率$\pi(\boldsymbol x_i) = \pi_i$として

$$\begin{align} \tag{3.1} \pi(\boldsymbol x_i) = \frac{\exp(\boldsymbol x_i^\top \boldsymbol \beta)}{1+\exp(\boldsymbol x_i^\top \boldsymbol \beta)} \end{align}$$

と定義する．ここで，$\boldsymbol\beta = (\beta_0, \beta_1, \ldots, \beta_m)^\top$であり，

$$\begin{align} \boldsymbol x^{\top} \boldsymbol\beta = \beta_0 + \beta_1 x_i + \cdots + \beta_m x_m \end{align}$$

である． 常に分母$>$分子が成り立つことと，指数関数の値は非負であることから，$0 \leq \pi(\boldsymbol x_i) \leq 1$が成り立つ．この式の右辺はロジスティック関数とも呼ばれ，図 3.1のようになる．

# プロットするためのベクトルを生成
f <- function(x){
  exp(0.5 + 1.0 * x) / {1 + exp(0.5 + 1.0 * x)}
}

x <- seq(-4, 4, 0.01)

dd <- data.frame(
  x = x,
  y = f(x) 
)

fig <- ggplot(dd) + 
  geom_line(aes(x=x, y=y)) +
  ggtitle("Logistic function") +
  ylab("p(x)")

plot(fig)

Figure 3.1: ロジスティック関数の形状

図の通り，$0~1$の範囲でS字カーブを描くような関数形をしている．特にこの関数はシグモイド関数とも呼ばれ

$$\begin{align} \frac{\exp(x)}{1 + \exp(x)} = \frac{1}{1 + \exp(-x)} \end{align}$$

と変形できることから後者の表現で定義されている場合もある．念のためこれが同じものであることを確かめておこう． $\exp(-x) = [\exp(x)]^{-1}$であることに注意すると

$$\begin{align} \frac{\exp(x)}{1 + \exp(x)} &= \left\{\left[ \frac{\exp(x)}{1 + \exp(x)} \right]^{-1} \right\}^{-1} \\ &= \left\{ \frac{\exp(-x)}{[1 + \exp(x)]^{-1}} \right\}^{-1} \\ &= \frac{1}{\exp(-x) (1 + \exp(x))} \\ &= \frac{1}{1 + \exp(-x)} \end{align}$$ とできる．

また，$\pi(\boldsymbol x_i)$に対して以下のように定義される変換をロジット変換と呼ぶ．

$$\begin{align} \tag{3.2} g(\boldsymbol x^{\top} \boldsymbol\beta) &= \log \left[ \frac{\pi(\boldsymbol x^{\top} \boldsymbol\beta)}{1 - \pi(\boldsymbol x^{\top} \boldsymbol\beta)} \right] \\ &= \boldsymbol x^{\top} \boldsymbol\beta \end{align}$$

ロジット変換について，$g(\boldsymbol x^{\top}) = \boldsymbol x^{\top} \boldsymbol\beta$となることは課題とする．

このモデルにおいては，それぞれのサンプルについて$y_i = 1$となる確率を観測することができないため，残差を求められず，最小二乗法を適用することができない． そのため，対数尤度と呼ばれる量を最大化する最尤法を用いる．

最尤法の手続きについて簡単に紹介する．回帰分析と同様に何らかの意味で最も良いパラメーターの推定値$\boldsymbol{\hat \beta}$を求めたい．このとき最尤法では，尤度関数という概念を用いて推定値を決定する．

尤度関数とは，確率関数$P(x;\theta)$または確率密度関数$f(x;\theta)$を$x$の関数ではなく，パラメーター$\theta$の関数と見たものである． 関数の形自体は変わらないが，変数の見方を変える．また$x$に関しては与えられたデータの値で固定しておく． このようにして$\theta$についての関数として最大化を検討する．

まず観測値が互いに独立であるという仮定の下でロジスティック回帰モデルの尤度関数$l(\boldsymbol\beta)$は

$$l(\boldsymbol \beta) = \prod_{i=1}^{n}\pi(\boldsymbol x_i^{\top} \boldsymbol\beta)^{y_i}[1-\pi(\boldsymbol x_i^{\top} \boldsymbol\beta)]^{1-y_i}$$

と表される．この関数の対数とったもの$L(\boldsymbol\beta)$を対数尤度関数とよぶ．すると$L(\boldsymbol\beta)$は

$$\begin{align} L(\boldsymbol\beta) &= \log(l(\boldsymbol\beta)) \\ &= \sum_{i=1}^{n} \left\{ y_i \log[\pi(\boldsymbol x_i^{\top} \boldsymbol\beta)] + (1-y_i) \log[1 - \pi(\boldsymbol x_i^{\top} \boldsymbol\beta)] \right\} \end{align}$$

となる．この関数を最大化するような$\boldsymbol\beta$を求めるために，$\boldsymbol\beta$で微分することで尤度方程式

$$\sum_{i=1}^{n}\boldsymbol x_i[y_i - \pi(\boldsymbol x_i^{\top} \boldsymbol\beta)] = \boldsymbol 0$$

を得る．この尤度方程式は非線形のため解析に解くことが難しいことが知られており，プログラムを用いて数値的に近似解を求めるのが一般的である．RやPythonでは多くの場合パラメーターを推定する関数などが開発者や第三者によって提供されているので分析者が個別に推定のためのプログラムを一から開発する必要はない．

大まかにはこのような考え方であるということを把握しておいてほしい．