Model II 回帰のMCMC

Model IIの回帰をMCMCでやってみる。

コード

library(rjags)

## data
set.seed(1234)

n <- 100	# sample size
sd.x <- 0.1
sd.y <- 0.1
b.x <- 0.5
x0 <- rnorm(n, 0, 1)
x <- x0 + rnorm(n, 0, sd.x)
y <- b.x * x0 + rnorm(n, 0, sd.y)

## model
model.file <- "model.bug"

m <- "## model II regression
var
  N, x[N], y[N], d[N], zero[N]
data {
  for (i in 1:N) {
    zero[i] <- 0
  }
}
model {
  for (i in 1:N) {
    zero[i] ~ dnorm(d[i], tau)
    d[i] <- (y[i] - beta.x * x[i] - beta) / sqrt(beta.x ^ 2 + 1)
  }
  beta ~ dnorm(0, 1.0E-4)
  beta.x ~ dnorm(0, 1.0E-4)
  tau ~ dgamma(1.0E-3, 1.0E-3)
  sigma <- 1 / sqrt(tau)
}
"
write(m, file = model.file)

## jags
model <- jags.model(model.file,
                    data = list(N = n,
                                x = x, y = y),
                    n.chain = 3, n.adapt = 2000)

## posterior
post <- coda.samples(model,
                     variable = c("beta", "beta.x",
                                  "sigma"),
                     n.iter = 30000, thin = 30)

summary(post)

## lmodel2
library(lmodel2)
lm2 <- lmodel2(y ~ x)
print(lm2)

MCMCの結果

> summary(post)

Iterations = 2030:32000
Thinning interval = 30 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean       SD  Naive SE Time-series SE
beta   0.01448 0.010955 0.0002000      0.0002219
beta.x 0.50971 0.010476 0.0001913      0.0001982
sigma  0.09438 0.006762 0.0001235      0.0001170

2. Quantiles for each variable:

           2.5%      25%     50%     75%   97.5%
beta   -0.00675 0.007192 0.01440 0.02166 0.03599
beta.x  0.48890 0.502642 0.50997 0.51683 0.53026
sigma   0.08192 0.089757 0.09382 0.09880 0.10888

> print(lm2)

Model II regression

Call: lmodel2(formula = y ~ x)

n = 100   r = 0.9794787   r-square = 0.9593786 
Parametric P-values:   2-tailed = 5.533767e-70    1-tailed = 2.766884e-70 
Angle between the two OLS regression lines = 0.967454 degrees

Regression results
  Method  Intercept     Slope  Angle (degrees)  P-perm (1-tailed)
1    OLS 0.01412038 0.5047313         26.78151                 NA
2     MA 0.01478723 0.5091002         26.98065                 NA
3    SMA 0.01573448 0.5153060         27.26232                 NA

Confidence intervals
  Method  2.5%-Intercept 97.5%-Intercept  2.5%-Slope 97.5%-Slope
1    OLS    -0.006982337      0.03522310   0.4839117   0.5255510
2     MA     0.011608416      0.01802050   0.4882742   0.5302828
3    SMA     0.012620794      0.01897651   0.4949068   0.5365461

Eigenvalues: 1.274834 0.008703556 

H statistic used for computing C.I. of MA: 0.0002781348 

結果をみるとヨサゲだが、これでよいのだろうか。

タグ：R MCMC