Taglibro de H

Kuboさんのところや足軽日記さんのところで取り上げられていたMCMCglmmを試してみた。


例によって架空のデータを用意。

set.seed(11)

i.logit <- function(x) {
	exp(x)/(1 + exp(x))
}

n <- 21
block <- as.factor(1:8)
m <- length(block)
x <- seq(-1, 1, length = n)
ranef.sd <- 0.2
ranef <- rnorm(m, 0, ranef.sd)

data <- expand.grid(x = x, block = block)
data$y <- rbinom(n * m, 1,
                 i.logit(0.7 + 2 * data$x + ranef[data$block]))

まずはglmmMLで。

## glmmML
library(glmmML)
m.ml <- glmmML(y ~ x, family=binomial, data=data, cluster = block)
print(m.ml)

結果

Call:  glmmML(formula = y ~ x, family = binomial, data = data, cluster = block) 


              coef se(coef)     z Pr(>|z|)
(Intercept) 0.7693   0.2250 3.420 6.27e-04
x           3.1594   0.4676 6.756 1.42e-11

Scale parameter in mixing distribution:  1.032e-07 gaussian 
Std. Error:                              0.3634 

Residual deviance: 144 on 165 degrees of freedom 	AIC: 150 

では、MCMCglmmで。

## MCMCglmm
library(MCMCglmm)

m.mcmc <- MCMCglmm(y ~ x, random = ~block, family="categorical",
                   data = data,
                   prior = list(R = list(n = 1, V = 1),
                                G = list(G1 = list(n = 1, V = 1))),
               nitt = 20000, thin = 10, burnin = 10000)
plot(m.mcmc$Sol)
summary(m.mcmc$Sol)

結果

Error :  外部関数の呼び出し(引数 1) 中に NA/NaN/Inf があります 
 追加情報:  Warning message:
step size truncated due to divergence 
Error in glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,  : 
  inner loop 1; cannot correct step size
 追加情報:  Warning message:
step size truncated due to divergence 

Iterations = 1:1000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean     SD Naive SE Time-series SE
(Intercept)  29.44  29.18   0.9226             NA
x           121.83 112.42   3.5552             NA

2. Quantiles for each variable:

              2.5%    25%   50%    75%  97.5%
(Intercept)  2.309  4.743 15.10  52.45  94.68
x           12.516 19.943 62.50 230.54 327.38

なんか変な値になっている。事前分布の与え方とかがよくわかっていないので(デフォルトだとこの場合エラーになる)、そのあたりで間違っている可能性もあるのだが。

ちなみにJAGSでは。

## rjags
library(rjags)
model <- jags.model("MCMCglmm_test.bug",
                    data = list(N = n * m, M = m,
                                y = data$y,
                                x = data$x,
                                b = data$block),
                    nchain = 3, n.adapt = 5000)
post <- coda.samples(model, c("beta", "beta.x", "sigma"),
                        n.iter = 10000, thin = 10)
plot(post)
summary(post)

MCMCglmm_test.bug

var
  N,	# number of data
  M,	# number of blocks
  x[N], y[N], p[N], b[N], e[M],
  beta, beta.x,	# parameters
  tau, sigma;	# hyperparameters

model {
  for (i in 1:N) {
    y[i] ~ dbern(p[i]);
    logit(p[i]) <- beta + beta.x * x[i] + e[b[i]];
  }
  for (j in 1:M) {
    e[j] ~ dnorm(0, tau);
  }

  # priors
  beta ~ dnorm(0.0, 1.0E-3);
  beta.x ~ dnorm(0.0, 1.0E-3);
  tau ~ dgamma(1.0E-2, 1.0E-2);
  sigma <- 1/sqrt(tau);
}

結果

Iterations = 5010:15000
Thinning interval = 10 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean     SD Naive SE Time-series SE
beta   0.7974 0.2660 0.004856       0.004881
beta.x 3.2919 0.4858 0.008870       0.009512
sigma  0.3000 0.2162 0.003947       0.005279

2. Quantiles for each variable:

          2.5%    25%    50%    75%  97.5%
beta   0.29698 0.6209 0.7857 0.9696 1.3427
beta.x 2.38036 2.9616 3.2678 3.6035 4.3134
sigma  0.07262 0.1503 0.2394 0.3866 0.8755