# Bayesian multiple regression issue in R

#### micdhack

##### New Member
I am attempting to generate a multiple regression model using JAGS on R and my results deviate substantially from MCMCregress and lm. I ran out of ideas checking what is it that I am missing. Any ideas why the results may be so different?

Code:
set.seed=06192014

dataset = data.frame(Y = c(58.6,58.7,52.6,50.7,55.1,69.9), M1 = c(38.3,33.1,36.8,37.8,35.7,36.4), M2 = c(0.48,0.52,0.493,0.498,0.497,0.502))
output = "Y"
predictor = c("M1","M2")

# Model parameters
n.simu <- 10000
n.burnin <- n.simu/2

D <- list(y=dataset[,output], x= dataset[,predictor], N=nrow(dataset), nPredictors = length(predictor))

packages = c("R.utils","R2jags","MCMCpack","coda","R2OpenBUGS","plyr","ggplot2","plyr","BEST")
for (package in packages){
if (require(package,character.only=TRUE)) {require(package,character.only=TRUE)} else {
install.packages(package)
library(package,character.only=TRUE) # This will ensure that the code fails if the installation failed
}
}

jags.bin <- function() {
for(i in 1:N) {
y[i] ~ dnorm(f[i], tau)
f[i] <- b0 + inprod(b[], x[i,])
}

# Priors
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
b0 ~ dnorm(0, 0.001)
for (j in 1:nPredictors) {
b[j] ~ dnorm(0,1.0E-12)
}

# R-squared
y.mean <- mean(y[])
for (i in 1:N) {
ss.res.temp[i] <- y[i] - f[i]
ss.res[i] <- pow(ss.res.temp[i], 2)
ss.reg.temp[i] <- f[i] - y.mean
ss.reg[i] <- pow(ss.reg.temp[i], 2)
ss.tot.temp[i] <- y[i] - y.mean
ss.tot[i] <- pow(ss.tot.temp[i], 2)
}
r.squared <- (sum(ss.reg[])) / (sum(ss.tot[]))
fraction <- (N-1)/(N-nPredictors-1)
adjusted.r.squared <- 1 - ((1 - r.squared) * (fraction))
}

write.model(jags.bin, "jags.txt")

m.jags <- jags.model("jags.txt", data = D, n.adapt = n.burnin, quiet = TRUE, n.chains = 4)
s <- coda.samples(m.jags, par, n.iter = n.simu - n.burnin, quiet = TRUE)

df = as.data.frame(as.matrix( s ))

for (i in (1:length(colnames(df)))) {
print(colnames(df)[i])
print(eval(parse(text = paste("mean(df[,",i,"])",sep=""))))
}

summary(lm(Y~.,data=dataset[,c(output,predictor)]))
summary(MCMCregress(Y~.,data=dataset[,c(output,predictor)], burnin = 1000, mcmc = 10000))