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)
ランダム効果の事前分布を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)
事後分布の中央値を地図にしてみます。
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)
一方、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
また、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
タグ:CARBayes
コメント 0