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))
par <- c("b0","b","r.squared","adjusted.r.squared")

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))