🐙

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

2022/12/14に公開約12,500字

はじめに

正規線形モデルの(パラメータの個数についての)仮説検定を考えます。仮説検定に使用する統計量は様々考えられますが、今回は逸脱度を用いて作成します。記事の前半では、仮説検定の方法を提示し統計量の標本分布を導出します。後半では、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_0H_1に対応するモデルM_1を考えます。M_0の方が単純なモデルであり、より一般的なモデルであるM_1の特別な場合となっています。もし、M_0M_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}

となります。


\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}


\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)です。

以上より

\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_ix_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)

data

有意水準と棄却域

有意水準は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%

参考文献

[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.

Discussion

ログインするとコメントできます