※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 06月 16日
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 = -2`、`b2 = 1` のリンク関数を使って事象の発生確率 `q` を求め、試行回数 `N = 20` の二項分布から値を生成しています。 なので、最尤推定すると、`b1 = -2`、`b2 = 1`、`s = 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期待通りにパラメータが推定されているのが分かります。