Stan: 多変量正規分布の空間自己相関モデリング：Taglibro de H：SSブログ

	ブログをはじめるログイン

Stan: 多変量正規分布の空間自己相関モデリング　[統計] [編集]

Hostetler and Chandler (in press)をよんでいたら、調査サイト間の空間自己相関を多変量正規分布でモデリングするというオプションが提示されていたので、やってみたテストです。

使用したデータは久保みどりぼん11章のものです。

ちなみに、StanによるCARモデリングはberoberoさんが実装されています。ついでに、RのhSDMパッケージによる例は自分で2年ばかり前に、CARBayesによる例は去年やっています。

ということで、コードです。

load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/Y.RData"))
N_site <- length(Y)

library(rstan)

stan_code <- "
data {
  int<lower=0> N_site;
  int<lower=0> Y[N_site];
}
transformed data {
  vector[N_site] mu;

  for (i in 1:N_site) {
    mu[i] <- 0;
  }
}
parameters {
  vector[N_site] r;
  real beta;
  real s[2];
}
transformed parameters {
  cov_matrix[N_site] Sigma;

  for (i in 1:N_site) {
    for (j in 1:N_site) {
      if (i == j) {
        Sigma[i, j] <- s[1];
      } else if (abs(i - j) == 1) {
        Sigma[i, j] <- s[2];
      } else {
        Sigma[i, j] <- 0;
      }
    }
  }
}
model {
  for (j in 1:N_site) {
    Y[j] ~ poisson(exp(beta + r[j]));
  }
  r ~ multi_normal(mu, Sigma);
  beta ~ normal(0, 100);
  s[1] ~ uniform(0, 100);
  s[2] ~ normal(0, 100);
}
generated quantities {
  vector[N_site] Y_hat;

  Y_hat <- exp(beta + r);
}
"

fit <- stan(model_code = stan_code,
            data = list(Y = Y, N_site = N_site),
            pars = c("beta", "r", "s", "Y_hat"),
            iter = 2000, warmup = 1000, thin = 1)

print(fit)

pars <- extract(fit)

## plot
plot(1:N_site, Y, type = "p", xlab = "Site", ylim = c(0, 25), las = 1)
lines(1:N_site, apply(pars$Y_hat, 2, mean))
lines(1:N_site, apply(pars$Y_hat, 2, quantile, probs = 0.9), lty = 2)
lines(1:N_site, apply(pars$Y_hat, 2, quantile, probs = 0.1), lty = 2)

サイト間のランダム効果rの事前分布が、平均mu、分散共分散行列Sigmaの多変量正規分布に従うとして、隣のサイト間の共分散がs[2]、その他のサイト間の共分散が0としてします。実行すると、警告が山のようにでますが、うまく収束します。

結果