R: MCMCglmm 二項分布

応答変数が0/1, y/n, あり/なしのデータ。MCMCglmmではfamily="categorical"とすればよいようだが。

# [R] MCMCglmm でポアソン回帰 | singular pointという記事があったので、とりあえずメモしておく。

データとして、MASSのbacteriaをつかってみる。

library(MASS)
data(bacteria)

lmer()では

library(lme4)
bac.lmer <- lmer(y ~ week + trt + (1|ID), family = binomial,
                 data = bacteria)
print(bac.lmer)

結果

Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ week + trt + (1 | ID) 
   Data: bacteria 
   AIC   BIC logLik deviance
 207.8 224.7 -98.89    197.8
Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 1.3146   1.1465  
Number of obs: 220, groups: ID, 50

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.1440     0.5290   5.943 2.79e-09 ***
week         -0.1437     0.0483  -2.975  0.00293 ** 
trtdrug      -1.3202     0.6252  -2.111  0.03473 *  
trtdrug+     -0.7956     0.6401  -1.243  0.21395    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
         (Intr) week   trtdrg
week     -0.552              
trtdrug  -0.628  0.071       
trtdrug+ -0.596  0.039  0.489

glmmML()では

library(glmmML)
bac.ml <- glmmML(I(y == "y") ~ week + trt, cluster = ID, family=binomial,
                 data = bacteria)
print(bac.ml)

結果

Call:  glmmML(formula = I(y == "y") ~ week + trt, family = binomial,      data = bacteria, cluster = ID) 


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

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

        LR p-value for H_0: sigma = 0:  0.007011 

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

MCMCglmmだと、

library(MCMCglmm)
bac.mcmc <- MCMCglmm(y ~ week + trt, random = ~ID, family = "categorical",
                     data = bacteria,
                     prior = list(R = list(V = 1, nu = 0),
                                  G = list(G1 = list(V = 1, nu = 0))),
                     nitt = 60000, thin = 50, burnin = 10000,
                     verbose = FALSE)
summary(bac.mcmc)

結果

 Iterations = 10001:59951
 Thinning interval  = 50
 Sample size  = 1000 

 DIC: 5.125227 

 G-structure:  ~ID

   post.mean  l-95% CI u-95% CI eff.samp
ID      5145 1.448e-12    17051    10.68

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     21231     3293    39949    129.4

 Location effects: y ~ week + trt 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   231.131  121.768  355.445    220.3 <0.001 ***
week          -10.284  -18.973   -2.641    501.6  0.006 ** 
trtdrug       -94.833 -187.894   -7.489    586.1  0.030 *  
trtdrug+      -54.059 -146.717   40.176   1000.0  0.206    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

何だかスケールがちがうような気がするので、priorのところで分散を1に固定してみる。

bac.mcmc2 <- MCMCglmm(y ~ week + trt, random = ~ID, family = "categorical",
                      data = bacteria,
                      prior = list(R = list(V = 1, nu = 1, fix = 1),
                                   G = list(G1 = list(V = 1, nu = 0))),
                      nitt = 60000, thin = 50, burnin = 10000,
                      verbose = FALSE)
summary(bac.mcmc2)

結果

 Iterations = 10001:59951
 Thinning interval  = 50
 Sample size  = 1000 

 DIC: 192.108 

 G-structure:  ~ID

   post.mean l-95% CI u-95% CI eff.samp
ID     2.874 0.002028    7.083    93.78

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: y ~ week + trt 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   3.90023  2.32113  5.55957    117.7 <0.001 ***
week         -0.17786 -0.31257 -0.06714    325.2  0.002 ** 
trtdrug      -1.62696 -3.40207 -0.06575    291.0  0.044 *  
trtdrug+     -0.93034 -2.56224  0.94574    501.9  0.276    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

どうするのが正解なのかはよくわからない。

タグ：R MCMCglmm