制限付き最尤法(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

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

土居正明『正準形1』分散分析資料
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・大数の法則』資料2 http://www012.upp.so-net.ne.jp/doi/math/anova/UMVU.pdf より)

土居正明『正準形2』分散分析資料
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の問題意識は、正規分布から得られたデータからの、分散の最尤推定量は不偏推定量にならない。
正準形を用いてデータを2つに分け、平均と分散の別々の尤度関数を立てて推定する方法を制限付き最尤法という。
これは、xについての尤度関数を、正準形にしてzについての尤度関数の積に分解したものと考えられる。

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

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

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

浅野・中村『計量経済学
Amazon CAPTCHA
の12章の12.2に「分散要素モデル」あり。あとで読む。

【学習予定】
REMLがまだよくわかっていない。ランダム効果モデル、混合モデルについて、もう少し詳しくなるべきで、
みどりぼんやアヒル本、『一般化線形モデル入門』や『統計モデル入門』を読もうと思う。

Rstudioとtidyverseの学習記録

【学習動機】

moratoriamuo.hatenablog.com

これのときに、気付いたのが、Hadley教徒になりたいけど、ぐぐって断片的に集めた知識じゃ、すぐ抜け落ちるってこと。というわけで、本で体系的に学ぼうということで、通称、宇宙本をやることにした。

 

松村 優哉『RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−』

http://amzn.asia/d/639ydJK

 

【学習記録】

<総括>

RStudioやtidyverseの入門として、とてもよかった。次に、各章の専門的に書かれた本に飛べばよいはず。

 

<忘れやすいショートカットめも>

Ctrl + Shift + M ;パイプ演算子%>%

Ctrl + Shift + C ;選択した行をコメント/非コメント

Ctrl + Alt + K :レンダリング実行(Knit)

Ctrl + Alt + I :Rチャンク挿入

 

c1:RStudioの基礎

・ファイルの読み込みは標準関数のread.csv()ではなく、readrのread_csv()を使うべし。

他に、data.tableのfread()も使えるが、ここでは扱わない。

標準関数のread.csv()は遅いし、読み込む型が不親切。

・読み込む列の型を指定したいなら、col_types引数を使う。

文字列[c]、整数[i]、実数[d]、TRUEFALSE[l]、日付[D]、時間[t]、日付時間[T]、数字以外の文字が含まれても数字として返す[n]、推測する[?]、列を読まない[_]

エンコーディングの指定。UTF-8で読み込まれるが、Shift-JIS(CP932)のファイルを読むと文字化けするので、local引数で指定。文字コードが分からないときは、guess_encoding()を使うといい。

・書き出しにはwrite_csv()を使う。標準関数のwrite.csv()より高速で、デフォルトでrow.names = FALSEなので、行名が入らなくてよい。

Excelファイルの読み込みにはreadxlのread_excel()を使えばよく、sheet引数でシートの指定ができる。

SASSPSS、STATAファイルの読み込みはhavenパッケージが使える。

 

c2:スクレイピングによるデータ収集

スクレイピングにはHTMLとCSSの知識必要。

・URLのHTML構造を抽出し、要素にアクセスし、その中身を取り出すか、属性にアクセスし、その値を取り出す。CSSセレクタXMLパスのような構文を利用したスクレイピングによってもドキュメントの要素や値を取得できる。

スクレイピングを実行するためにはrvestパッケージを使用。

まず、URLをread_htmlで読み込む。読み込むと、DOM(Document Object Model)という形で保存される。これはHTMLの要素を階層構造に変換したもの。HTMLの要素やクラスにあたる部分、DOMでノードという。DOMはブラウザでチェック可。

ページソースを読み解きやCSSセレクタXpathの文法暗記は困難なので、スクレイピングではデベロッパーツールというchromeの開発ツールを利用する。WindowsではFn+F12で使え、ノードが確認できる。ここから取得したいところにカーソルあてて、Copy XPathでコピー。クオーテーションマークはシングルとダブルを区別することに注意。

 

・ページを複数取りたい場合は、繰り返しを使う。URLが規則的な構造ならforループでなんとかなるが、そうでない場合はRSeleniumパッケージを使う。

(RSeleniumがCRANから消えている。2018/09/04現在。)

API(Application Programming Interface)は、つまり、ソフトウェアやサービスを外部から使用できるようにしたもの。無制限には使えず、単位時間での使用量決まっていること多い。

Twitterには、REST APIとSTREAMING APIがある。前者はタイムラインやDM、FF、リストなど取得できる。検索も可能。後者は流れるデータの1%を受信できる。

https://togetter.com/li/1250205

仕様変更していて、異なるところがあるみたい。

 

『Rによるスクレイピング入門』

 

で補強していくべき。

 

c3:dplyr/tidyrによるデータ前処理

 

Hadleyの提唱するtidy dataとは、

・1つの列が1つの変数を表す。

・1つの行が1つの観測を表す

・1つのテーブルが1つのデータセットだけを含む

 

(横長のワイド(wide)形式のデータを縦長のロング(long)形式のデータに変えるなど。パネルデータ分析の前処理で出てきたりする操作。)

ロング形式の方が、ベクトル演算しやすい。

データをtidy dataに変形するのに、tidyrのgather()とspread()が有用。

ワイドをロングに変形するのはgather()。

ロングをワイドに変形するのはspread()。

引数keyとvalueを指定、gather()なら並べる列名のkeyも指定、spread()はkeyが列名になるので指定不要。

 

列や行の絞り込みといったデータ操作はdplyrパッケージが有用、dplyrの主な関数は3つに大別でき、

・1つのデータフレームを操作する関数  行を絞るfilter()や列を絞るselect()など。

・2つのデータフレームを結合する関数  行同士で結合のjoin()系関数、集合演算intersect()やnion()など。

・ベクトルを操作する関数  ウィンドウ関数など。

 

・select()は引数に絞りたい列名を列挙。ヘルパ関数として、starts_with(), ends_with(), contains(), matches(), num_range(), one_of(),everything()が使える。

・filter()は引数に条件式を書いて絞り込み。論理式はいつも通り。

・mutate()は新しい列名を追加、transmute()は追加でなく別のtibble作成。

・arrange()が行を並べ替える。選択した列の昇順がデフォ。数値列なら-をつければ降順。そうでなければ、desc()でくくるといい。

・summarise()で集計計算。グループ化と組み合わせると強力。データを1行に要約する関数。

・group_by()を適用すると、データはグループ化されたデータフレーム(group_df)になり、mutate()やfilter()を使うと違う挙動を見せる。

 

セマンティクスに関しては、これ

https://speakerdeck.com/yutannihilation/dplyrfalseselecttomutatefalsesemanteikusufalsewei-i

 

・tidy dataに変換できないデータなら、複数の列を同時に処理するscoped functionが有用

・scoped functionは以上の関数の複数列バージョンで、サフィックスは_all, _if, _at。

これ、なかなか便利では!?

 

・2つのデータセットの結合と絞り込みで、_join()の関数を使う。

 

c4:ggplot2によるデータ可視化

カラーページで描かれていてよい。

パッケージggplot2の利点は、統一された文法、グラフの再現性、豊富な拡張パッケージ。

水準ごとの識別

・グループ化(mapping引数のaes()で指定。)

・ファセット(facet_wrap()とfacet_grid())

・統計的処理(stat())

水準ごとに自動で配色を割り当てる機能は可視化の効率化に大きく貢献する一方で、配布資料や印刷物でグラフを示す場合、モノクロやグレースケールが望ましいこともある。

背景真っ白なら、theme_classic()、目盛り線を残したいなら、theme_bw()、文字サイズやフォントの変更はtheme_()の中で変更。

 

 

c5:R Markdownによるレポート生成

ここのパートの著者が書いてるものが、

R Markdown入門

https://kazutan.github.io/kazutanR/Rmd_intro.html

 あと、これも役立つ。

www.amazon.co.jp

 

 

【学習予定】

 データハンドリングとスクレイピングをもう少し習得したい。

www.amazon.co.jp

 か

www.amazon.co.jp

 をやっていこうと思う。

dplyrの学習記録

【学習動機】

前の学習予定のデータ操作で便利らしいdplyrを学びたいというのが動機。

moratoriamuo.hatenablog.com

 

【学習記録】

dplyrはHadley WickhamというRの神様が作成したパッケージのひとつ。

dplyrは、でぃーぷらいあー、tidyrは、たいでぃあー、と読むらしい。

tidyverseをインストールすれば、dplyrの他、その他必要なパッケージもまとめてインストールすれば良い。

https://www.tidyverse.org

 

データの前処理に関して、下のスライドが読みやすく、よくまとまっていた。

www.slideshare.net

 

実際に使ってみながら

www.jaysong.net

 

この本が有用?これから読んでみる

https://www.amazon.co.jp/Rではじめるデータサイエンス-Hadley-Wickham/dp/487311814X

 

SQL経験のある人向けの記事

イメージ図付きで探している機能が見つけやすいかも。

qiita.com

qiita.com

 

qiita.com

 

それぞれの細かい説明(reshape2はdplyrなどの前身らしい。)

https://heavywatal.github.io/rstats/dplyr.html

https://heavywatal.github.io/rstats/tidyr.html

https://heavywatal.github.io/rstats/reshape2.html

 

チートシート(英語)

https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf

 

 

【学習予定】

 

松村 優哉 『RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界− 』

https://www.amazon.co.jp/R%E3%83%A6%E3%83%BC%E3%82%B6%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AERStudio-%E5%AE%9F%E8%B7%B5-%E5%85%A5%E9%96%80%E2%88%92tidyverse%E3%81%AB%E3%82%88%E3%82%8B%E3%83%A2%E3%83%80%E3%83%B3%E3%81%AA%E5%88%86%E6%9E%90%E3%83%95%E3%83%AD%E3%83%BC%E3%81%AE%E4%B8%96%E7%95%8C%E2%88%92-%E6%9D%BE%E6%9D%91-%E5%84%AA%E5%93%89/dp/4774198536

 

通称、宇宙本をやって体系的な理解を深めていきたい。

そんなあの子は透明少女

 『ペンギン・ハイウェイ』を見てきた。自分のなかにある記憶や物語やイメージがペンギンのように溢れだしてきて、たいへんな作品だった。溢れだしてきたものを思い出して、列挙してみる。順番は関係ない。『』すらつけない。輪るピングドラム、Highway chance、ポケモン、ハンマーセッション、セレビィ言の葉の庭サクラダリセット、コナン、君の名は、ぼくのなつやすみ宇宙よりも遠い場所、トトロ、GANTZ東京都美術館マリオサンシャインSHIROBAKO、雲のむこう約束の場所、千と千尋ハガレン夜は短し歩けよ乙女竹取物語ゼルダの伝説トワイライトプリンセスカリオストロの城化物語紅の豚エヴァあの日見た花の名前を僕達はまだ知らないトポロジー、ポニョ、25コ目の染色体思い出のマーニー時をかける少女、おーいでてこーい、いかにして問題をとくか。一通り、思い出せるだけでもこれくらいだろうか。

 少年の冒険はいつだって、心踊らされるテーマだ。ただし、この映画は、いつも出てくる無鉄砲で何も考えず突き進んでいく少年が主人公ではない。主人公は研究大好き少年なのだ。仮説を立て、調査し、実験し、観測し、仮説を検証する。複雑で難しいことは分割するし、理論として統合できることは統合する。問題点は冷静に反省し、次に活かしていく。そのような教育を父親から受けている。行き詰ったときは、一枚の紙に考えをまとめて書き、眺め、考え続ける、それでもダメなら、考えるのを一度やめてみるということすら教わっている。このような少年の冒険のストーリーは今まで一度も見たことがない。

 少年は大人ぶる。飲めないコーヒーを飲んで、背伸びをする。一方、自分は、あと半年もすれば、大人の仲間入りだ。モラトリアムの期間が終わろうとしている。最後の夏だ。平成最後の夏だと言われているが、人生最後の夏という気持ちで過ごしている。昨日より今日を、今日より明日を、より良く生きようという気持ちは少年と変わらない。

 以上で、触れなければならないはずの話に全く触れていない。そうだ、お姉さんだ。人類の永遠の謎だ。宇宙だ。母なる海だ。書きたいことを書くには理解と言葉と度胸が足りないので、ここでは触れたくないのだ。

 映画を見た後に購入したパンフレットにある森見登美彦へのインタビューによると、原作は少年時代の「原風景」を描いた作品だという。そこへの共感だろうか、強く引き込まれてしまって、心のなかに凝縮されていた何かが溢れだしたのだった。

Rでベイズ統計学の学習記録

【学習動機】

昔、MCMCやらなんやらは学んで、少し実装したことがあったが、

色々忘れていることも多いし、Rで手軽にできるようになっておきたいと思って、

J.アルバート『Rで学ぶベイズ統計学入門』

http://amzn.asia/d/hbKIgqd

をやっていくことにした。

 

【学習記録】

<総括>

LearnBayesパッケージにデータセットが豊富。

パッケージの重要な関数の作り方もあるのでわかりやすい。

練習問題。答えがあればなおよし。

密度や分布、周辺尤度あたりの言葉の使い方が気になった。和訳の問題かも。

やはり、BDA3を読むべき。

 

c1:R入門

LearnBayesパッケージ。

基本的な使い方の紹介。

 

c2:ベイズ的思考への誘い

・割合pのベイズ推定

事前分布の設定の仕方。

  • 離散事前分布の利用

LearnBayesのpdisc関数で計算。p discreteの意味か。

  • 連続事前分布(ベータ事前分布)の利用

共役の事前分布で、簡単に計算可能。

事後密度をカバーする区間でpの値のグリッドを選定する

このグリッド上で尤度と事前密度の積を計算する

それぞれの積を、それらの合計で割って正規化。これによりグリッド上の離散確率分布によって事後密度を近似している

Sample関数を使って、この離散分布からランダムな復元標本を抽出。

こうして得られたシミュレーション標本は事後分布からの標本を近似している。

 

予測

  • 離散事前分布

LearnBayesのpdiscp関数で予測密度計算。

  • ベータ事前分布

LearnBayesのpbetap関数で予測密度計算。

 

任意の事前分布から予測密度を計算できる便利な方法の一つがシミュレーションによるもの。

 

c3、c4 1パラメータモデル、複数パラメータモデル

 

c5:ベイズ計算入門

前2章では事後分布を要約する2タイプの方策。

ひとつは、サンプリング密度が指数分布族など周知の関数形ならば、パラメータを直接シミュレーションできる。もうひとつは周知の関数形でない場合、緻密な点のグリッド上で事後密度の値を計算して、これらの値が集中する離散事後密度によって、連続形の事後密度を近似できる(ブルートフォースと呼ばれる)。

本章では別のアプローチ。一般的なアプローチの一つが、モード周辺で事後分布の挙動に基づく方法で、多変量正規近似となり、最初の近似として役に立つ。続いて、事後分布の要約の計算にシミュレーションを利用する方法。事後分布が標準的な関数形ではない場合、適切に選ばれた提案密度の棄却サンプリング。重点サンプリングとSIR法のアルゴリズム積分計算し、一般の事後分布からシミュレーション。SIR法のアルゴリズムは、事前分布と尤度関数の変更に対して事後分布が敏感か調べられる。

 

事後分布の要約の多くは積分で表現できる。関数h(θ)の事後平均、h(θ)が集合Aに入る事後確率、関心のあるパラメータの周辺密度。これらは積分計算で求める場合、数値的に評価するが求積法は問題が限られる。

 

パラメータを再定式化して、近似しやすい形にする。

近似の一つは、事後モードにもとづく近似。多変量事後分布の要約を行う方法の一つが、そのモード周辺の密度の挙動を調べること

モード周りで多変量正規分布で近似するためにはモードを求める必要がある。

モードを求めるための汎用的な最適化アルゴリズムニュートン法やネルダー・ミード法。

 

棄却サンプリング、重点サンプリング、サンプリング重点リサンプリング。

小澄英男『ベイズ計算統計学

http://amzn.asia/d/6LIzKmB

C.P.ロバート、G.カセーラ『Rによるモンテカルロ法入門』

http://amzn.asia/d/e72cRHI

あたりで勉強しなおしたい。

c6:マルコフ連鎖モンテカルロ

棄却サンプリングは任意の事後分布からシミュレーションする一般的な方法であるが、適切な提案分布を構成する必要があるので、適用が困難。重点サンプリングとSIRもまた汎用的なアルゴリズムであるが、やはり提案分布を必要とし、特に高次元の問題の場合、難しい。そこで、MCMCによって事後分布を要約する方法。

 

c7:階層モデリング

 

心臓移植の死亡率データ。94の病院の心臓移植手術30日以内の死亡数データ。

死亡数y, 曝露数e。曝露1件あたりの死亡率λを推定したい。

すべての病院についての真の死亡率を同時推定したい。三つの方法。

一つ目の選択肢は、個々の死亡率で推定。 y_1/e_1,...,y_{94}/e_{94}

これは曝露数の少ない病院では不正確。ばらつきが大きい。

二つ目の選択肢は、真の死亡率はすべての病院で等しいという仮定で

\Sigma{y_j}/\Sigma{e_j}でプール推定。

三つ目の選択しは、一つ目と二つ目の複合。

(1-\lambda) y_i/e_i + \lambda \Sigma{y_j}/\Sigma{e_j}で、これは個別の推定値をプールされた推定値に縮める縮約推定値。

 

ベイズ分析では、モデルの仮定に変更を加えた場合の推論の感度を評価することが重要。

これにはサンプリング密度と事前密度についての仮定も含まれる。

5章のSIRは、ある事後分布からをシミュレーションした標本を、新しい分布に変換する便利な方法。

 

c8:モデル比較

あるパラメータについて二つの仮説を比較する方法として、ベイズファクター。

二つのベイズモデルの比較を行う。これは事前密度とサンプリング密度の選択に関わる。

この場合、ベイズファクターは二つのモデルの周辺密度の比で表される。

ベイズファクターは仮説の事前オッズに対する事後オッズの比で表される。

事前オッズ比は\pi_0 / \pi_1 = P(\theta_0 \in \Theta_0) / P(\theta_1 \in \Theta_1)

= \int_{\Theta_0}g(\theta)d\theta /  \int_{\Theta_1}g(\theta) d\theta

( g(\theta)は正則事前密度、g(\theta|y) \propto L(\theta)g(\theta)は事後密度 )

事後オッズは

p_0 / p_1 = P(\theta_0 \in \Theta_0 | y) / P(\theta_1 \in \Theta_1 | y)

= \int_{\Theta_0}g(\theta | y)d\theta / \int_{\Theta_1}g(\theta | y) d\theta

ベイズファクターは  BF = 事後オッズ/事前オッズ = (p_0/p_1) / (\pi_0/\pi_1)

統計量BFはデータによる証拠が仮説H0を支持する度合を示す基準。

仮説H0の事後確率はベイズファクターと仮説の事前確率の関数で表される。

 p_0 = \pi_0 BF / (\pi_0 BF + 1 -\pi_0)

 

仮説を比較するベイズの方法は、二つのモデル比較に一般化できる。

仮にyをデータベクトル、θをパラメータとすれば、ベイズモデルはサンプリング密度と事前密度を指定して構成される。

このモデルの周辺尤度は  m(y) = \int f(y|\theta)g(\theta)d\theta

M0とM1二つのベイズモデルを比較したいとする。

パラメータθの定義はモデル間で異なることもありうる。

このときモデルM0を支持するベイズファクターは、二つのモデルそれぞれに対するデータの周辺尤度の比になる。

 BF=m_0(y)/m_1(y)

とがモデルとそれぞれの事前確率を表すとすると、モデルの事前確率は以下で与えられる。

 P(M_0|y) = \pi_0 BF/(\pi_0 BF + \pi_1)

周辺尤度を近似する簡単な方法の一つがラプラス法。

\hat{\theta}が事後モードで、が対数事後密度のヘシアン(行列の2階微分)とする。このとき周辺尤度は以下で近似される。

ここでdはパラメータの数。

  

ベイズファクターの計算。

 

c9:回帰モデル

 

回帰モデル。

Zellnerのg事前分布によるモデル選択。

 

c10:ギブスサンプリング

これは別で学習するのでスルー

 

c11:RでWinBUGSを使うインターフェイス

BUGS言語、好きでないのでスルー

 

【学習予定】

『 StanとRでベイズ統計モデリング

www.amazon.co.jp

『実践ベイズモデリング

www.amazon.co.jp

あたりで、モデリングを習得していきたい。

 

 

 

 

 

 

 

一般化線形モデルの学習記録

【学習動機】

昔、『一般化線形モデル入門』

http://amzn.asia/fe0d2DC

を読んだけど、実際にコードは書いてないのと理解が浅いのがあって、

復習を含めて、Rでコードを書きながら学びなおそうと

『一般化線形モデル (Rで学ぶデータサイエンス 10)』

http://amzn.asia/7Q3Z8I9

をやっていくことにした。

 

【学習記録】

一般化線形モデルのセットアップ。

c1

・構成要素が線形予測子、リンク関数、誤差構造。

c2

・glm関数の基本

・最尤法

・Wald検定、スコア検定、尤度比検定

・尤度比検定

2*{(制約のない方のモデルの対数尤度) – (制約のある方のモデルの対数尤度)}が、制約の課されているモデルのパラメータの数が自由度となるカイ二乗分布に従う(Wilksの定理)

尤度比検定は便利であるが、モデルによって尤度を評価するのに使うデータ点が異なることがあるので注意が必要。

デビアンス(deviance)。飽和モデル(尤度が可能な最大となるモデル、つまり、データ点ごとにパラメータが異なる値になるモデル)と比べた時の、対数尤度の差の2倍。

c3:離散値データ

<割合を目的変数とするケース>

glmのformula引数に入れられるパターン。

・非集計値。つまり1か0か。

・説明変数ごとの集計値と全体値。

・説明変数ごとの割合(集計値/全体値)と全体値(引数weightsで指定)。

・データの与え方で飽和モデルが異なり、デビアンス変わるので注意。

・分離(separation)のケース。直線回帰では、説明変数にばらつきがないと回帰係数が計算できない。->この話、まだちゃんと理解できてない。

・多項ロジスティック回帰はパッケージnnetのmultinom関数または、multinomRobパッケージを使う。(VGAMパッケージも使えそう。また、この本には書いてないがIIA特性に注意。)

ポアソン回帰>

ポアソン回帰のリンク関数は、対数リンクがよく使われるが、恒等リンクと平方根リンクが適切なケースもある)

c4:過分散

過分散が起きるケース。

対策。

①どれだけ分散が大きいのか推定して、標準偏差の大きさを調整。Fisher情報量の大きさを調節すると考えると良い。

②線形予測子に平均が0で分散をもつ変数を追加して、過剰な分散を表現。ランダム効果。

ランダム効果は生物の個体間の差などデータ点が独立でない状況を表すのに有効。ただし、過少分散は扱えない。

③別の確率分布。ベータ二項分布。

 

GLMMをRでできるパッケージ

http://norimune.net/2365

glmmML, lme4

 

c5:擬似尤度

擬似尤度の拡張。一般化推定方程式(GEE)は擬似尤度の多変量版。

みどりぼんp159の注釈20に「以前は準尤度(quasi likelihood)を使って過分散のあるデータの統計モデリングをしていました。しかし、現在では使用する利点がないので使われなくなっています。」とある。

c6:ランダム効果の変数と混合モデル

Rでは、GLMMはlmer4(lmerTest)の関数glmerのほか、glmmMLやパッケージglmmAKの関数logpoissonRE、パッケージpolytomousの関数polytomous、パッケージsabreRの関数sabreなどで扱える

 

例:ブロック

データ点がどの単位(ブロック、層)に属するかにより、目的変数の値が異なる場合。

説明変数が目的変数に与える影響を見るには、同じブロックのデータ点を比較するべき。

ブロックの効果に対してランダム効果の変数を使うとGLMMの分析になる。

ブロックの効果があるのに、ブロックを無視して分析すると、シンプソンのパラドクス同様、係数の結果が正負逆になるということが起こりうる。

 

glmerの引数nAGQは数値積分の分点の数。ガウス・エルミート積分の一種。ラプラス近似よりも精度がよいことがおおい。1だとラプラス近似と同様。デフォはこれ。

 

マルチレベルモデル。複数のレベルを考える必要性。

第一に、あるレベルでの単位間の分散があるのにそれを無視すると、過分散があるのにないものとして分析した場合のように、目的変数のランダム効果による部分も、説明変数による変動に入ってしまい、説明変数の評価を誤ることがある。

第二に、あるレベルを無視してデータ点をプールすると、シンプソンのパラドクスに陥り、大きく異なる関係を検出してしまう可能性がある。

第三に、レベルによって変数間の関係が異なる場合には、見たいものとは別のところの関係を見誤ってしまうことがある。たとえば、アメリカでの識字率と移民であるかどうかの間には、個人ごとに見ると(1人が1データ点)負の相関があるが、州ごとに集計した識字率と移民の割合で見ると(1つの州が1データ点)正の相関がある、というものがある。

マルチレベルモデルでは、データ点を単純に1つの確率分布とi.i.d.ではうまく表せない時に表現する主な方法の一つ。もう一つの方法は、時系列データなどで使われるように、データ点間の相関を表す相関係数などの行列を与えるもの。

・識別可能性

2つ正規分布する変数の和では、2つの変数の期待値の和や2つの変数の分散の和は、正規分布する1つの変数の場合と同様に推定できる。しかし、2つの変数のそれぞれの期待値やそれぞれの分散は決まらない。このように使用可能なデータから決まらないことを、識別可能性(identifiability)がないという。GLM関係の例は、説明変数が1つで、identityリンク、誤差構造が等分散お正規分布という場合に、説明変数の値に正規分布する測定誤差があると、識別可能性がないと知られている(測定誤差の分散や測定誤差と目的変数の分散の比がわかれば推定できる)。ランダム効果を含む複雑なモデルでは、すぐには識別可能性が判断できないことも考えられるので、注意が必要。

・マルチレベルモデルと階層ベイズ

マルチレベルモデルでレベルの数が多いときなど、多くのランダム効果の変数がある場合には、より複雑な積分が必要。少なくとも現時点ではそれを計算して最尤推定に基づく分析をすることはあまり現実的とは言えない。MCMC法を使ってパラメータ推定する階層ベイズ法による分析では、この積分を避けることができる。

c7:交互作用

E[y] = β1x1 + β2x2 + γx1x2 + β0

のとき、γx1x2が交互作用の項で、γが交互作用項の係数。交互作用と区別して、それぞれの説明変数事態の効果を主効果(main effect)という。交互作用項を含むときの主効果の係数は、他の説明変数が0のときの問題の説明変数の効果の予測値というべきもの。

c8:パラメトリック・ブートストラップ(PB)

尤度比検定は、2つの仮説に対応するモデルでの尤度がそれぞれ求められれば使えるので、広い状況で使用できて強力。しかし、サンプルサイズが大きい必要があることと2つのモデル間には包含関係が必要で仮説が限定されることが弱点。これらはパラメトリック・ブートストラップを使うことで解決できる。

まだ全然よくわかってないので

http://amzn.asia/hHwDCtW

これで追加学習する。

理論はこっちで。

http://amzn.asia/fqqdrBa

 

c9:新しい誤差構造とリンク関数など

 

・tweedieモデル

正あるいは0の値をとる目的変数に対して、リンク関数にはべき乗リンクを使い、誤差構造は分散が期待値のp乗となる確率分布を使う。分散が平均とともに変化するデータは多い。Tweedieモデルでは目的変数の期待値μと線形予測子ηの間にμ^q = ηというべき乗リンクを考える。

・ベータ回帰

・二重一般化線形モデル(dual GLM, double GLM)

目的変数の期待値と線形予測子の間に、リンク関数(目的変数の期待値)=線形予測子、という式で表される関係に加えて、期待値まわりの分散の部分にも線形予測子に類似したものを考え、分散も何らかの変数にともない変化できるようにする。

説明変数の値と目的変数の分散の間にはっきりした関係があるときや目的変数の分散と似た挙動をする変数があるときに有効。

 

・Bradley-Terryモデル

https://www.gavo.t.u-tokyo.ac.jp/~mine/japanese/IT/2017/toukei171211.pdf

 

【学習予定】

 

『ブートストラップ入門 (Rで学ぶデータサイエンス 4)』

http://amzn.asia/hHwDCtW

でブートストラップの理解に励む。

その他、みどりぼん、あひる本、『Rで学ぶベイズ統計学入門』を読んでいきたい。

 

グループごとに通し番号を出てきた順にふりたいときの覚書

【用途】

階層ベイズモデルを組むとき、グループ名に対して、idを振り分けたい。

Rでfactor型に対して、as.numeric()を使うと、振り分けることができるのだけど、

as.numeric()はABC順で、番号を振り分けてしまうみたい。データの上から順番に

出てきたグループ順に数字を振ってくれればうれしい。

for文をゴリゴリ使えば、書けなくはないが、面倒くさい。楽をしたい。

そんなときのための技。

 

【内容】

{forcats}パッケージでカテゴリカル変数(factor型データ)をいじってみる

{forcats}パッケージのfct_inorder()で順番を変更してくれる。

 

統計ソフトRでグループ毎に通し番号を振るにはどうすればよいでしょ... - Yahoo!知恵袋

変更した後、as.numeric()を使えば、解決!!

 

【記憶の検索キーワード】

グループ, 通し番号, ふりたい, 出てきた順, 順番, 階層ベイズ, factor, level