R: glmmADMBによる負の二項分布のGLMM

glmmADMBをつかって、負の二項分布のGLMMへのあてはめをやってみる。

模擬データを生成して、それを解析させてみる。以下のようなコード。

library(glmmADMB)

## 模擬データの生成
set.seed(123)
n.data <- 100
n.block <- 10
## ランダム効果の標準偏差は0.5
ranef <- rep(rnorm(n.block, 0, 0.5), each = n.data / n.block)
## サイズパラメーターを5、平均をexp(2 + ランダム効果)と指定
y <- rnbinom(n.data, size = 5, mu = exp(2 + ranef))

## 平均と分散
mean(y)
var(y)

## ヒストグラム
h <- hist(y, breaks = 0:max(y), right = FALSE, plot = FALSE)
barplot(h$count, names = h$breaks[-length(h$breaks)],
        xlab = "Y", ylab = "Frequency", las = 1)

## データフレーム
data <- data.frame(y = y, 
                    b = as.factor(rep(1:n.block,
                                      each = n.data / n.block)))

## モデルあてはめ
fit1 <- glmmadmb(y ~ 1, random = ~ 1|b, family = "nbinom",
                data = data)
summary(fit1)

出力は以下のようになった。

平均と分散

> mean(y)
[1] 8.73
> var(y)
[1] 50.09808

ヒストグラム

結果のサマリー

Call:
glmmadmb(formula = y ~ 1, data = data, family = "nbinom", random = ~1 | 
    b)

AIC: 603.7 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.043      0.165    12.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=100, b=10 
Random effect variance(s):
Group=b
            Variance StdDev
(Intercept)   0.2336 0.4833

Negative binomial dispersion parameter: 4.2989 (std. err.: 1.0671)

Log-likelihood: -298.85 

タグ：glmm