霞と側杖を食らう

ほしいものです。なにかいただけるとしあわせです。[https://www.amazon.jp/hz/wishlist/ls/2EIEFTV4IKSIJ?ref_=wl_share]

制限付き最尤法(REML)の学習記録

学習記録

【学習動機】

  統計数理研究所の公開講座で、経時的反復測定データの取り扱い方で、MIxed models for repeated measures(MMRM法)について
習った。そのとき、パラメータ推定の仕方として、REML法が出てきたが、いまいちよく分かってなかったため。
制限付き最尤法(REML: REstricted Maximum Likelihood)（REMLは残差最尤法(REsidual Maximum Likelihood)ともいうみたい。罰則付き最尤法とは別物。）を、土居正明先生のpdf資料を用いて、学んでいく。
直交行列、Helmert変換、多変量正規分布、正準形あたりが前提知識となっており、それらも同様にして、準備していく。
その後、まだ理解が不足しているので、追加的に随時、補足する。

【学習記録】

土居正明『直交行列と Helmert 変換』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/orthogonal_mtx.pdf

直交行列の定義。 P’P=I_n
命題。n×n行列Pが直交行列であることはPが \mathbb{R}^nのある正規直交基底であることは同値。
性質。 O(n)={P: n×n行列|P’P|=I_n}として、
P \in O(n)に対して、群の性質 (a)-(c)に加えて、(d)と(e)がある。

(a) I_n \in O(n)
(b) P,Q \in O(n) \Rightarrow PQ \in O(n)
(c) P^{-1} \in O(n)
(d) P’ = P^{-1}
(e) |P| = \pm1

土居正明『一次独立と正規直交基底』講義資料
http://www012.upp.so-net.ne.jp/doi/math/anova/linear_indep.pdf
Parsevalの等式は射影と相性が良い。正準形の理論でも役立つ。
\mathbb{R}^n において正規直交基底\{ e_1, ... , e_n \} をとると、]
\mathbb{R}^n中の任意のベクトルaに対して、以下の(i)(ii)が成り立つ。
(i) a = (a \cdot e_2)e_1 + (a \cdot e_2)e_2 + ... + (a \cdot e_n)e_n
(ii) || a ||^2 = {(a \cdot e_1)}^2 + {(a \cdot e_2)}^2 + ... + {(a \cdot e_n)}^2

Parsevalの等式は射影と相性が良い。正準形の理論でも役立つ。

特殊な直交行列として、Helmert行列の作り方。
(i) ある手順で \mathbb{R}^nの直交基底を作る。
(ii) (i)で作った直交基底を正規直交基底に正規化する。
(iii) (ii)で作った正規直交基底を横に並べて直交行列にする。
(iv) (iii)で作った行列を転置する。

こうして作られたHelmert行列 H_nを用いた変換をHelmert変換という。
(エルミート行列(Hermitian matrix)のことか？と思ったけど、違うね。)
https://en.wikipedia.org/wiki/Helmert_transformation

これは、 1) \bar{X} と {s}^2 は互いに独立である。 2) n{s}^2/{\sigma}^2 ~{\chi}^2(n-1) を同時に示すのに使える。

c.f. 竹村数理統計p67,68

土居正明『多変量正規分布』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/m_normal.pdf

多変量正規分布は多変量標準正規分布を一次変換して、定数ベクトルを加えたもの。
「定理４：標準正規分布の標本平均と標本分散の独立性」
「系２：一般の正規分布の標本平均と標本分散の独立性」
上の1)と2)のこと。

土居正明『正準形１』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/canonical_form1.pdf

群の数を Iとし、 n_i (i=1,2,...,I)とおき、 \Sigma{n_i} = Nとする。
y= X \beta + \epsilon, \epsilon ~ N(0, \sigma ^2 I_N)
ただし、 \beta = (\mu_1, \mu_2, ..., \mu_I )

このとき、適当な直交行列G(Helmert行列で考えたのと同様)でyをzに変換すると
z_1 ~ N(\sqrt{n_1}\mu_1, \sigma ^2), z_I ~ N(\sqrt{n_I}\mu_I, \sigma ^2)
z_{I+1}, ..., z_N ~ N(0, \sigma ^2)
で、すべて独立で、
\hat{\mu_1} = z_1 / \sqrt{n_1},..., \hat{\mu_I} = z_I / \sqrt{n_I},
\hat{\sigma ^2} = \frac{1}{N-I} \Sigma _ {i=I+1} ^{N} z_i ^2
となり、 \hat{\mu_1}, \hat{\mu_2}, ..., \hat{\mu_I}, \hat{\sigma ^2}はすべて独立。
\hat{\mu_i} ~ N( \mu_i, \frac{1}{n_i} \sigma ^2)
(N-I) \frac{\hat{\sigma ^2}} {\sigma ^2} = \Sigma _ {i=I+1} ^{N} \frac{z_i ^2}{\sigma ^2} ~ \chi ^2 (N-I)
となる。

\beta と \sigma ^2 のUMVUである\hat{\beta}と \hat{\sigma ^2}が互いに独立であり、それぞれ従う分布が確認できた。

(UMVUとは)
\theta の不偏推定量 \hat{\theta} が一様最小分散不偏推定量(UMVU : Uniformly Minimum Variance Unbiased Estimator)であるとは \theta の任意の不偏推定量 \hat{\theta}に対して、
V_{\theta} (\hat{\theta}) \leq V_{\theta} (\hat{\theta}), \forall \theta
が成り立つときをいう。
(土居正明『不偏推定量・UMVU・大数の法則』資料２ http://www012.upp.so-net.ne.jp/doi/math/anova/UMVU.pdf より)

土居正明『正準形２』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/canonical_form2.pdf
分散分析のF検定の分子と分母の定数倍が帰無仮説のもとで独立なchi ^2分布に従うことを示して、
F統計量が帰無仮説のもとでF分布に従うことを示している。

土居正明『制限付き最尤法(REML)』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/REML.pdf REMLの問題意識は、正規分布から得られたデータからの、分散の最尤推定量は不偏推定量にならない。
正準形を用いてデータを２つに分け、平均と分散の別々の尤度関数を立てて推定する方法を制限付き最尤法という。
これは、xについての尤度関数を、正準形にしてzについての尤度関数の積に分解したものと考えられる。

長島健悟『線形推測論』
http://nshi.jp/contents/doc/glm/glm2016_10.pdf

ここまでの説明ではランダム効果の説明が入っていない。

長島先生がREMLの導出のメモをしているものを見つけました。
https://nshi.jp/contents/stat/reml/

REML以外にも良い資料がたくさんありました。
https://nshi.jp/contents/

南美穂子『制限付き最尤推定法(REML推定法)』 https://www.jstage.jst.go.jp/article/jappstat1971/25/2/25_2_73/_article/-char/ja
REMLに関する覚書き。この方は岩波DSで野球のバントの因果推論をしていたような気がします。

増田豊『線形混合モデルにおける分散成分の推定』分散成分の推定と歴史
http://www.obihiro.ac.jp/~suzukim/masuda2/stat/review_variance.html
帯広畜産大学のwebサイトのリニューアルで消えたらしいのですが、増田豊先生がミラーではてなブログに転載しているのを見つけました。

yutakamasuda.hatenablog.com

上の記事は2007年に書かれたもので、2019年の状況を踏まえた注釈の記事も書かれていました。
yutakamasuda.hatenablog.com

増田豊『制限付き最尤法』
https://slidesplayer.net/slide/11509221/
混合モデルのことを分散成分モデル(variance componets model)と呼んでいた。 REMLの推定のアルゴリズム
・導関数を使用しない方法（Derivative Free REML）
・1次導関数を使用する方法（EM REML）
・2次導関数を使用する方法（AI REML）
・AI REML以外にも，いくつかの準ニュートン法による分散成分の推定法が提案されている
・REML以外のアプローチとして、ベイズ。

黒木玄さんの不偏分散の分母がn-1の理由に関するマストドンのpdf https://genkuroki.github.io/documents/mathtodon/2017-05-16%20%E4%B8%8D%E5%81%8F%E5%88%86%E6%95%A3%E3%82%84%E4%B8%8D%E5%81%8F%E5%85%B1%E5%88%86%E6%95%A3%E3%81%AE%E5%AE%9A%E7%BE%A9%E3%81%AB%E3%81%8A%E3%81%84%E3%81%A6%E3%81%A9%E3%81%86%E3%81%97%E3%81%A6n%E3%81%A7%E5%89%B2%E3%82%89%E3%81%9A%E3%81%ABn%E2%88%921%E3%81%A7%E5%89%B2%E3%82%8B%E3%81%8B.pdf
n-1になるのは線形代数の直交射影で空間の次元が落ちるからという話。統計学の教科書はなかなか教えてくれない、自由度が何たるかを理解する助けになる。REMLの話も、上で見てきたように、同様の線形代数なので、この話もとても有用。

浅野・中村『計量経済学』
https://amzn.to/4dX2toj
の12章の12.2に「分散要素モデル」あり。あとで読む。

【学習予定】

REMLがまだよくわかっていない。ランダム効果モデル、混合モデルについて、もう少し詳しくなるべきで、
みどりぼん(https://amzn.to/4cZkMbL)やアヒル本(https://amzn.to/4dDaijm)、『一般化線形モデル入門』(https://amzn.to/4dD9lax)や『統計モデル入門』(https://amzn.to/4fVTpSx)を読もうと思う。