TIM Labs

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

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、サンプルデータを用意して、実際に一般化線形混合モデルを最尤推定してみるコードを用意しました。 コードはRで書きました。

library('glmmML')

dataSize <- 100
N <- 20

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

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 2)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  id <- 1:dataSize
  data.frame(x, y, id)
}

data <- createData()
glmmML(cbind(y, N - y) ~ x, data = data, family = binomial, cluster = id)

コードを実行するには、glmmMLパッケージをインストールする必要があります。

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータは、0 から 5 までの一様分布から x を生成し、標準偏差 s = 2 の正規分布から r を生成し、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、最尤推定すると、b1 = -2b2 = 1s = 2 に近い値を推定してくれるはずです。

コードを実行すると以下のような出力が出力されます。

Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = data,      cluster = id)


              coef se(coef)      z Pr(>|z|)
(Intercept) -2.163   0.4570 -4.732 2.22e-06
x            1.210   0.1604  7.539 4.73e-14

Scale parameter in mixing distribution:  2.144 gaussian
Std. Error:                              0.1654

        LR p-value for H_0: sigma = 0:  5.114e-131

Residual deviance: 311.6 on 97 degrees of freedom       AIC: 317.6

期待通りにパラメータが推定されているのが分かります。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/646

コメントする

このブログ記事について

このページは、sato.yuichiroが2017年6月16日 12:00に書いたブログ記事です。

ひとつ前のブログ記事は「書評:『自然言語処理と深層学習《C言語によるシミュレーション》』」です。

次のブログ記事は「ChainerがMicrosoft Azure, Windows に対応」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。