Stanでordered_logistic

Stanではordered_logisticという分布がつかえるので、ためしてみた。

以前drmを使ってみたときのデータを使用。

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

set.seed(123)
plot <- 1:30
n.plot <- length(plot)
time <- 1:3
n.time <- length(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.plot * n.time, 4, invlogit(logit.c1)) + 1

Stanのコード。Stanのマニュアルを参考にした。

model.text <- "
data {
  int<lower=2> K;
  int<lower=0> N;
  int<lower=0> Np;
  int<lower=0> Nt;
  int<lower=1,upper=K> y[N];
  int<lower=1,upper=Np> p[N];
  int<lower=1,upper=Nt> t[N];
  real x[N];
}
parameters {
  real beta;
  real ep[Np];
  real et[Nt];
  ordered[K-1] c;
  real<lower=0> sigma[2];
}
model {
  for (i in 1:N) {
  	y[i] ~ ordered_logistic(x[i] * beta + ep[p[i]] + et[t[i]], c);
  }
  ep ~ normal(0, sigma[1]);
  et ~ normal(0, sigma[2]);
  beta ~ normal(0, 100);
  sigma[1] ~ uniform(0, 100);
  sigma[2] ~ uniform(0, 100);
}
"

RStanで実行。

## Stan
library(rstan)

model <- stan_model(model_code = model.text)
fit <- sampling(model,
                data = list(K = 5, N = nrow(data), Np = n.plot, Nt = n.time,
                            y = data$c1, p = data$plot, t = data$time,
                            x = data$slope),
                chains = 4, iter = 2000, warmup = 1000)

結果

> print(fit, digits = 2)
Inference for Stan model: model.text.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

            mean se_mean   sd    2.5%     25%     50%     75%  97.5% n_eff Rhat
beta       -0.90    0.01 0.18   -1.29   -1.02   -0.89   -0.78  -0.60   526 1.01
ep[1]       0.48    0.02 1.10   -1.67   -0.24    0.43    1.17   2.74  4000 1.00
ep[2]      -1.15    0.02 1.47   -4.52   -1.99   -1.00   -0.18   1.39  4000 1.00
ep[3]       0.35    0.02 1.58   -2.63   -0.64    0.26    1.28   3.80  4000 1.00
ep[4]       1.65    0.03 1.07   -0.28    0.89    1.59    2.37   3.81  1066 1.00
ep[5]      -0.38    0.02 1.58   -3.74   -1.26   -0.32    0.60   2.57  4000 1.00
ep[6]       1.68    0.04 1.18   -0.36    0.83    1.60    2.46   4.14   947 1.00
ep[7]       0.67    0.02 1.51   -2.09   -0.31    0.54    1.52   3.99  4000 1.00
ep[8]       0.34    0.03 1.63   -2.73   -0.68    0.27    1.28   3.83  4000 1.00
ep[9]      -1.08    0.02 1.53   -4.57   -1.88   -0.91   -0.07   1.59  4000 1.00
ep[10]     -1.67    0.04 1.10   -3.96   -2.38   -1.63   -0.92   0.40   904 1.01
ep[11]     -2.05    0.04 1.28   -4.84   -2.84   -1.95   -1.15   0.15   819 1.01
ep[12]     -0.22    0.02 1.06   -2.32   -0.90   -0.20    0.45   1.91  4000 1.00
ep[13]     -0.33    0.02 1.01   -2.38   -0.99   -0.29    0.32   1.60  4000 1.00
ep[14]      1.01    0.04 1.34   -1.49    0.13    0.93    1.83   3.88  1353 1.00
ep[15]      1.41    0.04 1.12   -0.59    0.66    1.37    2.11   3.82  1024 1.00
ep[16]     -0.44    0.02 1.03   -2.53   -1.12   -0.41    0.26   1.49  4000 1.00
ep[17]     -2.01    0.05 1.28   -4.77   -2.85   -1.94   -1.09   0.27   764 1.01
ep[18]     -0.09    0.02 0.98   -2.05   -0.73   -0.09    0.54   1.85  4000 1.00
ep[19]     -0.55    0.02 1.24   -3.12   -1.33   -0.53    0.25   1.97  2753 1.00
ep[20]     -0.34    0.02 1.07   -2.49   -1.03   -0.30    0.35   1.76  4000 1.00
ep[21]      0.36    0.03 1.66   -2.66   -0.70    0.22    1.30   4.02  4000 1.00
ep[22]     -0.78    0.03 1.27   -3.45   -1.59   -0.71    0.07   1.56  1970 1.00
ep[23]      0.81    0.02 1.50   -1.81   -0.17    0.70    1.65   4.18  4000 1.00
ep[24]     -0.96    0.02 1.00   -3.07   -1.58   -0.93   -0.27   0.92  1807 1.00
ep[25]      2.23    0.05 1.46   -0.16    1.23    2.05    3.09   5.48   938 1.00
ep[26]     -2.17    0.05 1.21   -4.62   -2.97   -2.14   -1.33   0.02   627 1.01
ep[27]      0.96    0.03 1.08   -0.95    0.22    0.86    1.62   3.25  1746 1.00
ep[28]      0.93    0.02 1.00   -0.89    0.25    0.86    1.57   3.07  1740 1.00
ep[29]      1.32    0.03 1.08   -0.61    0.56    1.23    2.02   3.62  1372 1.00
ep[30]      0.32    0.02 0.95   -1.49   -0.29    0.31    0.92   2.28  4000 1.00
et[1]       1.26    0.35 5.47   -9.04    0.30    1.27    2.37  10.68   246 1.01
et[2]      -0.90    0.35 5.45  -11.51   -1.85   -0.73    0.22   8.55   245 1.01
et[3]      -0.72    0.35 5.45  -11.22   -1.67   -0.58    0.39   8.63   245 1.01
c[1]      -17.15    0.37 6.11  -29.21  -19.63  -16.95  -14.26  -6.68   272 1.01
c[2]      -14.99    0.37 6.02  -27.27  -17.28  -14.77  -12.30  -4.63   268 1.01
c[3]      -13.47    0.36 5.95  -25.73  -15.69  -13.23  -10.93  -3.02   267 1.01
c[4]      -10.90    0.36 5.85  -22.57  -12.94  -10.68   -8.55  -0.48   262 1.01
sigma[1]    1.70    0.04 0.59    0.65    1.30    1.66    2.06   2.97   240 1.01
sigma[2]    5.59    0.61 9.52    0.58    1.42    2.56    5.32  33.78   243 1.01
lp__     -108.76    0.59 8.43 -124.61 -114.26 -109.07 -103.96 -89.61   207 1.02

Samples were drawn using NUTS(diag_e) at Fri Jan 10 05:01:14 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

JAGSのときもそうだったが、drmの結果とはややちがう。推定方法の差によるものか。

データと期待値とをかさねてみる。

## plot
beta <- mean(extract(fit, pars = c("beta"))$beta)
c <- apply(extract(fit, pars = c("c"))$c, 2, mean)

x <- seq(5, 30, 0.1)
expected <- 5 * invlogit(beta * x - c[4]) +
            4 * (invlogit(beta * x - c[3]) -
                invlogit(beta * x - c[4])) +
            3 * (invlogit(beta * x - c[2]) -
                invlogit(beta * x - c[3])) +
            2 * (invlogit(beta * x - c[1]) -
                invlogit(beta * x - c[2])) +
            1 * (1 - invlogit(beta * x - c[1]))

col <- rainbow(n.plot)
plot(c1 ~ slope, data = data, pch = time, col = col[plot],
     las = 1, ylab = "cover", xlim = c(5,25))
lines(x, expected)

タグ：STAn