霞と側杖を食らう

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

VSCodeやGIt・Github、Markdownの整備の覚書

【用途】 

エディタの選択、EmacsとかAtomとかあるけど、VSCode使いやすそうではということでVSCode整備。今のところ、pythonmarkdownとかが書ければいいかなくらい。RはRStudioで良いし。このへんを整備するが、 また整備する機会のために個人的メモとして、記録しとく。あと、整備で詰まったところとか書き残してる。

 

【内容】 

VSCode

元からsurfaceに入ってたので、それを使用。

コマンドパレットは Ctrl + Shift + P

 

python

python anaconda3

www.messyer813.com

 

分かりやすい。

 

Jupyter をVSCode上で展開したくて、後程、拡張でそれに関係するものを入れる。

 

<Git・Github

Gitが入ってなかったのと、Githubのアカウントは昔作ったまま放置してた。

これに従って入れて、整備していくことにした。

employment.en-japan.com

 

vimのcondfigの終わらせ方が分からなくて一瞬戸惑った。

qiita.com

 

hacknote.jp

 

qiita.com

 

リポジトリのクローンで

Permission denied (publickey)ってエラーが出て、困った。

hacknote.jp

基本的にこれに従えばよいのだが、「`ssh-agent`」を「"ssh-agent"」と打っていたせいで、上手くいかなかった。"と`は別物なのね。。。

 

コミットしようとしたら、who are you?って聞かれた

qiita.com

これでなんとかなった。

コミットでコメント残したと思ったら、aborting commit due to empty commit messageと言われる。

qiita.com

どうやらこういうことらしい。あと#の後にスペースが要る。

  

<拡張>

このへんから欲しい拡張機能を選んでインストールした

ics.media

qiita.com

 

vscode-icons

アイコンついて見やすい。

・Clipboard Ring

コピーを複数記憶できる。

 ・Bracket Pair Colorizer

カッコが色付きになって判別しやすくなる。

・Path Autocomplete

パスの入力が楽に。

・Rainbow CSV

csvファイルを見るときに列で色別になって見やすくなる。

・Partial Diff

選択範囲の差分を調べられる。写経のときにスペルミスとかで溶かしてた時間を削減できるはず。

・Trailing Spaces

スペースを強調してくれる。変なスペースで溶かしてきた時間ともこれでお別れ。

 

Markdown

 マークダウン関連の拡張

 ・Markdown All in One

mdファイルを開いた際に自動的にプレビューを起動(デフォルトでは無効)。

ショートカットキーでMarkdownのタグを入力可能(ctrl+bで太字、ctrl+shift+]でインデントの追加など)。

・Auto-Open Markdown Preview

自動でプレビュー起動。

Markdown PDF

mdからpdfにエクスポート。mdファイル保存して開いて右クリックで選択。

これの問題点として、数式入っているmdファイルをpdfにしても、数式を変換してくれない。なので、数式使うときは下のMarkdown Preview Enhancedを使用する。

Markdown Preview Enhanced

プレビューがつよくなってる。数式入りMarkdownもpdfにできる。直接吐き出す方法は分からないけれど、ブラウザで開く(open in browser)で、chromeで開き、印刷でpdfで保存とすれば出来るみたい。

Markdown + Math

数式入れるため。$$内に数式を書いてctrl+shift+.と打って、数式対応のプレビューが表示。

qiita.com

次に、数式を複数行書くとき、=を揃えて書こうとしてつまづいた。

\begin{aligned}

a &=b\\

&=c

\end{aligned}

で、できた。KaTeXだと、{align}じゃなくて、{aligned}だった。

tera-kun.hatenadiary.com

・Paste Images

 画像をコピーしたら、Ctrl+Alt+V でMarkdownファイル内に貼り付けられ、同一ディレクトリに保存。(保存先は設定ファイルで変更可能。)

・Jupyter

Jupyterを使うため。

VS Code Jupyter Notebook Previewer

プレビューに使うため。

・ Tables Generator

excelの表を変換してくれるサービス。markdownで表作るのってめんどくさいけど、これがあれば楽になるはず。

www.procrasist.com

 

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

VSCode, GIt, GIthub, Markdown, 拡張, 整備, エディタ

 

Rによる多重代入法の学習記録

【学習動機】

Kaggleで欠測値を埋めなきゃということがある。他にも埋めなきゃいけない状況がある。
欠測値の理論は、因果推論の一般的な枠組みとして捉えるのに役立つ。
あと、ベイズと関わってくるので、勉強しておきたいと思った。実際に、Rで実行できるように。

【学習記録】

高橋・渡辺『欠測データ処理 Rによる単一代入法と多重代入法』 https://amzn.to/3WxKJKN

で学んでいくこととする。

<総括>

Rで実際にどうやって欠測値に多重代入法を用いて対処するかが分かったのは、とても良い。
今出ている日本語書籍では、理論を学ぶことはできるが、どれをどうやって使うのかがすぐに分からなくて困っていた。
データセットは書籍HPに用意すべき。データは書面に記載しているってのはモダンじゃない。
norm2とmiceとAmelia。代入・分析・統合。

c1:Rによるデータ解析

Rの使い方と欠測あるときの反応を少しチェックしている。

c2:不完全データの統計解析

多くの統計ソフトウェアでは、リストワイズ除去で、欠測値を含む観測対象の行をすべて削除し、擬似的「完全」な状態にして分析することが多い。
欠測パターンと欠測メカニズムで分類される。
欠測パターンには、単変量欠測パターン、単調欠測パターン、計画的欠測パターン、一般的な欠測パターン。多重代入法では、いずれの欠測パターンにも有効で、欠測パターンはそんなに重視されない。
欠測メカニズムには、完全に無作為な欠測(MCAR:Missing Completely At Random)と条件付きで無作為な欠測(MAR:Missing At Random)と無作為ではない欠測(NMAR:Not Missing At Random)がある。 Dをn×pのデータセット D = {  D _ {obs}, D _ {miss} }。Kは回答指示行列でn×p行列。 K _ {i j}=1で観測、  K _ {ij}=0で欠測
・MCAR
 Pr(K|D) = Pr(K) であり、欠測する確率が、データと独立。
・MAR
 Pr(K|D) = Pr(K|D _ {obs}) で、観測データを条件とし、欠測確率の分布が非観測データから独立。欠測データがMARならば、欠測を除去する分析は偏っている。この偏りは共変量(covariate)または補助変数(auxiliary variable)による代入で是正可能。 ・NMAR
 Pr(K|D) \neq Pr(K|D _ {obs}) で、欠測する確率が欠測した変数の値自体に依存している。代入法で偏りを是正できるとは限らず、選択モデルやパターン混合モデルを用いた分析が必要で、感度分析が有用。

分析の目的が母集団の推定なら、多重代入法を使用すべき。

c3:単一代入法

従来、公的統計では、調査データの集計が主目的で、確定的単一代入法を用いることが通例。方法として、よく用いられるのが、回帰代入法、比率代入法、平均値代入法、ホットデック法。
・確定的単一代入法はガウス・マルコフの仮定あるとよい。
・比率代入法は不均一分散のとき使える。
・平均値代入法は百害あって一利なし。
・ホットデック法(hot deck imputation)は、ノンパラな代入法。ホットデックでは、欠測している人の補助変数の情報を利用して、回答者から似通った値のものを選び出し、その回答者の値を代入値として採用。何を基準にするかは議論あり。最近隣法が多い。
HotDeckImputationパッケージのimpute.NN_HD()で実行できる。
・確率的回帰代入法
各々の代入値に誤差項を加える。誤差項は、平均が0、分散が回帰モデルの残差分散の正規分布から無作為抽出するなどのやり方があり、miceのmice()で引数methをnorm.nobとすれば実行できる。

c4:多重代入法の概要

ガウス・マルコフの仮定とMARの過程がすべて満たされた場合、多重代入法によるパラメータ推定量はBLUEになる。
代入、分析、統合の3ステージ
統合方法。  \tilde{\theta_m}をパラメータ \thetaのm番目の代入済みデータセットに基づいた推定値とする。
統合した点推定値は
 \bar{\theta_M} = \frac{1}{M} \Sigma \tilde{\theta_m}
多重代入法の代入値の分散は代入間分散(between)と代入内分散(within)に分解でき、
代入内分散は、
 \bar {W _ m}  = \frac{1}{M}  \Sigma var(\tilde{\theta_m})
代入間分散は、
 \bar{B _ M} = \frac{1}{M-1} \Sigma (\tilde {\theta _ m} -  \bar{\theta_M}) ^2
 \bar{\theta_M}の分散 T_Mは、これらを合算して
 T_M = \bar{W_m} + (1 + \frac{1}{M} ) \bar{B_M}

この方法をRubin(1987)のルールという。

代入間の繰り返しをどれだけ行えば収束するかは、様々な要素が絡んでくるが、とりわけ欠測情報の比率が重要。
適切な多重代入法(proper multiple imputation)でないと、推測の結果が妥当でなくなるおそれあり。
代入モデルと分析モデルが同一の変数を持ち、同じ数のパラメータを推定するとき、2つのモデルには適合性(congeniality)があるという。適合していない場合、多重代入法のパラメータ推定量の一致性は保証されない。
代入モデルよりも大きな分析モデルを使うと、リストワイズ除去より悪い結果になりうる。
また、多重代入済みデータ数は多い方が望ましいとされる。

c5:多重代入法のアルゴリズム

3つのアルゴリズム
・DAアルゴリズム。norm2でできる。しかし、norm2が消えてる。
FCSアルゴリズム。miceでできる。これが主流か?
・EMBアルゴリズム。Ameliaでできる。 (EMアルゴリズムについてはまた別の機会で詳しく学ぶ)

norm2パッケージがCRANから消えてるので

url <- "https://cran.r-project.org/src/contrib/Archive/norm2/norm2_2.0.1.tar.gz"
pkgFile <- "norm2_2.0.1.tar.gz"
download.file(url = url, destfile = pkgFile)
install.packages(pkgs=pkgFile, type="source", repos=NULL)
unlink(pkgFile)

あと、ホットデックのコードを動かした後、Ameliaパッケージの多重代入するとエラー起きるかもしれないらしい。
Amelia、norm2、miceと見比べて、すべての変数が量的な連続変数なら、どれも変わらず、miceが時間的に効率悪いか。
質的データや時系列データになると、アルゴリズムの使い分けが必要。

c6:多重代入モデルの診断

MARが正しいと仮定したとき、代入の結果がMARの仮定と整合性があるかを診断する。
MCARの場合を除くと、観測データの分布と代入データの分布は一致せず、むしろ、異なっていることが期待される。診断で、分布の違いを発見し、MARの仮定と整合性があるか確認する。分布が極端に違う場合も、どの変数に問題があるか見つけるのに診断ツールが役立つ。 Ameliaとmiceの診断ツールを用いる。

c7,8:量的データの多重代入法

・多重代入法を用いた平均値のt検定 AmeliaとMKmiscを使うか、miceのpool.scalar()の返り値から直接計算。
欠測データによって、自由度の計算が大変。分散の比率λと代入の不確実性の分散の相対的な増加rを計算する。

このへん、本のコードが変。本が書かれた後にアップデートされて関数が変わったのか知らないが。 miceのpool.scalar()は引数にmethodはとらないし、返り値にlambdaはない。
rは返ってくるので、 \lambda = r/(1+r) で自分で計算しないとかと。

・重回帰分析
回帰診断で、ジャック・ベラの正規性検定(JB検定)で誤差項の正規性を検証、ブルーシュ・ペイガン検定(BP検定)で不均一分散の検証、分散拡大要因(VIF:Variance Inflation Factor)で多重共線性の検証、ボンフェロー二外れ値検定で外れ値の検証。

c9,10:質的データの多重代入法

・ダミー変数のある重回帰分析 質的データの代入法に関して、決定打はないようだ。ここでは、応用研究で評判が高いとされるmiceによる方法とhot.deckの多重ホットデックによる方法が紹介されている。

・ロジスティック回帰
miceでいつも通り代入して、glm.mids()の引数familiyをbinomialとすればよさげ。

c11:時系列データの多重代入法

一般的に、時系列データにおける欠測値は、補完(interpolation)やカルマンフィルタを用いて処理されることが多い。
一方で、多重代入法を用いて欠測値を処理することで不確実性を反映した分析を可能にする。

c12:パネルデータの多重代入法

Ameliaで実行可能。amelia()の引数csに横断面の単位、引数tsに時系列の単位、引数polytimeに3以内の多項式を指定、引数lagsにラグ付き従属変数を指定。未来の値から過去の値を予測することで代入モデルの精度を向上させたいなら引数leadsを指定してもよい。
Ameliaに入っているafricaのデータを用いて、plmでパネルデータ分析をする。
p151のコードを回してみたが、変量効果モデルの推定がうまくいかない。特異行列が出て、solve()で逆行列を求めるところがうまくいってないみたいだが。。。

c13:感度分析:NMARの統計分析

SensMice(SensMiceDA)パッケージを使えば、感度分析は可能。感度分析の目的は、欠測データに関する設定を変えて、分析結果にどう影響するかを調べること。デルタ修正のパラメータデルタをどう設定するかが重要。デルタの設定方法として、変数の標準偏差を基準とする方法がある。標準偏差の-1倍, -0.5倍, 0倍, 0.5倍, 1倍した5つのデルタ値を用意して、どう影響するかを見て、ドメイン知識と照合して検討。

c14:事前分布の導入

Ameliaならamelia()の引数boundsに(変数の位置、最小値、最大値)を行にした行列を指定することができる。リッジ事前分布を使ってモデルの安定性を実現したいなら、amelia()の引数empriで指定。
norm2ならば、emNorm()のprior="ridge",prior.df=0.5などとする。miceはsqueeze()を用いる。

【学習予定】

高井・星野・野間『欠測データの統計科学』

https://www.amazon.co.jp/%E6%AC%A0%E6%B8%AC%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E7%B5%B1%E8%A8%88%E7%A7%91%E5%AD%A6%E2%80%95%E2%80%95%E5%8C%BB%E5%AD%A6%E3%81%A8%E7%A4%BE%E4%BC%9A%E7%A7%91%E5%AD%A6%E3%81%B8%E3%81%AE%E5%BF%9C%E7%94%A8-%E8%AA%BF%E6%9F%BB%E8%A6%B3%E5%AF%9F%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E3%81%AE%E5%AE%9F%E9%9A%9B-%E7%AC%AC1%E5%B7%BB-%E6%98%9F%E9%87%8E-%E5%B4%87%E5%AE%8F/dp/400029847Xwww.amazon.co.jp

阿部『欠測データの統計解析』

https://www.amazon.co.jp/%E6%AC%A0%E6%B8%AC%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E7%B5%B1%E8%A8%88%E8%A7%A3%E6%9E%90-%E7%B5%B1%E8%A8%88%E8%A7%A3%E6%9E%90%E3%82%B9%E3%82%BF%E3%83%B3%E3%83%80%E3%83%BC%E3%83%89-%E9%98%BF%E9%83%A8-%E8%B2%B4%E8%A1%8C/dp/4254128592/ref=pd_lpo_sbs_14_t_2?_encoding=UTF8&psc=1&refRID=3NQKK07ZT4K6GGH8AC3Ewww.amazon.co.jp

あたりで理論を補強していきたい。

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

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

長島先生が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://www.amazon.co.jp/%E8%A8%88%E9%87%8F%E7%B5%8C%E6%B8%88%E5%AD%A6-y21-%E6%B5%85%E9%87%8E-%E7%9A%99/dp/4641163367
の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

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