記事の目的

変量効果モデル(線形混合モデル)のパラメータの検定・信頼区間推定の具体的なアルゴリズムに関して、数理的な導出の解説を残す。
変量効果モデルのパラメータ推定方法は、最尤推定法、制限付き最尤推定法、ベイズ推定などいくつかある。この記事では、変量効果モデルの各パラメータについてそれ以外のパラメータに依存しない検定統計量を導き、単純帰無仮説の仮説検定を基にした信頼区間推定の方法を解説する。仮説検定を基にした信頼区間推定の方法の一般論については私の過去の記事（リンク）で解説している。

問題設定

変量効果モデルを以下のように定義する。

y = X \boldsymbol{\beta} + B \omega + \varepsilon, \quad \omega \sim N(\boldsymbol{0}_m, \rho \sigma^2 I_m),  \quad \varepsilon \sim N(\boldsymbol{0}_n, \sigma^2 I_n)

y \in \mathbb{R}^{n}, X \in \mathbb{R}^{n \times p}, \boldsymbol{\beta} \in y \in \mathbb{R}^{p},B \in \mathbb{R}^{n \times m}, \sigma^2 > 0, \rho \ge 0

ただし、$I_n$は$n$次元の単位行列、$\boldsymbol{0}_n$は$n$次元の零ベクトルを表す。$y$がサイズ$n$の観測、$X$は各観測の固定効果を表す既知の非確率的な行列で、$\boldsymbol{\beta}$は次元$p$の固定効果のパラメータである。
$B$は変量効果の構造を表す既知の非確率的な行列であり、$\rho$は測定誤差$\varepsilon$の分散と変量効果の確率変数$\omega$の分散の比を表す。
変量効果モデルでは、$\omega, \varepsilon$は独立と仮定する。観測された$y, X, B$から、$\boldsymbol{\beta}, \sigma, \rho$を推定することが、変量効果モデルの目的となる。
$\epsilon = B \omega + \varepsilon$とすれば、$\omega, \varepsilon$は独立な正規分布に従うので、

\epsilon \sim N(\boldsymbol{0}_n, \sigma^2 (I_n + \rho B B^T))

となる。$y = X \boldsymbol{\beta} + \epsilon$だから、

y \sim N(X \boldsymbol{\beta}, \sigma^2 (I_n + \rho B B^T))

である。
ここでは、簡単のため$X, B$に関してそれぞれ冗長性がない、つまり

\displaylines{
\mathrm{rank}(X) = p \\
\mathrm{rank}(B) = m \\
}

を仮定する。これは、$X$と$B$の列ベクトルが線型独立であることを意味する。これを仮定しないとパラメータが不定となって、有界な信頼区間を持たない場合がある。また、

\displaylines{
\mathrm{rank}([X \quad B]) = p + m - l \\
}

としておく。$0 \le l \le \min(p, m)$である。ここでは固定効果の$X$と変量効果の$B$が区別されるよう$0 \le l < \min(p, m)$を仮定する。

パラメータの検定統計量・信頼区間の導出

σ^2に関する検定統計量・信頼区間

$\beta, \rho$が未知の場合に、$\sigma^2$についての単純帰無仮説に関する検定統計量を導く。

\mathrm{rank}(
\begin{bmatrix} 
X & B
\end{bmatrix}
) = p + m - l

の仮定から、線型写像の次元定理より

\mathrm{dim} \left( \mathrm{Ker}\left(
\begin{bmatrix} 
X^T \\
B^T
\end{bmatrix}
\right) \right) = n - m - p + l

である。ただし、$\mathrm{Ker}$はベクトルへ左から積をとる線型写像として行列を捉えた場合の線型写像の零空間を表す。つまり、

\mathrm{Ker}\left(
\begin{bmatrix} 
X^T \\
B^T
\end{bmatrix}
\right) = \left\{v \in \mathbb{R}^{n} \middle| \begin{bmatrix} 
X^T \\
B^T
\end{bmatrix} v = \boldsymbol{0}_{m+p}\right\}

を意味する。この$n - m - p + l$次元の$\mathbb{R}^{n}$の線型部分空間の$n-m-p+l$本の正規直交基底を列ベクトルとして並べた行列を$V \in \mathbb{R}^{n \times (n-m-p+l)}$とする。$V$は$[X \quad B]$を特異値分解することなどで具体的に計算できる。特異値分解を使った零空間の正規直交基底の計算方法は私の過去の記事（リンク）に解説を載せている。
この$V$を得る方法は$\beta, \rho$に依存せず、$B, X$のみに依存する。この方法で得た$V$は、

\begin{align}
X^T V &= O_{p, n-m-p+l} \\
B^T V &= O_{m, n-m-p+l} \\
V^T V &= I_{n-m-p+l}
\end{align}

の性質を持つ。ただし、$O_{p, n-m-p}$は$p \times (n-m-p)$の零行列を表す。
これにより、

\begin{align}
V^T y = V^T X \boldsymbol{\beta} + V^T \epsilon = V^T\epsilon
\end{align}

であり、$\sigma^2V^T(I_n + \rho B B^T)V = \sigma^2 I_{n-m-p+l}$なので、

V^T y \sim N(\boldsymbol{0}_{n-m-p+l}, \sigma^2 I_{n-m-p+l})

となる。従って、$\sigma^2$が真の値であるという仮定の下で、$\beta, \rho$によらず。

\frac{1}{\sigma^2} y^T V V^T y \sim \chi^2 (n-m-p+l)

となる。
これを用いることで、任意の値$\hat{\sigma}^2$について単純帰無仮説$\sigma^2 = \hat{\sigma}^2$の有意水準$\alpha$の両側検定が観測$y, X, B$の下で棄却されない条件は、

{\chi^2}^{(n-m-p+l)}_{\frac{\alpha}{2}} < \frac{1}{\hat{\sigma}^2} y^T V V^T y < {\chi^2}^{(n-m-p+l)}_{1 - \frac{\alpha}{2}}

となる。従って、

\left( \frac{y^T V V^T y}{{\chi^2}^{(n-m-p+l)}_{1 - \frac{\alpha}{2}} }, \, \frac{y^T V V^T y}{{\chi^2}^{(n-m-p+l)}_{\frac{\alpha}{2}}} \right)

は$\sigma^2$の$(1-\alpha)100\%$信頼区間となる。
また、

\frac{1}{\sigma^2} E[y^T V V^T y] = n-m-p+l

であるため、$\frac{1}{n-m-p+l} y^T V V^T y$は$\sigma^2$の不偏推定量になる。

ρに関する検定統計量・信頼区間

$\mathrm{rank} (X) = p$の仮定と線型写像の次元定理から、$\mathrm{dim}(\mathrm{Ker}(X^T)) = n - p$である。$V$の列ベクトルは$\mathrm{Ker}(X^T)$に含まれるので、$V$の$n-m-p+l$本の列ベクトルとそれらに直交する$m-l$本の互いに直交する大きさ$1$のベクトルで$\mathrm{Ker}(X^T)$の正規直交基底を構成できる。その$m-l$本のベクトルを列ベクトルとして並べた行列を$U \in \mathbb{R}^{n \times (m-l)}$とする。$U$は、

\begin{align}
X^T U &= O_{p, m-l} \\
U^T V &= O_{m-l, n-m-p+l} \\
U^T U &= I_{m-l}
\end{align}

の性質を持つ。

$U$は以下の方法で$B, X$から具体的に計算できる。まず、$\mathrm{rank} ([X \quad V]) = n-m+l$である。なぜなら、$\mathrm{rank} ([X \quad V]) < n-m+l$の場合、$[X \quad V]$の列ベクトルが線型独立でないので、

Xw_X + V w_V = 0_n

となる$w_X \in \mathbb{R}^{p}$, $w_X \in \mathbb{R}^{n-m-p+l}$, $[w_X, w_V]^T \neq 0_{n-m+l}$が存在する。これに前から$V^T$をかけることで

V^T (Xw_X + V w_V) = w_V = 0_{n-m-p+l}

となる。よって、$w_X \neq 0_p$が存在して$Xw_X = 0_n$となるが、これは$\mathrm{rank}(X) = p$に矛盾する。従って、$\mathrm{rank} ([X \quad V]) = n-m+l$となる。
線型写像の次元定理から

\mathrm{dim} \left( \mathrm{Ker} \left(\begin{bmatrix} X^T \\ V^T\end{bmatrix}\right) \right) = m-l

となるので、この零空間の$m-l$本の正規直交基底を列ベクトルとして並べた行列を$U$とすれば、$U$の性質を全て満たす。

ここで、

\begin{align}
U^T y = U^T X \boldsymbol{\beta} + U^T \epsilon = U^T \epsilon
\end{align}

であり、$U^T (I_n + \rho B B^T) U = I_m + \rho U^T B B^T U$なので、

U^T y \sim N(\boldsymbol{0}_{m-l}, \sigma^2 (I_{m-l} + \rho U^T B B^T U^T))

となる。後述する理由により$(I_{m-l} + \rho U^T B B^T U)$は正定値行列であり正則なので、

\frac{1}{\sigma^2} y^T U (I_{m-l} + \rho U^T B B^T U)^{-1} U^T y \sim \chi^2 (m-l)

となる。
$U^T V = O_{m-l, n-m-p}$なので、$U^T y$と$V^T y$は独立である。従って、真の$\beta, \sigma^2$に依存せず、$\rho$が真の値であれば、

\frac{y^T U (I_{m-l} + \rho U^T B B^T U)^{-1} U^T y / (m-l)}{y^T V V^T y / (n-m-p+l)} \sim F(m-l, n-m-p+l)

となる。これを用いることで、任意の値$\hat{\rho}$について単純帰無仮説$\rho = \hat{\rho}$の有意水準$\alpha$の両側検定が観測$y, X, B$の下で棄却されない条件は、

F^{(m-l, n-m-p)}_{\frac{\alpha}{2}} < \frac{y^T U (I_{m-l} + \hat{\rho} U^T B B^T U)^{-1} U^T y / (m-l)}{y^T V V^T y / (n-m-p)} < F^{(m-l, n-m-p)}_{1 - \frac{\alpha}{2}}

となる。これを満たす$\hat{\rho}$の範囲が観測$y, X, B$の下での$\rho$の$(1-\alpha)100\%$信頼区間である。陽には解けないが、下記の理由で$F$検定統計量が$\rho$に対して単調減少なので、不等式の両端についてNewton法などで解けば、有界な信頼区間を得られる。

I_m + ρ U^T B B^T Uが正定値行列である理由

任意の$v \in \mathbb{R}^{m-l}$について、

v^T (U^T B B^T U) v = (B^T U v)^T (B^T U v) = \| B^T U v \|^2 \ge 0

となるので、$U^T B B^T U \in \mathbb{R}^{(m-l) \times (m-l)}$は半正定値行列である。さらに、$\mathrm{rank}(B^T U)=m-l$であれば、$U^T B B^T U$は正定値行列である。
$\mathrm{rank}(B^T U) < m-l$を仮定すると、$v \neq \boldsymbol{0}_{m-l}$が存在して$B^T U v = \boldsymbol{0}_m$となる。$U$の性質として、

X^T U v = O_{p, m-l} v = \boldsymbol{0}_p

が成り立つので、

U v \in \mathrm{Ker}\left(
\begin{bmatrix} 
X^T \\
B^T
\end{bmatrix}
\right)

となる。よって、$U v = V w$を満たす$w \in \mathbb{R}^{n-m-p+l}$, $w \neq \boldsymbol{0}_{n-m-p+l}$が存在する。これは$U$の列ベクトルと$V$の列ベクトルが線型独立であることに矛盾する。従って、$\mathrm{rank}(B^T U) = m-l$であり、$U^T B B^T U$は正定値行列。
$\tau$を$U^T B B^T U$の任意の固有値とし、その固有ベクトルを$v$とすると、正定値性から$\tau>0$であり、

(I_m + \rho U^T B B^T U)v = (1 + \rho \tau) v

となるから、$1 + \rho \tau$は$I_m + \rho U^T B B^T U$の固有値。$\rho, \tau > 0$なので、$I_m + \rho U^T B B^T U$は正定値行列である。

F検定統計量がρに対して単調減少の理由

$U^T B B^T U$が正定値行列であるので、$U^T B B^T U$の$m-l$個の固有値は全て正である。それを$\tau_1, ..., \tau_{m-l} > 0$として、固有ベクトルを列ベクトルとして並べた行列を$W \in \mathbb{R}^{(m-l) \times (m-l)}$とする。$U^T B B^T U$が対称行列なので、$W$は直交行列。$U^T B B^T U$の対角化は

U^T B B^T U = W \mathrm{diag}(\tau_1, ..., \tau_{m-l}) W^T

と表すことができ、$W$は直交行列なので、$I_m = W W^T$であるから、

I_{m-l} + \hat{\rho} U^T B B^T U = W \mathrm{diag}(1+\rho \tau_1, ..., 1+\rho \tau_{m-l}) W^T

となる。
これにより、

y^T U (I_{m} + \hat{\rho} U^T B B^T U)^{-1} U^T y= y^T U W \mathrm{diag}\left( \frac{1}{1+\rho \tau_1}, ..., \frac{1}{1+\rho \tau_{m-l}} \right) W^T U^T y

となる。$W^T U^T y = [y'_1, ...,y'_m]$とすれば、

\frac{y^T U (I_{m} + \rho U^T B B^T U)^{-1} U^T y / (m-l)}{y^T V V^T y / (n-m-p+l)} = \frac{n-m-p+l}{(m-l) y^T V V^T y} \sum_{j=1}^{m-l} \frac{y'^2_j}{1+\rho \tau_j}

となり、$\tau_1, ..., \tau_{m-l} > 0$なので、これは$\rho$に対して単調減少する。

βに関する検定統計量・信頼区間

$\mathrm{rank} (B) = m$を仮定しているので、線型写像の次元定理から、$\mathrm{dim}(\mathrm{Ker}(B^T)) = n - m$である。$V$の列ベクトルは$\mathrm{Ker}(B^T)$に含まれるので、$V$の$n-m-p+l$本の列ベクトルとそれらに直交する$p-l$本の互いに直交する大きさ$1$のベクトルで$\mathrm{Ker}(B^T)$の正規直交基底を構成できる。その$p-l$本のベクトルを列ベクトルとして並べた行列を$S \in \mathbb{R}^{n \times (p-l)}$とする。$S$は、

\begin{align}
B^T S &= O_{m, p-l} \\
S^T V &= O_{p-l, n-m-p+l} \\
S^T S &= I_{p-l}
\end{align}

の性質を持つ。ちょうど$\rho$についての検定統計量を導出した際の$U$について、$X$と$B$、$p$と$m$を入れ替えて導いているのが$S$である。

$S$を具体的に計算する方法も、$U$を具体的に計算する方法について$X$と$B$、$p$と$m$を入れ替えることで得られる。

$y = X \boldsymbol{\beta} + B \omega + \varepsilon$に対して、右から$S^T$をかける。$S^T B = O_{p-l, m}$なので、

S^T y = S^T X \boldsymbol{\beta} + S^T \varepsilon

となり$\rho$を消去できる。$S^T S = I_{p-l}$なので、$S^T y \sim N(S^T X \boldsymbol{\beta}, \sigma^2 I_{p-l})$となる。従って、$\beta$の推定は$S^T X$を独立変数、$S^T y$を観測とする線形モデルに帰着する。
ただし、$\mathrm{rank}(B^T U) = m-l$と同様の議論により、$\mathrm{rank}(S^T X) = p-l$であるから、$\mathrm{rank}(X^T S S^T X) = p - l$である。$l = 0$であれば$X^T S S^T X$は逆行列を持つが、$l \ge 0$ならば逆行列を持たないので擬似逆行列で議論を進めることとなる。$X^T S S^T X$の擬似逆行列を$(X^T S S^T X)^{-}$と表記する。
ある$H \in \mathbb{R}^{p \times r}$, $d \in \mathbb{R}^r$に対して、$\boldsymbol{\beta}$の真の値について線形制約$H^T \boldsymbol{\beta} = d$が成立しているかを検定するための検定統計量を導出する。この検定が可能な条件は$\mathrm{rank}(H)=r$と$H = (X^T S S^T X) A$となる$A \in \mathbb{R}^{n \times r}$が存在することである。これが成り立つ時、

\begin{align}
(X^T S S^T X) (X^T S S^T X)^{-} H &= (X^T S S^T X) (X^T S S^T X)^{-} (X^T S S^T X) A \\
&= (X^T S S^T X) A \\
&= H
\end{align}

となる。
$S^T y = S^T X \boldsymbol{\beta} + S^T \varepsilon$についての正規方程式$(X^T S S^T X) \boldsymbol{\beta} = X^T S S^T y$を満たす$\boldsymbol{\beta}$を$\tilde{\boldsymbol{\beta}}$とする。$l = 0$であれば$\tilde{\boldsymbol{\beta}}$は唯一つに定まるが、$l \ge 0$では一つに定まらない。
ここで、$\boldsymbol{\beta}' \in \mathbb{R}^p$として$H^T \boldsymbol{\beta}' = d$を見たすベクトルを任意に1つ取ると、

\begin{align}
H^T \tilde{\boldsymbol{\beta}} &= H^T (X^T S S^T X)^{-} (X^T S S^T X) \tilde{\boldsymbol{\beta}} \\
&= H^T (X^T S S^T X)^{-} X^T S S^T y\
\end{align}

であり、$S^T y \sim N(S^T X \boldsymbol{\beta}, \sigma^2 I_{p-l})$と、

\displaylines{
H^T (X^T S S^T X)^{-} X^T S S^T X \boldsymbol{\beta} = H^T \boldsymbol{\beta} \\
H^T (X^T S S^T X)^{-} X^T S (H^T (X^T S S^T X)^{-} X^T S)^T = H^T (X^T S S^T X)^{-} H
}

から、

H^T \tilde{\boldsymbol{\beta}} - d \sim N(H^T \boldsymbol{\beta} - d, H^T (X^T S S^T X)^{-} H)

になる。$H^T (X^T S S^T X)^{-} H$は、

\begin{align}
\mathrm{rank}(H^T (X^T S S^T X)^{-} H) &= \mathrm{rank}(H^T (X^T S S^T X)^{-} (X^T S S^T X) (X^T S S^T X)^{-}H) \\
&= \mathrm{rank}(S^T X (X^T S S^T X)^{-}H) \\
&= \mathrm{rank}(X^T SS^T X (X^T S S^T X)^{-}H) \\
&= \mathrm{rank}(H) = r
\end{align}

となり正則である。従って、帰無仮説$H^T \boldsymbol{\beta} = d$が真ならば、

\frac{1}{\sigma^2} (H^T \tilde{\boldsymbol{\beta}} - d)^T (H^T (X^T S S^T X)^{-} H) ^{-1} (H^T \tilde{\boldsymbol{\beta}} - d)^T \sim \chi^2(r)

となる。また、$S^TB=O_{p-l,m}$, $S^T V=O_{p-l,n-m-p+l}$, $X^T V=O_{p,n-m-p+l}$であり、$y \sim N(X \boldsymbol{\beta}, \sigma^2 (I_n + \rho B B^T))$なので、

\begin{align}
E\left[(H^T \tilde{\boldsymbol{\beta}} - d) (V^T y)^T\right] &= E \left[ H^T (X^T S S^T X)^{-} X^T S (S^T y - S^T X \boldsymbol{\beta}') y^T V \right] \\
&=H^T (X^T S S^T X)^{-} X^T S S^T E[y y^T] V \\
&\quad - H^T (X^T S S^T X)^{-} X^T S S^T X \boldsymbol{\beta}' E[y]^TV\\
&=H^T (X^T S S^T X)^{-} X^T S S^T (\sigma^2(I_n + \rho B B^T) + X \boldsymbol{\beta} \boldsymbol{\beta}^T X^T) V \\
&\quad - H^T (X^T S S^T X)^{-} X^T S S^T X \boldsymbol{\beta}' \boldsymbol{\beta}^T X^T V\\
&= O_{r, n-m-p+l}
\end{align}

となる。つまり、$H^T \tilde{\boldsymbol{\beta}} - d$と$V^T y$は独立である。
従って、真の$\boldsymbol{\beta}$に対して、線形制約$H^T \boldsymbol{\beta} = d$が真の場合、

\frac{(H^T \tilde{\boldsymbol{\beta}} - d)^T (H^T (X^T S S^T X)^{-} H) ^{-1} (H^T \tilde{\boldsymbol{\beta}} - d)^T / r}{y^T V V^T y / (n-m-p+l)} \sim F(r, n-m-p+l)

となる。

この$F$分布に従う検定統計量を用いることで、$\boldsymbol{\beta}$の個々の要素の信頼区間を求めることができる。
例えば、$\boldsymbol{\beta} = [\beta_1,...,\beta_p]^T$とした場合の$\beta_1$の$(1-\alpha)100\%$信頼区間を、$\beta_2,...,\beta_p, \sigma^2, \rho$に依存することなく求めたい場合、

r=1, \quad H = \begin{bmatrix}
1 \\ 0 \\ \vdots \\0
\end{bmatrix} \in \mathbb{R}^{p}, \quad d=b_1

として、$H^T \boldsymbol{\beta} = d \Leftrightarrow \beta_1 = b_1$についての有意水準$\alpha$の検定が棄却されない$b_1$の範囲を求めればよい。この範囲は、

\frac{(H^T (X^T S S^T X)^{-} H) ^{-1} (H \tilde{\boldsymbol{\beta}} - b_1)^2}{y^T V V^T y / (n-m-p+l)} < F_{1 - \alpha} (1, n-m-p+l)

を満たす$b_1$の範囲となる。$H \tilde{\boldsymbol{\beta}}$は$\tilde{\boldsymbol{\beta}}$の第1成分である。これを$\tilde{\beta}_1$としておく。

\tilde{\beta}_1 \sim N(\beta_1, H^T (X^T S S^T X)^{-} H)

であり、$\tilde{\beta}_1$は$\beta_1$の不偏推定量。

分母の$\frac{1}{n-m-p+l} y^T V V^T y$は$\sigma^2$の不偏推定量である。これを$\hat{\sigma^2}$とおく。

$H^T (X^T S S^T X)^{-1} H$は$p \times p$行列$(X^T S S^T X)^{-}$の第$(1, 1)$成分である。これを$[(X^T S S^T X)^{-}]_{(1,1)}$とする。

$X^T S S^T X$が半正定値行列なので、$(X^T S S^T X)^{-}$も半定値行列である。そのため、その対角成分である$[(X^T S S^T X)^{-}]_{(1,1)}$は非負である。

$H$の条件から$H^T (X^T S S^T X)^{-1} H$が正則なので、$[(X^T S S^T X)^{-}]_{(1,1)}$は正である。$\beta_1 = b_1$についての有意水準$\alpha$の検定が棄却されない$\hat{\beta}_1$の範囲は、

F^{(1, n-m-p+l)}_{\frac{\alpha}{2}} < \frac{(\tilde{\beta}_1 - b_1)^2}{[(X^T S S^T X)^{-}]_{(1,1)}\hat{\sigma^2}} < F^{(1, n-m-p+l)}_{1 - \frac{\alpha}{2}}

となる。

自由度$n-m-p$の$t$分布に従う確率変数を二乗した確率変数は自由度$(1, n-m-p+l)$の$F$分布に従うので、

\left(t^{(n-m-p+l)}_{1-\frac{\alpha}{2}}\right)^2 = F^{(1, n-m-p+l)}_{1 - \alpha}

であるから、$\beta_1$の$(1-\alpha)100\%$信頼区間は、

\left( \tilde{\beta}_1 - t^{(n-m-p+l)}_{1-\frac{\alpha}{2}} \sqrt{[(X^T S S^T X)^{-1}]_{(1,1)}\hat{\sigma^2}}, \, \tilde{\beta}_1 + t^{(n-m-p+l)}_{1-\frac{\alpha}{2}} \sqrt{[(X^T S S^T X)^{-1}]_{(1,1)}\hat{\sigma^2}} \right)

となる。$\beta_2,...,\beta_p$についても同様に、$\tilde{\boldsymbol{\beta}}$の各成分と$(X^T S S^T X)^{-1}$の対角成分から$t$分布の分位点に基づいた信頼区間を計算できる。