はじめに

次のような混合モデルを考えます。

y=X\beta+Zu+e

u \sim MVN(0,K\sigma_{g}^2), e \sim MVN(0,I\sigma_{e}^2)

この時、このモデルのβとuのBLUE及びBLUPは次の方程式を解くことで求めらることが分かっています(混合モデル方程式)。

\begin{bmatrix} X^{T}X & X^{T}Z \\ Z^{T}X & Z^{T}Z+G^{-1} \\ \end{bmatrix}\begin{bmatrix} \beta \\ u \\ \end{bmatrix}=\begin{bmatrix} X^{T}y \\ Z^{T}y \\ \end{bmatrix}

また、この時のパラメータの分散も次のように求まることが知られています。

\begin{bmatrix} Var(\beta) & Cov(\beta, u)\\ Cor(\beta, u) & Var(u) \\ \end{bmatrix}=\begin{bmatrix} X^{T}X & X^{T}Z \\ Z^{T}X & Z^{T}Z+G^{-1} \\ \end{bmatrix}^{-1}\sigma_{e}^2

非常に覚えやすい形をしていて個人的に気に入っているのですが、今回の私はこれに関する話題です。

一般に線形回帰 $ y_{i}=\beta_{0}+\beta_{1}x_{i}+e_{i}$ の解は最小二乗法を解くことにより求まります。

そしてこの最小二乗解は

\beta=(\beta_{0}, \beta_{1})^T , X= \begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}

としたとき、
$ X^{T}X \beta = X^{T}y$ を解くことにより求まるわけですがこれは混合モデル方程式のuが無いときに対応しています。

この解であるβの分散も知ることができれば傾きや切片自体に対して統計的検定(例えば傾きが有意に0と異なるか、など)を行えるのですが、この分散はやや分かりにくい形をしています。

Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}

Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}

さて、最小二乗解が混合モデル方程式のuが無いバージョンに対応しているように、この分散も混合モデルの分散を求める式のuが無いところ$Var(\beta)=(X^{T}X)^{-1}\sigma_{2}^2$に対応していると考えられます。
そこで今回は、この線形回帰の分散を混合モデルの分散を求める方程式から求めてみようと思います。

問題設定

$ Var(\beta)=(X^{T}X)^{-1}\sigma_{2}^2$　から　$ Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}, Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}$　を求める。

計算

まず、

X= \begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}

であることを利用して $ X^{T} X$を計算します

X^{T}X=\begin{bmatrix} 1 & 1 & \dots & 1\\  x_{1} & x_{2} & \dots & x_{n}\end{bmatrix}\begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}=\begin{bmatrix} n & \Sigma x_{i} \\ \Sigma x_{i} & \Sigma x_{i}^2 \end{bmatrix}

2×2の行列の逆行列はすでに公式が知られているのでこれに照らし合わせると

(X^{T}X)^{-1}=\begin{bmatrix} n & \Sigma x_{i} \\ \Sigma x_{i} & \Sigma x_{i}^2 \end{bmatrix}^{-1}=\frac{1}{n\Sigma x_{i}^2-(\Sigma x_{i})^2}\begin{bmatrix}  \Sigma x_{i}^2 & -\Sigma x_{i} \\  -\Sigma x_{i} & n  \end{bmatrix}

ここでよく知られた分散の公式から $ n\Sigma x_{i}^2-(\Sigma x_{i})^2 = n^2(\Sigma x_{i}^2/n-(\Sigma x_{i}/n)^2=n\Sigma(x_{i}-\bar{x})$ であるので、

(X^{T}X)^{-1}=\frac{1}{n\Sigma(x_{i}-\bar{x})}\begin{bmatrix}  \Sigma x_{i}^2 & -\Sigma x_{i} \\  -\Sigma x_{i} & n  \end{bmatrix}

よって、

Var(\beta)=\frac{\sigma_{e}^2}{n\Sigma(x_{i}-\bar{x})}\begin{bmatrix}  \Sigma x_{i}^2 & -\Sigma x_{i} \\  -\Sigma x_{i} & n  \end{bmatrix}

ここで以下が成り立つ

Var(\beta)=\begin{bmatrix} Var(\beta_{0}) & Cor(\beta_{0}, \beta_{1}) \\ Cor(\beta_{0}, \beta_{1}) & Var(\beta_{1}) \end{bmatrix}

よって、

 Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}

Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}

まとめ

今回、マトリックスを分解することで、混合モデルの分散の式から線形回帰のパラメータの分散を求めることできるとわかりました。
混合モデルの分散の式を覚えるのはそんな難しくはないので、パラメータの分散の式を忘れたときは、今後はここに立ち戻って計算しようと思います