This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more of my R tutorials visit http://mikewk.com/statistics.

Contents

Just Another Gibbs Sampler

Load the rjags package.

library(rjags)

Read in the data - thesis data with generic variable names

data <- read.csv('generic.csv', header=TRUE)
head(data)
##   X x1 x2 x3 z1 z2 z3 w1 w2 w3 y1  y2  y3 y4 y5 y6 y7
## 1 1  2  2  1  2  3  4  2  4  3  2 0.0 0.0  0  0  1  0
## 2 2  1  1  1  5  5  5  2  2  2  0 0.0 0.0  0  0  0  0
## 3 3  1  1  1  3  4  5  4  5  1  2 1.0 0.0  2  0  0  0
## 4 4  2  1  1  3  4  4  2  2  2  2 0.5 1.5  0  0  0  0
## 5 5  4  3  3  1  1  1  5  5  1  3 3.0 2.5  3  3  3  1
## 6 6  2  1  1  3  3  3  1  3  3  1 1.5 1.0  1  0  0  0

Specify model variables

N <- nrow(data)
epsilon <- rnorm(N, 0, 1)

data$yVar <- with(data, (y1+y2+y3+y4+y5)/5)
data$xVar <- with(data, (x1+x2+x3)/3)
data$wVar <- with(data, (w1+w2+w3)/3)
data$zVar <- with(data, (z1+z2+z3)/3)

yVar <- with(data, (y1+y2+y3+y4+y5)/5)
xVar <- with(data, (x1+x2+x3)/3)
wVar <- with(data, (w1+w2+w3)/3)
zVar <- with(data, (z1+z2+z3)/3)

Example 1: Simple Linear Regression

Specify model and save as .bug file.

rjags.ex.1.bug <-'
model {
   for (i in 1:N){
     yVar[i] ~ dnorm(yVar.hat[i], tau)
     yVar.hat[i] <- a + b * xVar[i]
   }
   a ~ dnorm(0, .0001)
   b ~ dnorm(0, .0001)
   tau <- pow(sigma, -2)
   sigma ~ dunif(0, 100)
}
'

Run model

jags1 <- jags.model('rjags.ex.1.bug',
            data = list('x' = xVar,
                'y' = yVar,
                'N' = N),
            n.chains = 4,
            n.adapt = 100)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 944
## 
## Initializing model
# update(jags1, 1000)

Generate posterior samples

coefs1 <- jags.samples(jags1,
  c('a', 'b'),
  2000)

Extract estimates - means from posterior samples and results from frequentist lm() model

intercept <- mean(coefs1$a[])
b <- mean(coefs1$b[])

lm1 <- lm(yVar ~ xVar, data=data)

Compare the results

# results from frequentist lm() model
coefficients(lm1)
## (Intercept)        xVar 
##   0.5242269   0.1009760
# results from bayesian jags.samples() model
data.frame(intercept, b)
##   intercept        b
## 1  0.523091 0.101616

Example 2: Linear Regression with Multiple Predictors

Specify model with multiple predictors as .bug file

rjags.ex.2.bug <- '
model {
  for (i in 1:N){
    y[i] ~ dnorm(y.hat[i], tau)
    y.hat[i] <- a + 
      b1 * x[i] + 
      b2 * w[i] +
      b3 * z[i]
  }
  a ~ dnorm(0, .0001)
  b1 ~ dnorm(0, .0001)
  b2 ~ dnorm(0, .0001)
  b3 ~ dnorm(0, .0001)
  
  tau <- pow(sigma, -2)
  sigma ~ dunif(0, 100)
}
'

Run model

jags2 <- jags.model('rjags.ex.2.bug',
                   data = list('x' = xVar,
                               'w' = wVar,
                               'z' = zVar,
                               'y' = yVar,
                               'N' = N),
                   n.chains = 4,
                   n.adapt = 100)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 2192
## 
## Initializing model
#update(jags2, 1000)

Generate posterior samples and compare means with estimates from lm() model.

coefs2 <- jags.samples(jags2,
                      c('a', 'b1', 'b2', 'b3'),
                      1000)

intercept <- mean(coefs2$a[])
b.xVar <- mean(coefs2$b1[])
b.wVar <- mean(coefs2$b2[])
b.zVar <- mean(coefs2$b3[])

# results from lm()
lm2 <- lm(yVar ~ xVar + wVar + zVar, data=data)

# means from jags.sample()
coefficients(lm2)
## (Intercept)        xVar        wVar        zVar 
##   0.5400022  -0.0491332   0.3226198  -0.2022058

Example 3: Linear Regression with Interaction Term

Specify model as .bug file

rjags.ex.3.bug <-'
model {
    for (i in 1:N){
y[i] ~ dnorm(y.hat[i], tau)
y.hat[i] <- a + 
b1 * x[i] + 
b2 * w[i] +
b3 * z[i] +
b4 * x[i] * w[i]
}
a ~ dnorm(0, .0001)
b1 ~ dnorm(0, .0001)
b2 ~ dnorm(0, .0001)
b3 ~ dnorm(0, .0001)
b4 ~ dnorm(0, .0001)

tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
'

Run model

jags3 <- jags.model('rjags.ex.3.bug',
                    data = list('x' = xVar,
                                'w' = wVar,
                                'z' = zVar,
                                'y' = yVar,
                                'N' = N),
                    n.chains = 4,
                    n.adapt = 100)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 2295
## 
## Initializing model
#update(jags3, 1000)

Extract and compare estimates

coefs3 <- jags.samples(jags3,
                       c('a', 'b1', 'b2', 'b3', 'b4'),
                       3000)

intercept <- mean(coefs3$a[])
b.xVar <- mean(coefs3$b1[])
b.wVar <- mean(coefs3$b2[])
b.zVar <- mean(coefs3$b3[])
b.xwIx <- mean(coefs3$b4[])

# frequentist estimates
lm3 <- lm(yVar ~ xVar + wVar + zVar + xVar*wVar, data=data)
round(coefficients(lm3), 3)
## (Intercept)        xVar        wVar        zVar   xVar:wVar 
##       0.754      -0.164       0.258      -0.203       0.035
# bayesian estimates
round(data.frame(intercept, b.xVar, b.wVar, b.zVar, b.xwIx),3)
##   intercept b.xVar b.wVar b.zVar b.xwIx
## 1     0.723 -0.161   0.26 -0.198  0.035

Plots

Some quick and easy plots

b1.mcmc <- as.mcmc.list(coefs3$b1)
b4.mcmc <- as.mcmc.list(coefs3$b4)
plot(b1.mcmc, smooth=TRUE)

plot(b4.mcmc, smooth=TRUE)

geweke.plot(b1.mcmc)

traceplot(b1.mcmc)

acfplot(b4.mcmc)

autocorr.plot(b4.mcmc)

crosscorr.plot(b4.mcmc)