対応のあるt検定を線形モデルのパラメータについてのF検定から導出する

## 記事の概要
対応のある$t$検定は、対応する観測同士の差を1つの観測値として扱い、「対応する観測同士の差」の分布の母平均が$0$であるか否かを検定するような方法である。素朴な方法であるが、これは本来$2n$個あった観測を$n$個の観測として扱うような方法であり、個々の観測が持つ誤差の取り扱い方が自明でない。この記事では線形モデルのパラメータについてのF検定から導く。
一般的な対応のある$t$検定の検定統計量は対応のある観測に対して固定効果を設定する固定効果モデルから導かれる。そのため、観測誤差は正規分布であることを仮定するが、個々の観測に紐づいた固定効果については正規分布を仮定しない。

## 問題設定
$n$個の$2$つの観測値が対応するような観測を
```math
(y^{(1)}_1, y^{(2)}_1),(y^{(1)}_2, y^{(2)}_2),...,(y^{(1)}_n, y^{(2)}_n)
```
とする。観測は真の値に平均$0$の分散が等しい正規分布に従う誤差を加えた値であると仮定して、
```math
\displaylines{
\begin{align}
y^{(1)}_{i} &= \beta_{i} + \varepsilon^{(1)}_{i} \\
y^{(2)}_{i} &= \beta_{i} + \delta + \varepsilon^{(2)}_{i}
\end{align} \\
\varepsilon^{(1)}_{1},...,\varepsilon^{(1)}_{n},\varepsilon^{(2)}_{1},...,\varepsilon^{(2)}_{n} \sim N(0, \sigma^2) \quad i.i.d.
}
```
というモデルを設定する。$\beta_{1}, \beta_{2},...,\beta_{n}$は$n$個の対応のある測定の各1つ目の測定に対応する真の値であり、これはそれぞれ独立した未知の定数で分布の仮定はしない。$\delta$は1つ目の測定の真の値と2つ目の測定の真の値の差であり、$n$個の対応のある測定で共通と仮定する。$\sigma^2$は測定誤差の共通の分散。帰無仮説$\delta=0$を検定するのが今回の目的となる。
このモデルを行列で表現するために、
```math
\boldsymbol{y}_1 = 
\begin{bmatrix}
y^{(1)}_1 \\
\vdots \\
y^{(1)}_n
\end{bmatrix}, \quad
\boldsymbol{y}_2 = 
\begin{bmatrix}
y^{(2)}_1 \\
\vdots \\
y^{(2)}_n
\end{bmatrix}, \quad
\boldsymbol{y} = 
\begin{bmatrix}
\boldsymbol{y}_1 \\
\boldsymbol{y}_2
\end{bmatrix}, \quad
\boldsymbol{X} = 
\begin{bmatrix}
\boldsymbol{0}_{n} & I_n \\
\boldsymbol{1}_{n} & I_n
\end{bmatrix}, \quad
\boldsymbol{\beta} = 
\begin{bmatrix}
\delta \\
\beta_1 \\
\vdots \\
\beta_n
\end{bmatrix}
```
とすれば、モデルは、
```math
\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \quad N(\boldsymbol{0}_{2n}, \sigma^2 I_{2n})
```
と表現できる。ただし、$\boldsymbol{0}_n, \boldsymbol{1}_n$はすべての要素が$0, 1$の$n$次元の列ベクトル。
```math
\boldsymbol{H} = \begin{bmatrix}
1 \\ \boldsymbol{0}_n
\end{bmatrix}
```
とすれば、帰無仮説$\delta = 0$の検定は、$\boldsymbol{\beta}$についての線形制約$\boldsymbol{H}^T \boldsymbol{\beta} = 0$の成立を帰無仮説とする検定と同値である。

## 具体的な導出
一般論として$\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}$の$\boldsymbol{\beta}$についての線形制約
```math
\boldsymbol{H}^T \boldsymbol{\beta}=\boldsymbol{d}, \quad \boldsymbol{H} \in \mathbb{R}^{p \times r}, \quad \boldsymbol{\beta} \in \mathbb{R}^{p}, \quad \boldsymbol{d} \in \mathbb{R}^{r}
```
の成立を帰無仮説とする$F$検定の検定統計量は、
```math
F = \frac{(\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d})^T (\boldsymbol{H}^T (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{H})^{-1} (\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d}) / r}{\boldsymbol{y}^T (I_{n} - \boldsymbol{X}(\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T) \boldsymbol{y} / (n-p)}
```
である。ただし、$\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{y}$で、$\boldsymbol{\beta}$の最小二乗解。
帰無仮説$\boldsymbol{H}^T \boldsymbol{\beta}=\boldsymbol{d}$が真ならば、この$F$は自由度$(r, n-p)$の$F$分布に従う。$\boldsymbol{H}^T \boldsymbol{\beta}=\boldsymbol{d}$について有意水準$\alpha$で検定する場合、帰無仮説棄却の条件は、
```math
\frac{(\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d})^T (\boldsymbol{H}^T (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{H})^{-1} (\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d})}{\boldsymbol{y}^T (I_{n} - \boldsymbol{X}(\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T) \boldsymbol{y}} \ge F^{(r, n-p)}_{1-\alpha}
```
である。ただし、$F^{(r, n-p)}_{1-\alpha}$は自由度$(r, n-p)$の$F$分布の上側$\alpha$の分位点を表す。

これを今回の線形モデル
```math
\boldsymbol{y}_1 = 
\begin{bmatrix}
y^{(1)}_1 \\
\vdots \\
y^{(1)}_n
\end{bmatrix}, \quad
\boldsymbol{y}_2 = 
\begin{bmatrix}
y^{(2)}_1 \\
\vdots \\
y^{(2)}_n
\end{bmatrix}, \quad
\boldsymbol{y} = 
\begin{bmatrix}
\boldsymbol{y}_1 \\
\boldsymbol{y}_2
\end{bmatrix}, \quad
\boldsymbol{X} = 
\begin{bmatrix}
\boldsymbol{0}_{n} & I_n \\
\boldsymbol{1}_{n} & I_n
\end{bmatrix}, \quad
\boldsymbol{\beta} = 
\begin{bmatrix}
\delta \\
\beta_1 \\
\vdots \\
\beta_n
\end{bmatrix}
```
と線形制約
```math
\boldsymbol{H}^T \boldsymbol{\beta} = d, \quad \boldsymbol{H} = \begin{bmatrix}
1 \\ \boldsymbol{0}_n
\end{bmatrix}, \quad d=0
```
に対して具体的に計算する。まず、$\boldsymbol{1}_n^T \boldsymbol{1}_n = n$なので、
```math
\begin{align}
(\boldsymbol{X}^T \boldsymbol{X})^{-1} = \begin{bmatrix}
n & \boldsymbol{1}_{n}^T \\
\boldsymbol{1}_{n} & 2 I_n
\end{bmatrix}^{-1} = \begin{bmatrix}
\frac{2}{n} & -\frac{1}{n} \boldsymbol{1}_{n}^T \\
-\frac{1}{n} \boldsymbol{1}_{n} & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T
\end{bmatrix}
\end{align}
```
である。
また、
```math
\boldsymbol{1}_n^T \boldsymbol{y}_1 = \sum_{i=1}^{n_1} y_{j}^{(1)} = n_1 \bar{y}_1
```
であるから、
```math
\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{y} \\
&=
\begin{bmatrix}
\frac{2}{n} & -\frac{1}{n} \boldsymbol{1}_{n}^T \\
-\frac{1}{n} \boldsymbol{1}_{n} & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{0}_{n}^T & \boldsymbol{1}_{n}^T \\
I_n & I_n
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{y}_1 \\
\boldsymbol{y}_2
\end{bmatrix}\\
&= 
\begin{bmatrix}
\frac{2}{n} & -\frac{1}{n} \boldsymbol{1}_{n}^T \\
-\frac{1}{n} \boldsymbol{1}_{n} & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T
\end{bmatrix}
\begin{bmatrix}
n \bar{y}_2 \\
\boldsymbol{y}_1 + \boldsymbol{y}_2
\end{bmatrix}\\
&=
\begin{bmatrix}
\bar{y}_2 - \bar{y}_1 \\
\frac{1}{2} \left( \boldsymbol{y}_1 + \boldsymbol{y}_2  - (\bar{y}_2 - \bar{y}_1)\boldsymbol{1}_{n} \right)
\end{bmatrix}
\end{align}
```
となる。より、$\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - d = \bar{y}_2 - \bar{y}_1$である。
```math
\begin{align}
(\boldsymbol{H}^T (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{H})^{-1} &= \left( \begin{bmatrix}
1 & \boldsymbol{0}_n^T
\end{bmatrix} 
\begin{bmatrix}
\frac{2}{n} & -\frac{1}{n} \boldsymbol{1}_{n}^T \\
-\frac{1}{n} \boldsymbol{1}_{n} & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T
\end{bmatrix}
\begin{bmatrix}
1 \\
\boldsymbol{0}_n
\end{bmatrix} \right)^{-1} = \frac{n}{2}
\end{align}
```
なので、$F$検定統計量の分子は、
```math
(\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d})^T (\boldsymbol{H}^T (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{H})^{-1} (\boldsymbol{H}^T \hat{\boldsymbol{\beta}} - \boldsymbol{d}) = \frac{n}{2} (\bar{y}_2 - \bar{y}_1)^2
```
となる。
次に、分母について、
```math
\begin{align}
I_{2n} - \boldsymbol{X} (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T 
&= I_{2n} - \begin{bmatrix}
\boldsymbol{0}_{n} & I_n \\
\boldsymbol{1}_{n} & I_n
\end{bmatrix}
\begin{bmatrix}
\frac{2}{n} & -\frac{1}{n} \boldsymbol{1}_{n}^T \\
-\frac{1}{n} \boldsymbol{1}_{n} & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{0}_{n}^T & \boldsymbol{1}_{n}^T \\
I_n & I_n
\end{bmatrix} \\
&=
I_{2n} - 
\begin{bmatrix}
\boldsymbol{0}_{n} & I_n \\
\boldsymbol{1}_{n} & I_n
\end{bmatrix}
\begin{bmatrix}
-\frac{1}{n_1} \boldsymbol{1}_{n_1}^T & \frac{1}{n_1} \boldsymbol{1}_{n_1}^T \\
\frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & \frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
\end{bmatrix} \\
&=
I_{2n} - 
\begin{bmatrix}
\frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & \frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
\frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
\end{bmatrix} \\
&=
\begin{bmatrix}
\frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & - \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
- \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & \frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
\end{bmatrix} 
\end{align}
```
であるから、
```math
\begin{align}
&\boldsymbol{y}^T (I_{2n} - \boldsymbol{X} (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T) \boldsymbol{y} \\
&= 
\begin{bmatrix}
\boldsymbol{y}_1^T & \boldsymbol{y}_2^T
\end{bmatrix} 
\begin{bmatrix}
\frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & - \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
- \frac{1}{2} I_n + \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T & \frac{1}{2} I_n - \frac{1}{2n} \boldsymbol{1}_{n} \boldsymbol{1}_{n}^T \\
\end{bmatrix} 
\begin{bmatrix}
\boldsymbol{y}_1 \\ \boldsymbol{y}_2
\end{bmatrix} \\
&= \frac{1}{2} \boldsymbol{y}_1^T \boldsymbol{y}_1 - \frac{n}{2} \bar{y}_1^2 - \boldsymbol{y}_1^T \boldsymbol{y}_2 - n \bar{y}_1 \bar{y}_2 + \frac{1}{2} \boldsymbol{y}_2^T \boldsymbol{y}_2 - \frac{n}{2} \bar{y}_2^2 \\
&= \frac{1}{2} (\boldsymbol{y}_1 - \boldsymbol{y}_2)^T (\boldsymbol{y}_1 - \boldsymbol{y}_2) - \frac{n}{2} (\bar{y}_1 - \bar{y}_2)^2
\end{align}
```
となる。$\boldsymbol{y}_1 - \boldsymbol{y}_2$は対応する観測に関して1つ目の観測と2つ目の観測の差のベクトルである。ここで、対応する観測の1つ目の観測と2つ目の観測の差を
```math
\Delta y_i = y^{(2)}_{i} - y^{(1)}_{i}, \quad i=1,...,n
```
と表す。$\Delta y_i, i=1,...,n$の標本平均$\bar{\Delta y}$は$\bar{\Delta y} = \bar{y}_2 - \bar{y}_1$であり、
```math
\begin{align}
\boldsymbol{y}^T (I_{2n} - \boldsymbol{X} (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T) \boldsymbol{y} &= \frac{1}{2} \sum_{i=1}^{n} \Delta y_i^2 - \frac{n}{2}  \bar{\Delta y}^2 \\
&= \frac{n-1}{2} S_{\Delta y}
\end{align}
```
となる。ただし、$S_{\Delta y}$は$\Delta y_1,...,\Delta y_n$の不偏分散。
$F$検定統計量の$n,p,r$は$n=2n, p=n+1, r=1$に対応するため、$F$検定統計量は、
```math
F = \frac{\frac{n}{2} (\bar{y}_2 - \bar{y}_1)^2}{\frac{n-1}{2} S_{\Delta y} / (2n - (n+1))} = \frac{\bar{\Delta y}^2}{ \frac{1}{n} S_{\Delta y}}
```
になる。帰無仮説$\delta=0$が$F$検定で棄却される条件は、
```math
\frac{\bar{\Delta y}^2}{ \frac{1}{n} S_{\Delta y}} \ge F^{(1, n-1)}_{1-\alpha}
```
となる。自由度$n-1$のstudentの$t$分布に従う確率変数を2乗した確率変数は自由度$(1, n-1)$の$F$分布に従うため、
```math
\left( t^{n-1}_{1-\frac{\alpha}{2}} \right)^2 = F^{(1, n-1)}_{1-\alpha}
```
従って、帰無仮説$\beta_1 = 0$が棄却される$F$検定の条件は、
```math
\frac{|\bar{\Delta y} |}{ \sqrt{\frac{1}{n} S_{\Delta y}}} \ge t^{n-1}_{1-\frac{\alpha}{2}}
```
となる。これは、$n$個の対応する観測について1つ目と2つ目の差を観測値として、観測値の母平均が$0$であることを帰無仮説とする$t$検定の帰無仮説棄却の条件と等しい。