Fit a simple bayesian regression model using the iris dataset

iris_bayesian_lm.R

library(rjags)

data(iris)

independent <- iris$Petal.Length
dependent <- matrix(iris$Petal.Width)

lm.m <- lm(independent~dependent)

dependent <- cbind(1, dependent)
n <- nrow(dependent)
p <- ncol(dependent)

data <- list("y" = independent, "X" = dependent, "n" = n, "p" = p)
inits <- list(beta = coefficients(lm.m), sigma.sq = 4000)

modelstring <- "
  model{
  
   for(i in 1:n){
    y[i] ~ dnorm(mu[i], tau.sq)
   }
  
   mu <- X%*%beta
  
   for(i in 1:p){
    beta[i] ~ dnorm(0, 0.000001)
   }	
  
   tau.sq <- 1/sigma.sq
   sigma.sq ~ dgamma(4, 0.001)
  
   sigma.sq.prior.check ~ dgamma(4, 0.001)
  }
"

jags.m <- jags.model(file = textConnection(modelstring), data = data, inits = inits, n.chains = 3, n.adapt = 1000)

params <- c("beta", "sigma.sq", "sigma.sq.prior.check")
samps <- coda.samples(jags.m, params, n.iter = 40000)

plot(samps, density = FALSE)