So-net無料ブログ作成

累積ロジットモデルのベイズ推定(2) [統計]

前回のつづき。今度はランダム効果の入ったデータを推定してみる。

まずはデータを用意。

# inverse logit
invlogit <- function(x) exp(x)/(1+exp(x))

set.seed(12345)
n.rank <- 5
plot <- 1:20
n.plot <- length(plot)
time <- 1:5
n.time <- length(time)
n <- n.plot * n.time
data <- expand.grid(plot = plot, time = time)
slope <- rgamma(n.plot, shape = 15^2/5^2, rate = 15/5^2)
data$slope <- rep(slope, n.time)
ranef.plot <- rnorm(n.plot, 0, 1)
ranef.time <- rnorm(n.time, 0, 1)
logit.c1 <- 7.5 - 0.5 * slope +
            ranef.plot[data$plot] + ranef.time[data$time]
data$c1 <- rbinom(n, n.rank - 1, invlogit(logit.c1)) + 1

累積データを作成して、JAGSにかける。

cum <- matrix(NA, nrow = n.plot * n.time, ncol = n.rank - 1)
cum[,1] <- as.numeric(data$c1 >= 2)
cum[,2] <- as.numeric(data$c1 >= 3)
cum[,3] <- as.numeric(data$c1 >= 4)
cum[,4] <- as.numeric(data$c1 >= 5)

library(rjags)
model <- jags.model("cumlogit_MCMC2.bug",
                    data = list(N = n, M = n.rank - 1,
                                K = n.plot, L = n.time,
                                s = data$slope, c = cum,
                                plot = data$plot, time = data$time),
                    nchain = 3, n.adapt = 10000)
post <- coda.samples(model,
                     variable = c("b", "beta", "sigma.p", "sigma.t"),
                     n.iter = 20000, thin = 20)

BUGSコード

var
  N,  # number of data
  M,  # number of ranks
  K,  # number of plots
  L,  # number of mesurements
  beta,
  b[M],
  c[N, M],
  s[N],
  plot[N], e.plot[K],
  time[N], e.time[L];

model {
  for (j in 1:M) {
    for (i in 1:N) {
      c[i, j] ~ dbern(p[i, j]);
      logit(p[i, j]) <- b[j] + beta * s[i] +
                        e.plot[plot[i]] + e.time[time[i]];
    }
    b[j] ~ dnorm(0.0, 1.0E-4);
  }

  for (i in 1:K) {
    e.plot[i] ~ dnorm(0, tau.p);
  }

  for (i in 1:L) {
    e.time[i] ~ dnorm(0, tau.t);
  }

  beta ~ dnorm(0.0, 1.0E-4);
  tau.p ~ dgamma(1.0E-3, 1.0E-3);
  tau.t ~ dgamma(1.0E-3, 1.0E-3);
  sigma.p <- 1/sqrt(tau.p);
  sigma.t <- 1/sqrt(tau.t);
}

結果

> summary(post)

Iterations = 10020:30000
Thinning interval = 20 
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
b[1]    14.4082 2.6237 0.047902       0.195126
b[2]    11.6440 2.4110 0.044018       0.180754
b[3]    10.2122 2.3343 0.042619       0.172398
b[4]     7.4226 2.1665 0.039554       0.152169
beta    -0.7306 0.1376 0.002513       0.009946
sigma.p  2.0650 0.5855 0.010690       0.019606
sigma.t  1.6044 0.8677 0.015841       0.021022

2. Quantiles for each variable:

           2.5%     25%    50%     75%   97.5%
b[1]     9.8577 12.6381 14.235 15.9341 19.9900
b[2]     7.4117 10.0158 11.449 13.0504 16.8387
b[3]     6.1724  8.5996 10.039 11.5374 15.2372
b[4]     3.5728  5.9820  7.287  8.6834 12.0245
beta    -1.0236 -0.8163 -0.720 -0.6385 -0.4859
sigma.p  1.1528  1.6539  1.988  2.3900  3.4480
sigma.t  0.6625  1.0580  1.399  1.8622  3.8623

同じデータをdrmで解析してみる。

library(drm)
fit <- drm(c1 ~ slope + cluster(plot) + Time(time),
                data = data,
                start = c(0, 1, 2, 3, -1),
                print.level = 0)

結果

> summary(fit)
Call:
drm(formula = c1 ~ slope + cluster(plot) + Time(time), data = data, 
    start = c(0, 1, 2, 3, -1), print.level = 0)

Coefficients:
                 Value Std. Error   z value
(Intercept)1 -8.707243 1.09793339 -7.930575
(Intercept)2 -7.003733 0.94661067 -7.398747
(Intercept)3 -6.134059 0.88513704 -6.930067
(Intercept)4 -4.475916 0.79272635 -5.646230
slope        -0.444208 0.05900017 -7.528927

Residual deviance: 216.881 	AIC: 226.881 

Convergence code 1 in 32 iterations (See help(nlm) for details)

Correlation of Coefficients:
             (Intercept)1 (Intercept)2 (Intercept)3 (Intercept)4
(Intercept)2        0.948                                       
(Intercept)3        0.933        0.966                          
(Intercept)4        0.892        0.909        0.924             
slope               0.957        0.957        0.953        0.921

なんだか微妙に結果が違うが、何が原因なのだろう? 推定法の違いによるものなのか、どっか間違えているのか。切片の符号が逆なのは、モデルでの符号の与え方の違いによるものなので、それは問題ない。

データに期待値の線を重ねてプロットしてみる。

Rplot1.png
ベイズ推定によるもの

Rplot2.png
drmによるもの


nice!(0)  コメント(0)  トラックバック(1) 
共通テーマ:学問

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1