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

【学習動機】

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

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で学ぶベイズ統計学入門』を読んでいきたい。