データ解析のための統計モデリング入門 一般化線形混合モデル(GLMM) 読書メモ4


2017年 05月 26日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。
今回は第7章、一般化線形混合モデル(GLMM)についてのまとめの四回目です。

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。
まずは二項分布と正規分布を混ぜた混合分布について解説されています。
なので、混合分布をグラフにしてみたいところですが、やっぱり変数がたくさんあって複雑なので、まずは感覚をつかむために数値を羅列してくれるコードを用意しました。
コードはRで書きました。

linkFunction <- function(b1, b2, x) {
  function(r) {
    1 / (1 + exp(-b1 - b2 * x - r))
  }
}

l <- linkFunction(-0.5, 0.3, 2)

cat("r\ty=0\t\ty=1\t\ty=2\t\ty=3\n")
sum <- vector(length = 4)
for (r in seq(-5, 5, by = 0.5)) {
  q <- l(r)
  cat(r, "\t")
  for (y in 0:3) {
    b <- dbinom(y, 3, q)
    n <- dnorm(r, mean = 0, sd = 2)
    mixture <- b * n
    sum[y + 1] <- sum[y + 1] + mixture * 0.5
    cat(mixture, "\t")
  }
  cat("\n")
}
cat("sum\t", sum[1], "\t", sum[2], "\t", sum[3], "\t", sum[4], "\n")

リンク関数は xr を変数に含んでいます。
x2 に固定した場合について計算しました。
なので、あとは r を与えれば、二項分布している事象の発生確率 q が求まります。
r は正規分布する変数です。
今回はr は平均 0、標準偏差 2 の正規分布から得られるとしました。
事象の試行回数は 3 にして、0 回から 3 回事象が発生する場合の確率を計算して出力しています。
r-5 から 5 まで、0.5 刻みで計算しました。
コードを実行すると、次の出力が得られます。

r        y=0             y=1             y=2             y=3
-5         0.008571241     0.0001914794    1.425867e-06    3.539279e-09    
-4.5     0.01529937      0.0005635068    6.918364e-06    2.831304e-08    
-4         0.02542036      0.00154367      3.124683e-05    2.108318e-07    
-3.5     0.03909264      0.003913947     0.0001306212    1.453086e-06    
-3         0.05514584      0.009102904     0.0005008711    9.186514e-06    
-2.5     0.07038014      0.01915423      0.001737632     5.254481e-05    
-2         0.07963943      0.03573468      0.005344786     0.0002664708    
-1.5     0.07772425      0.05749969      0.01417925      0.00116552  
-1         0.06325714      0.0771553       0.031369        0.004251228     
-0.5     0.04148674      0.08342818      0.05592358      0.01249557  
0         0.02138051      0.07088734      0.07834263      0.02886066  
0.5     0.008601664     0.04701976      0.08567559      0.05203704  
1         0.002741934     0.02471168      0.07423798      0.07434107  
1.5     0.0007137071    0.01060504      0.05252712      0.08672284  
2         0.0001570974    0.003848653     0.03142876      0.08555086  
2.5     3.018185e-05    0.001219082     0.0164134       0.07366188  
3         5.187402e-06    0.0003454491    0.007668262     0.0567399   
3.5     8.116416e-07    8.911395e-05    0.003261413     0.03978732  
4         1.169644e-07    2.1173e-05      0.001277585     0.02569661  
4.5     1.564146e-08    4.66824e-06     0.0004644167    0.01540073  
5         1.950226e-09    9.596395e-07    0.0001574019    0.008605787     
sum         0.2548242       0.2235203       0.2303399      0.2828235 

r の増加に従って q も増加するので、r が小さいうちは事象の発生回数は少ない方が確率が高く、r が大きくなると多い方が確率が高くなるのが分かります。
混合分布では、r についての積分がとられるので、最後の行にざっくりとした積分の近似値を出力しておきました。
二項分布と正規分布の混合分布は、だいたいこういう数字になります。