R: glmmBUGS

Brown P, Zhou L (2010) MCMC for Generalized Linear Mixed Models with glmmBUGS. R Journal 2(1): 13–17

WinBUGSを使って(呼び出しは手動だが)、GLMMの計算をさせるパッケージ。Gaussian Markov random fieldにも対応している。

例題を自分でもためしてみた。

MASSパッケージのbacteriaデータを使用する。

library(MASS)
data(bacteria)

応答変数を整数型に(N -> 0, Y -> 1)。

bacteria$yInt <- as.integer(bacteria$y == "y")

bacteria$trtの水準名を変更。記号を含まないように。

levels(bacteria$trt) <- c("placebo", "drug", "drugplus")

glmmBUGSパッケージのglmmBUGS()関数を使用する。これによりmodel.bugとgetInits.Rが生成される。

library(glmmBUGS)
bac.rag <- glmmBUGS(yInt ~ trt + week, data = bacteria,
                    effects = "ID",
                    modelFile = "model.bug",
                    family = "bernoulli")

生成されたmodel.bug。

model{

for(DID in 1:NID) {

  RID[DID] ~ dnorm(meanID[DID], TID)
  meanID[DID] <- intercept +inprod2(betaID[] , XID[DID,])

  for(Dobservations in SID[DID]:(SID[DID+1]-1)){

    yInt[Dobservations] ~ dbern(meanobservations[Dobservations])
    logit(meanobservations[Dobservations]) <- RID[DID] + betaobservations * Xobservations[Dobservations]

  }#observations
}#ID


# priors

intercept ~ dflat()
betaID[1] ~ dflat()
betaID[2] ~ dflat()
betaobservations ~ dflat()

TID <- pow(SDID, -2)
SDID ~ dunif(0, 25)

} # model

getInits.Rを読み込み、startingValues変数をセット。

source("getInits.R")
startingValues <- bac.rag$startingValues

WinBUGSを呼び出す。

library(R2WinBUGS)
bac.post <- bugs(bac.rag$ragged, getInits,
                 model.file = "model.bug",
                 n.chain = 3, n.iter = 21000, n.burnin = 1000,
                 parameters.to.save = names(getInits()),
                 n.thin = 20,
                 useWINE = TRUE)

結果を表示する。

bac.params <- restoreParams(bac.post, bac.rag$ragged)
bac.sum <- summaryChain(bac.params)
print(bac.sum$betas)

結果。

                   mean       pval         sd      0.5%      2.5%        5%     50%
observations -0.1609507 0.00100000 0.05485962 -0.305804 -0.268710 -0.252005 -0.1614
trtdrug      -1.4794127 0.02833333 0.82326786 -3.968010 -3.189450 -2.895100 -1.4480
trtdrugplus  -0.9397532 0.11233333 0.82739272 -3.440025 -2.674225 -2.308100 -0.9125
                   95%     97.5%       99.5%
observations -0.070245 -0.053756 -0.01911455
trtdrug      -0.204680  0.052824  0.66198950
trtdrugplus   0.366200  0.621550  1.24301500

ちなみに、glmmMLではこうなる。

library(glmmML)
bac.glmm <- glmmML(yInt ~ trt + week, family = binomial,
                   data = bacteria, cluster = ID,
                   method = "Laplace")
summary(bac.glmm)

結果

Call:  glmmML(formula = yInt ~ trt + week, family = binomial, data = bacteria,      cluster = ID, method = "Laplace") 


               coef se(coef)      z Pr(>|z|)
(Intercept)  3.1440  0.61916  5.078 3.82e-07
trtdrug     -1.3202  0.64199 -2.056 3.97e-02
trtdrugplus -0.7955  0.65172 -1.221 2.22e-01
week        -0.1437  0.05092 -2.822 4.77e-03

Scale parameter in mixing distribution:  1.147 gaussian 
Std. Error:                              0.3782 

Residual deviance: 197.8 on 215 degrees of freedom 	AIC: 207.8 

タグ：R MCMC glmm