My First Bayesian (Markov Chain Monte Carlo) Simulation

# My First Bayesian (Markov Chain Monte Carlo) Simulation

# I know very little about Baysian methods and this post will probably not reveal much information information.  However, it should hopefully lead to many more fruitful posts in the future.

# You will need the package MCMCglmm (Markov Chain Monte Carlo)
require(MCMCglmm)

# Let's set the random seed
set.seed(121)

nobs = 100

x=rnorm(nobs)
e=rnorm(nobs)
y=10+2*x+e*20

d <- data.frame(x=x,e=e,y=y)

mc <- MCMCglmm(y~x, data=d)

plot(mc$Sol)

# This will give the estimates of the coefficients that we are typically looking for.
posterior.mode(mc$Sol)

# Probability that the intercept is greater than 0
table(mc$Sol[,1] > 0)/1000
# 100% of the estimates of the intercept indicate that it is greater than zero.

# Probability that the coefficient on x is greater than 0
table(mc$Sol[,2] > 0)/1000
# 888% of estimate results indicate that the coefficient on x is greater than 0.

# The 95% confidence interval of the MCMC estimates can be found:
HPDinterval(mc$Sol, 0.95)

summary(lm(y~x,data=d))
# We can see that the standard linear model fails to reject the null for the coefficient on x even though it is large.

# To read a small tutorial on the use of MCMC written by the package's author J.D. Hatfield
# http://cran.uvigo.es/web/packages/MCMCglmm/vignettes/Tutorial.pdf