JAGSで相関分析

独習のために、JAGSで相関分析したので、備忘録を残します。

参考にしたサイトは以下です。
他にもっと良いやり方や、コードにミスがあったら教えてくださいm(_ _)m。



まずは相関が0.75を想定したシミュレーションデータの生成です。
library(MASS)set.seed(1234) 
mu <- c(0,0) 
Sig <- matrix(c(1,0.75,0.75,1),2,2) 
x <- mvrnorm(100,mu,Sig) 
cor.test(x[,1],x[,2])
##
## Pearson's product-moment correlation
##
## data: x[, 1] and x[, 2] 
## t = 10.824, df = 98, p-value < 2.2e-16 
## alternative hypothesis: true correlation is not equal to 0 
## 95 percent confidence interval: 
## 0.6332983 0.8160617 
## sample estimates: 
## cor 
## 0.737925
通常の方法で相関を求めると0.73でした。




次にJAGSで相関をベイズ推定します。
library(runjags)# data 
data = list(x = x, N = nrow(x)) 
# JAGS code 
corrModel <- "model{
# likelihood
for(i in 1:N){
 x[i,1:2] ~ dmnorm(mu[], prec[ , ]) } 

# corresponding the precision matrix
 prec[1:2, 1:2] <- inverse(cov[,]) 

# constructing the covariance matrix
 cov[1,1] <- sigma[1]^2
 cov[1,2] <- sigma[1]*sigma[2]*rho
 cov[2,1] <- sigma[1]*sigma[2]*rho
 cov[2,2] <- sigma[2]^2 

# priors rho ~ dunif(-1,1) 
 mu[1] ~ dnorm(0, 1.0e-4) 
 mu[2] ~ dnorm(0, 1.0e-4)
 sigma[1] ~ dt(0, 1, 1)T(0,)
 sigma[2] ~ dt(0, 1, 1)T(0,) 
}"
fit <- run.jags(corrModel, data=data, monitor=c("mu", "sigma","rho"), method="rjparallel", n.chains=4)

print(fit)

## 
## JAGS model summary statistics from 40000 samples (chains = 4; adapt+burnin = 5000): 
## 
## Lower95 Median Upper95 Mean SD Mode MCerr 
## mu[1] -0.36156 -0.16063 0.037017 -0.16099 0.10178 -- 0.00090675 
## mu[2] -0.33062 -0.13158 0.064364 -0.13202 0.1003 -- 0.00089935 
## sigma[1] 0.88731 1.017 1.1685 1.0211 0.07206 -- 0.00071626 
## sigma[2] 0.86396 0.99901 1.1418 1.003 0.071438 -- 0.00070274 
## rho 0.63097 0.72913 0.81503 0.7257 0.04763 -- 0.00047156 
## 
## MC%ofSD SSeff AC.10 psrf 
## mu[1] 0.9 12599 -0.0094185 1.0002 
## mu[2] 0.9 12438 -0.0031534 1.0002 
## sigma[1] 1 10121 0.01468 1.0001 
## sigma[2] 1 10334 0.01241 1.0003 
## rho 1 10202 0.0071601 1.0002 
## 
## Total time taken: 2.5 minutes
Rhatを確認すると、きちんと収束しています。

相関は0.72と推定され、通常の方法の結果とほとんど同じになりました。

図も確認します。
plot(fit)
NewImage
NewImageNewImageNewImage
NewImage