正規線形モデルの仮説検定

この記事は RICORA Advent Calendar 2022 15日目の記事になります

はじめに

正規線形モデルの（パラメータの個数についての）仮説検定を考えます。仮説検定に使用する統計量は様々考えられますが、今回は逸脱度を用いて作成します。記事の前半では、仮説検定の方法を提示し統計量の標本分布を導出します。後半では、Pythonを用いて、人工的に作成したデータを使って検出力を計算します。

表記や前提

データ数: $N$
目的変数: $\mathbf{Y} = (Y_1, \ldots, Y_N)^T$
計画行列: $\mathbf{X} = (\mathbf{x}_1, \ldots, \mathbf{x}_N)^T$であり、$\mathbf{X}$の各列は線型独立とします。
標本: $(\mathbf{x}_1, y_1), \ldots (\mathbf{x}_N, y_N)$を同一母集団からの無作為標本とします。

理論

仮説検定

$\mathbf{Y} \sim N(\mathbf{X\beta}, \sigma^2\mathbf{I})$という正規線形モデルに関して、$q < p < N$として、以下の帰無仮説$H_0$と対立仮説$H_1$を考えます。

$$\begin{aligned} \begin{cases} H_0: \mathbf{\beta} = \mathbf{\beta}_0 = (\beta_1, \ldots, \beta_q)^T \\ H_1: \mathbf{\beta} = \mathbf{\beta}_1 = (\beta_1, \ldots, \beta_p)^T \end{cases} \end{aligned}$$

$H_0$に対応するモデル$M_0$と$H_1$に対応するモデル$M_1$を考えます。$M_0$の方が単純なモデルであり、より一般的なモデルである$M_1$の特別な場合となっています。もし、$M_0$が$M_1$と同じくらいデータに適合するのであれば、節約の観点から単純な方が選択されて、$H_0$が採択されます。逆に、$M_1$の方が有意にデータに適合するのであれば、$H_0$は棄却されます。この後、適合度合いを測る統計量の標本分布を導出します。

対数尤度比統計量

推定されうるパラメータの最大個数を含んだ飽和モデルと、考えたいモデルの尤度比を考えます。飽和モデルの計画行列を$N \times N$の行列$\mathbf{X}_{max}$、パラメータを$\mathbf{\beta}_{max} = (\beta_1, \ldots, \beta_N)^T$とします。また、$\mathbf{\beta}, \mathbf{\beta}_{max}$の最尤推定量を$\hat{\mathbf{\beta}}, \hat{\mathbf{\beta}}_{max}$とします。この時、尤度比は以下になります。

$$\begin{aligned} \lambda = \frac{L(\hat{\mathbf{\beta}}_{max}; \mathbf{y})}{L(\hat{\mathbf{\beta}}; \mathbf{y})} \end{aligned}$$

この後、尤度比の対数を使用します。

$$\begin{aligned} log\lambda = l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) - l(\hat{\mathbf{\beta}}; \mathbf{y}) \end{aligned}$$

逸脱度

逸脱度を次のように定義します。

$$\begin{aligned} D &= 2log\lambda \\ &= 2[l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) - l(\hat{\mathbf{\beta}}; \mathbf{y})] \end{aligned}$$

正規分布の対数尤度関数は、

$$\begin{aligned} l(\mathbf{\beta}; \mathbf{y}) = -\frac{1}{2\sigma^2}(\mathbf{y} - \mathbf{X}\mathbf{\beta})^T(\mathbf{y} - \mathbf{X}\mathbf{\beta}) - \frac{1}{2}Nlog(2\pi \sigma^2) \end{aligned}$$

と表されます。飽和モデルでは、$\mathbf{X}_{max}\hat{\mathbf{\beta}}_{max} = \mathbf{y}$となるので、対数尤度関数の最大値は

$$\begin{aligned} l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) = - \frac{1}{2}Nlog(2\pi \sigma^2) \end{aligned}$$

となります。パラメータ数が$p < N$であるような他のモデルに対して最尤推定量は、

$$\begin{aligned} \hat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y} \end{aligned}$$

であることより、対応する対数尤度関数の最大値は以下になります。

$$\begin{aligned} l(\hat{\mathbf{\beta}}; \mathbf{y}) = -\frac{1}{2\sigma^2}(\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}})^T(\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}}) - \frac{1}{2}Nlog(2\pi \sigma^2) \end{aligned}$$

これより、逸脱度は

$$\begin{aligned} D &= 2[l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) - l(\hat{\mathbf{\beta}}; \mathbf{y})] \\ &= \frac{1}{\sigma^2}(\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}})^T(\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}}) \\ &= \frac{1}{\sigma^2}[\mathbf{y} - \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}]^T[\mathbf{y} - \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}] \\ &= \frac{1}{\sigma^2}[(\mathbf{I} - \mathbf{H})\mathbf{y}]^T[(\mathbf{I} - \mathbf{H})\mathbf{y}] \quad (\mathbf{H} \equiv \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T) \\ &= \frac{1}{\sigma^2}\mathbf{y}^T(\mathbf{I} - \mathbf{H})\mathbf{y} \quad (\because 補足1) \end{aligned}$$

となります。

補足1
$\mathbf{I} - \mathbf{H}$は射影行列、つまり

$$\begin{aligned} (\mathbf{I} - \mathbf{H})^T &= \mathbf{I}^T - \mathbf{H}^T = \mathbf{I} - \mathbf{H} \\ (\mathbf{I} - \mathbf{H})^2 &= \mathbf{I} - 2\mathbf{H} + \mathbf{H}^2 = \mathbf{I} - \mathbf{H} \quad (\because \mathbf{H}^2 = \mathbf{H}) \end{aligned}$$

であるので、

$$\begin{aligned} (\mathbf{I} - \mathbf{H})^T(\mathbf{I} - \mathbf{H}) = (\mathbf{I} - \mathbf{H})^2 = \mathbf{I} - \mathbf{H} \end{aligned}$$

となります。

---

$\mathbf{Y} \sim N(\mathbf{X\beta}, \sigma^2\mathbf{I})$のとき、

$$\begin{aligned} (\mathbf{I} - \mathbf{H})\mathbf{y} &\sim N((\mathbf{I} - \mathbf{H})\mathbf{X}\mathbf{\beta}, \sigma^2(\mathbf{I} - \mathbf{H})^T(\mathbf{I} - \mathbf{H})) \\ &= N(\mathbf{0}, \sigma^2(\mathbf{I} - \mathbf{H})) \quad (\because (\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{0}, 補足1) \end{aligned}$$

補足2と補足3より、

$$\begin{aligned} &[(\mathbf{I} - \mathbf{H})\mathbf{y}]^T(\sigma^2(\mathbf{I} - \mathbf{H}))^-[(\mathbf{I} - \mathbf{H})\mathbf{y}] \sim \chi^2(N-p) \\ \Rightarrow \quad &\frac{1}{\sigma^2}\mathbf{y}^T(\mathbf{I} - \mathbf{H})\mathbf{y} \sim \chi^2(N-p) \\ \Rightarrow \quad &D \sim \chi^2(N-p) \end{aligned}$$

---

補足2
n次元の確率変数$\mathbf{y} = N(\mathbf{\mu}, \mathbf{V})$について、$\mathbf{V}$のrankが$k \lt n$で特異だとします。このとき、$\mathbf{V}$の一般逆行列を$\mathbf{V}^-$と表すと、確率変数$\mathbf{y}^T\mathbf{V}^-\mathbf{y}$は、自由度$k$、非心パラメータ$\mathbf{\mu}^T\mathbf{V}^-\mathbf{\mu}$の非心カイ二乗分布$\chi^2(k, \mathbf{\mu}^T\mathbf{V}^-\mathbf{\mu})$に従います。

補足3
射影行列$\mathbf{A}$について、$tr(\mathbf{A}) = rank(\mathbf{A})$が成り立ちます。いま、$\mathbf{I} - \mathbf{H}$は射影行列であるので、$tr(\mathbf{I}) = N, tr(\mathbf{H}) = tr((\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{X}) = p$より

$$\begin{aligned} rank(\mathbf{I} - \mathbf{H}) = tr(\mathbf{I} - \mathbf{H}) = tr(\mathbf{I}) - tr(\mathbf{H}) = N - p \end{aligned}$$

が成り立ちます。

---

$\mathbf{Y} \sim N(\mu, \sigma^2\mathbf{I}), \mu \ne \mathbf{X\beta}$のとき、

$$\begin{aligned} (\mathbf{I} - \mathbf{H})\mathbf{y} &= N((\mathbf{I} - \mathbf{H})\mu, \sigma^2(\mathbf{I} - \mathbf{H})) \quad (\because (\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{0}) \end{aligned}$$

補足2と補足3より、

$$\begin{aligned} &[(\mathbf{I} - \mathbf{H})\mathbf{y}]^T(\sigma^2(\mathbf{I} - \mathbf{H}))^-[(\mathbf{I} - \mathbf{H})\mathbf{y}] \sim \chi^2(N-p, [(\mathbf{I} - \mathbf{H})\mathbf{\mu}]^T(\sigma^2(\mathbf{I} - \mathbf{H}))^-[(\mathbf{I} - \mathbf{H})\mathbf{\mu}]) \\ \Rightarrow \quad &\frac{1}{\sigma^2}\mathbf{y}^T(\mathbf{I} - \mathbf{H})\mathbf{y} \sim \chi^2(N-p, \frac{1}{\sigma^2}\mathbf{\mu}^T(\mathbf{I} - \mathbf{H})\mathbf{\mu}) \\ \Rightarrow \quad &D \sim \chi^2(N-p, \frac{1}{\sigma^2}\mathbf{\mu}^T(\mathbf{I} - \mathbf{H})\mathbf{\mu}) \end{aligned}$$

---

適合度合いを図る統計量の標本分布

$M_0, M_1$に対応する計画行列を$\mathbf{X}_0, \mathbf{X}_1$、最尤推定量を$\hat{\mathbf{\beta}_0}, \hat{\mathbf{\beta}_0}$、逸脱度を$D_0, D_1$とします。このとき、

$$\begin{aligned} D_0 &= \frac{1}{\sigma^2}(\mathbf{y} - \mathbf{X}_0\hat{\mathbf{\beta}_0})^T(\mathbf{y} - \mathbf{X}_0\hat{\mathbf{\beta}_0}) \\ &= \frac{1}{\sigma^2}(\mathbf{y}^T\mathbf{y} - \mathbf{y}^T\mathbf{X}_0\hat{\mathbf{\beta}}_0 - \hat{\mathbf{\beta}}_0^T\mathbf{X}_0^T\mathbf{y} + \hat{\mathbf{\beta}}_0^T\mathbf{X}^T\mathbf{X}_0\hat{\mathbf{\beta}}_0) \\ &= \frac{1}{\sigma^2}(\mathbf{y}^T\mathbf{y} - 2\hat{\mathbf{\beta}}_0^T\mathbf{X}_0^T\mathbf{y} + \hat{\mathbf{\beta}}_0^T\mathbf{X}^T\mathbf{X}_0\hat{\mathbf{\beta}}_0) \\ &= \frac{1}{\sigma^2}(\mathbf{y}^T\mathbf{y} - \hat{\mathbf{\beta}}_0^T\mathbf{X}_0^T\mathbf{y}) \\ D_1 &= \frac{1}{\sigma^2}(\mathbf{y} - \mathbf{X}_1\hat{\mathbf{\beta}_1})^T(\mathbf{y} - \mathbf{X}_1\hat{\mathbf{\beta}_1}) \\ &= \frac{1}{\sigma^2}(\mathbf{y}^T\mathbf{y} - \hat{\mathbf{\beta}}_1^T\mathbf{X}_1^T\mathbf{y}) \end{aligned}$$

となります。、$M_1$はより一般的なモデルであったので、データによく適合すると仮定します。ゆえに、$D_1 \sim \chi^2(N - p)$となります。一方、$D_0$はデータによく適合する場合には、$D_0 \sim \chi^2(N - q)$となり、このとき補足4より$D_0 - D_1 \sim \chi^2(p-q)$です。適合しない場合には、$D_0 \sim \chi^2(N - q, \mu^T(\mathbf{I} - \mathbf{H})\mu$)となり、このとき補足4より$D_0 - D_1 \sim \chi^2(p-q, \mu^T(\mathbf{I} - \mathbf{H})\mu)$です。

補足4
Cochranの定理より、次の定理が成り立ちます。2つの確率変数$X_1^2 \sim \chi^2(m), X_2^2 \sim \chi^2(k), m \gt k$を考えます。このとき、$X^2 = X_1^2 - X_2^2 \ge 0$ならば、$X^2 \sim \chi^2(m - k)$が成り立ちます。

---

$M_0$は$M_1$の特別な場合であったことから、

$$\begin{aligned} & l(\hat{\mathbf{\beta}_1}; \mathbf{y}) - l(\hat{\mathbf{\beta}_0}; \mathbf{y}) \ge 0 \\ \Rightarrow \quad &2[l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) - l(\hat{\mathbf{\beta}_0}; \mathbf{y})] - 2[l(\hat{\mathbf{\beta}}_{max}; \mathbf{y}) - l(\hat{\mathbf{\beta}_1}; \mathbf{y})] \ge 0 \\ \Rightarrow \quad &D_0 - D_1 \ge 0 \end{aligned}$$

となります。また、$q \lt p$より$N-q \gt N-p$です。よって、$D_0$が（中心）カイ二乗分布に従うとき$D_0 - D_1 \sim \chi^2(p-q)$が成り立ちます。

---

$D_0$が非心カイ二乗分布に従うときも同様の結果が成り立ち、$D_0 - D_1 \sim \chi^2(p-q, \mu^T(\mathbf{I} - \mathbf{H})\mu)$となります。詳しくは、Rao(1973)の第3章をご覧ください。

以上より

$$\begin{aligned} F^* &= \frac{D_0 - D_1}{p - q} / \frac{D_1}{N - p} \\ &= \frac{\hat{\mathbf{\beta}}_1^T\mathbf{X}_1^T\mathbf{y} - \hat{\mathbf{\beta}}_0^T\mathbf{X}_0^T\mathbf{y}}{p-q} / \frac{\mathbf{y}^T\mathbf{y} - \hat{\mathbf{\beta}}_1^T\mathbf{X}_1^T\mathbf{y}}{N-p} \end{aligned}$$

は、$H_0$が正しい場合には（中心）F分布$F(p-q, N-p)$に従い、$H_0$が正しくない場合には非心F分布$F(p-q, N-p, \mu^T(\mathbf{I} - \mathbf{H})\mu)$に従います。

棄却域

非心F分布の値は、同じ自由度を持つ（中心）F分布で期待される値より大きくなります。よって、有意水準を$\alpha\%$とすると、$F(p-q, N-p)$の上側$\alpha\%$点より$F^*$が大きくなった時に$H_0$を棄却すれば良いです。

Pythonを用いて検出力を考える

import文

from matplotlib import pyplot as plt
import numpy as np
from scipy.stats import f

データ数とパラメータ数

$$\begin{aligned} \begin{cases} H_0: \mathbf{\beta} = \mathbf{\beta}_0 = (\beta_1, \beta_2)^T \\ H_1: \mathbf{\beta} = \mathbf{\beta}_1 = (\beta_1, \beta_2, \beta_3)^T \end{cases} \end{aligned}$$

N = 30 # データ数
q = 2 # M0のパラメータ数
p = 3 # M1のパラメータ数
ND = 100000 # 生成するデータセット数

モデル

$$\begin{aligned} \begin{cases} M_0: Y_i \underset{i.i.d.}{\sim} N(\mu_i, \sigma^2), \mu_i = \beta_1 + \beta_2x_i \\ M_1: Y_i \underset{i.i.d.}\sim N(\mu_i, \sigma^2), \mu_i = \beta_1 + \beta_2x_i + \beta_3x_i^2 \end{cases} \end{aligned}$$

データ

以下のようなデータを作成します。$\mu_i$が$x_i$の2次関数になっていることから、$H_1$が真です。

$$\begin{aligned} \begin{cases} x_i \underset{i.i.d.}{\sim} U(0, 1) \\ Y_i \underset{i.i.d.}{\sim} N(\mu_i, 0.25), \mu_i = 1 - x_i + 2x_i^2 \end{cases} \end{aligned}$$

ここで、$U(0, 1)$は一様分布です。例えば、以下のようなデータが生成されます。

描画コード

np.random.seed(0)
es = np.random.normal(0, 0.25, N)
xs = np.random.random(N)
ys = 1 - xs + 2*xs**2 + es
xs_line = np.linspace(0, 1, 1000)
true_ys = 2*xs_line**2 - xs_line + 1

plt.scatter(xs, ys)
plt.plot(xs_line, true_ys)

有意水準と棄却域

有意水準は$5\%$とします。

print(f"F({p-q}, {N-p})の上側5%点: {f(p-q, N-p).ppf(0.95)}")

F(1, 27)の上側5%点: 4.210008468359754

この値より$F^*$が大きい場合、$H_0$を棄却します。

検出力

cnt = 0 # 棄却域に入ったデータセットの数
F_list = [] # F^*の値を記録
f95 = f(p-q, N-p).ppf(0.95) # F(p-q, N-p)の上側5%点

for i in range(ND):
    np.random.seed(i)

    es = np.random.normal(0, 0.25, N)
    xs = np.random.random(N)
    ys = 1 - xs + 2*xs**2 + es

    X0 = np.array([np.array([1, xs[i]]) for i in range(N)])
    b0 = np.dot(np.dot(np.linalg.inv(np.dot(X0.T, X0)), X0.T), ys)
    X1 = np.array([np.array([1, xs[i], xs[i]**2]) for i in range(N)])
    b1 = np.dot(np.dot(np.linalg.inv(np.dot(X1.T, X1)), X1.T), ys)

    tmp0 = np.dot((ys - np.dot(X0, b0)).T, ys - np.dot(X0, b0))
    tmp1 = np.dot((ys - np.dot(X1, b1)).T, ys - np.dot(X1, b1))
    F_star = ((tmp0 - tmp1) / (p-q)) / (tmp1 / (N-p))

    F_list.append(F_star)
    cnt += F_star > f95

print(f"検出力: {100 * cnt / ND}%")

検出力: 83.348%

$$\begin{aligned} \hat{F} = \frac{D_0}{N - q} / \frac{D_1}{N - p} \end{aligned}$$

検定統計量を上記の$\hat{F}$と$F^*$にして、検出力を比べることを考えます。$F^*$は第一自由度が1、第二自由度が27です。$\hat{F}$は第一自由度が28、第二自由度が27となります。以下のグラフでは、横軸は$\nu$です。縦軸は$F(n_1, n_2)$の上側$5\%$点を、$F(n_1, n_2, \nu)$の累積分布関数の引数にしたものです。

描画コード

from matplotlib import pyplot as plt
import numpy as np
from scipy.stats import f, ncf

nus = np.linspace(0, 100, 1000)
a, b = [], []
for nu in nus:
    a.append(ncf(1, 27, nu).cdf(f(1, 27).ppf(0.95)))
    b.append(ncf(28, 27, nu).cdf(f(28, 27).ppf(0.95)))
plt.plot(nus, a, label='F*')
plt.plot(nus, b, label='$\hat{F}$')
plt.grid()
plt.legend()

これは、$F(n_1, n_2, \nu)$に従う確率変数が棄却域に入らない確率を表しているので、小さいほどいいです。青い線の$F^*$はオレンジ色の線の$\hat{F}$よりも、常に小さいことがグラフから見てとれます。よって、検出力を考えると、検定統計量として$F^*$が$\hat{F}$より優れていることがわかります。

参考文献

[1] Dobson, A. J and Barnett, A. G. (2018). An introduction to generalized linear models,
4th Edition. Chapman and Hall/CRC.
[2] Rao, C. R. (1973). Linear Statistical Inference and Its Applications, 2nd Edition. Wiley,
New York.