SSブログ

CARBayes 4.0をつかってみるテスト [統計]

CARをつかった空間モデリングのできるRパッケージCARBayesがバージョン4.0になっていたので、つかってみるテストです。

CARBayes 4.0では以下のような関数で推定ができるようになっています。

S.independent
Independent and identically distributed random effects.
S.CARiar
Spatially correlated random effects modelled by the intrinsic autoregressive (IAR) model proposed by Besag et al (1991).
S.CARbym
Spatially correlated random effects modelled by the Besag-York-Mollie (BYM) model proposed by Besag et al (1991).
S.CARleroux
Spatially correlated random effects modelled by the proper CAR model proposed by Leroux (1999).
S.CARdissimilarity
Spatially correlated random effects modelled by the localised spatial smoothing approach proposed by Lee and Mitchell (2012).
S.CARcluster
Spatially correlated random effects modelled by the localised spatial smoothing approach proposed by Lee and Sarran (2014).

ということで、ためしてみました。使用したデータは、CARBayesで2次元のCARおよび 2次元でのStanによるCARモデルのテストで つかったものです。

## Data
## http://cse.ffpri.affrc.go.jp/hiroki/stat/Qglauca.csv
qg <- read.csv("Qglauca.csv")

X座標, Y座標, カウント数という形式でデータが格納されています。

> head(qg)
  X   Y N
1 0   0 0
2 0  -5 0
3 0 -10 2
4 0 -15 1
5 0 -20 0
6 0 -25 1

地図にしてみます。

## Plot
library(ggplot2)
p.map <- ggplot(qg) +
    geom_tile(aes(x = X + 2.5, y = Y - 2.5, fill = N)) +
    xlab("X") + ylab("Y") +
    coord_fixed(ratio = 1)
print(p.map)

Rplot001.png

ランダム効果の事前分布をIntrinsic CARとして、切片だけのモデルをあてはめてみます。

##
## CARBayes 4.0
##
library(CARBayes)
library(coda)

## Adjacent matrix
n.x <- length(levels(as.factor(qg$X)))
n.y <- length(levels(as.factor(qg$Y)))
n.sites <- n.x * n.y
W <- matrix(0, nrow = n.sites, ncol = n.sites)
for (x in 0:(n.x - 1)) {
  for (y in 1:n.y) {
    if (x > 0)       W[x * n.y + y, (x - 1) * n.y + y] <- 1
    if (x < n.x - 1) W[x * n.y + y, (x + 1) * n.y + y] <- 1
    if (y > 1)       W[x * n.y + y, x * n.y + y - 1] <- 1
    if (y < n.y)     W[x * n.y + y, x * n.y + y + 1] <- 1
  }
}

## Intrinsic Conditional Autoregressive model
fit.iar <- S.CARiar(N ~ 1, family = "poisson", data = qg, W = W,
                    burnin = 1000, n.sample = 41000, thin = 20)
geweke.diag(fit.iar$samples$beta)
print(fit.iar)

結果です。

> geweke.diag(fit.iar$samples$beta)

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

 var1 
-1.11 

> print(fit.iar)

#################
#### Model fitted
#################
Likelihood model - Poisson (log link function) 
Random effects model - Intrinsic CAR
Regression equation - N ~ 1

############
#### Results
############
Posterior quantiles and DIC

             Median    2.5%   97.5% n.sample % accept
(Intercept) -0.5193 -0.7963 -0.2915     2000     60.6
tau2         1.3035  0.7023  2.3333     2000    100.0

DIC =  488.7306       p.d =  53.76614 

データと、事後分布の中央値および95%信用区間をかさねてプロットしてみます。/p>

y.est <- apply(fit.iar$samples$fitted, 2, quantile,
               probs = c(0.025, 0.5, 0.975))
df.iar <- data.frame(x = 1:n.sites,
                     y = qg$N,
                     ymed = y.est[2, ],
                     ymax = y.est[1, ],
                     ymin = y.est[3, ])
p.iar <- ggplot(df.iar) +
    geom_point(aes(x = x, y = y)) +
    geom_pointrange(aes(x = x, y = ymed, ymax = ymax, ymin = ymin),
                    colour = "red", alpha = 0.5, shape = 3) +
    xlab("Index") + ylab("N") + ylim(0, 9)
print(p.iar)

Rplot.iar.png

事後分布の中央値を地図にしてみます。

p.map.iar <- ggplot(df.iar) +
    geom_tile(aes(x = qg[x, ]$X + 2.5,
                  y = qg[x, ]$Y - 2.5, fill = ymed)) +
    xlab("X") + ylab("Y") +
    coord_fixed(ratio = 1)
print(p.map.iar)

Rplot.iar.map.png

一方、Independent and identically distributed random effectsだと以下のようになります。

## Independent
fit.indep <- S.independent(N ~ 1, family = "poisson", data = qg,
                           burnin = 1000, n.sample = 41000, thin = 20)
geweke.diag(fit.indep$samples$beta)
print(fit.indep)

## Plot
y.est <- apply(fit.indep$samples$fitted, 2, quantile,
               probs = c(0.025, 0.5, 0.975))
p.indep <- ggplot(data.frame(x = 1:n.sites,
                             y = qg$N,
                             ymed = y.est[2, ],
                             ymax = y.est[1, ],
                             ymin = y.est[3, ])) +
    geom_point(aes(x = x, y = y)) +
    geom_pointrange(aes(x = x, y = ymed, ymax = ymax, ymin = ymin),
                    colour = "red", alpha = 0.5, shape = 3) +
    xlab("Index") + ylab("N") + ylim(0, 9)
print(p.indep)

結果です。

> print(fit.indep)

#################
#### Model fitted
#################
Likelihood model - Poisson (log link function) 
Random effects model - Independent
Regression equation - N ~ 1

############
#### Results
############
Posterior quantiles and DIC

             Median    2.5%   97.5% n.sample % accept
(Intercept) -0.3673 -0.5966 -0.1475     2000     59.7
sigma2       0.8024  0.4473  1.3087     2000    100.0

DIC =  529.9479       p.d =  74.13365 

Rplot.indep.png

また、Proper CARだと以下のようになりました。

## Proper conditional autoregressive prior
fit.leroux <- S.CARleroux(N ~ 1, family = "poisson", data = qg, W = W,
                          burnin = 1000, n.sample = 41000, thin = 20)
geweke.diag(fit.leroux$samples$beta)
print(fit.leroux)
plot(fit.leroux$samples$beta)

## Plot
y.est <- apply(fit.leroux$samples$fitted, 2, quantile,
               probs = c(0.025, 0.5, 0.975))
df.leroux <- data.frame(x = 1:n.sites,
                     y = qg$N,
                     ymed = y.est[2, ],
                     ymax = y.est[1, ],
                     ymin = y.est[3, ])
p.leroux <- ggplot(df.leroux) +
    geom_point(aes(x = x, y = y)) +
    geom_pointrange(aes(x = x, y = ymed, ymax = ymax, ymin = ymin),
                    colour = "red", alpha = 0.5, shape = 3) +
    xlab("Index") + ylab("N") + ylim(0, 9)
print(p.leroux)

結果

> print(fit.leroux)

#################
#### Model fitted
#################
Likelihood model - Poisson (log link function) 
Random effects model - Leroux CAR
Regression equation - N ~ 1

############
#### Results
############
Posterior quantiles and DIC

             Median    2.5%   97.5% n.sample % accept
(Intercept) -0.4837 -0.7389 -0.2637     2000     59.4
tau2         1.3634  0.7105  2.3508     2000    100.0
rho          0.9488  0.8010  0.9955     2000     59.8

DIC =  491.7704       p.d =  56.1957 

Rplot.leroux.png


タグ:CARBayes
nice!(1)  コメント(0)  トラックバック(0) 

nice! 1

コメント 0

コメントを書く

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

Facebook コメント

トラックバック 0