Webinar Part 4: Multilevel Modeling

webinar.june.2014.part4.R

#----------------------------------------------
# R code for Bayesian Computation Webinar
# Jim Albert - June 12, 2014
# albert@bgsu.edu
#----------------------------------------------

##################################################
# PART IV:  JAGS for Multilevel Modeling
# Require R packages Lahman, rjags, ggplot2
##################################################

library(Lahman)

# create a data frame with variables
# y = log(W/L)
# x = log(R/RA)
# yearID = season (between 1951 and 2010)

data <- subset(Teams, yearID>=1951 & yearID<=2010)
data <- data[, c("yearID", "W", "L", "R", "RA")]
data$y <- with(data, log(W/L))
data$x <- with(data, log(R/RA))

# some initial plots

years <- seq(1951, 2010, by=7)
data2 <- subset(data, yearID %in% years)

library(ggplot2)
ggplot(data2, aes(log(R/RA), log(W/L))) + 
  geom_point() +
  stat_smooth(method = "lm") +
  facet_wrap(~ yearID, ncol=3) +
  theme(axis.text = element_text(size = rel(2))) +
  theme(axis.title = element_text(size = rel(2))) +
  theme(strip.text = element_text(size = rel(2)))

# load in the rjags package

library(rjags)

# description of multilevel model

modelString = "
model {
for(i in 1:N){
mu.y[i] <- beta[j[i]] * x [i]
y[i] ~ dnorm(mu.y[i], tau[1])
}
for (p in 1:J){
beta[p] ~ dnorm(mu, tau[2])
}
mu ~ dnorm(0, .0001)
for(p in 1:2){
tau[p] <- pow(sigma[p], -2)
sigma[p] ~ dunif(0, 10)
}
}
"
# save model description to a file

writeLines(modelString, con="normalexchmodel.bug")

# defining data variables

N <- 1456
j <- data$yearID - 1950
y <- data$y
x <- data$x
J <- 60

# create object representing Bayesian model

jags <- jags.model('normalexchmodel.bug',
                   data = list('y' = y, "j" = j, 'x' = x,
                               "N" = N, "J" = J),
                   n.chains = 1,
                   n.adapt = 1000)

# perform 5000 iterations of MCMC

update(jags, 5000)

# take 10,000 iterations, storing simulated draws of beta
# and sigma

posterior <- coda.samples(jags, c("beta", "sigma"),
                          n.iter=10000)

# find summaries of parameters

S <- summary(posterior)

# extract posterior means of beta's

Multilevel <- data.frame(yearID = 1951:2010,
                         Slope = S$statistics[1:60, "Mean"],
                         Type = "Multilevel")

# find individual slopes

library(dplyr)

Individual <- summarize(group_by(data, yearID), 
              Slope = coef(lm(log(W / L) ~ 0 + log(R / RA))))
Individual$Type = "Individual"

# display the individual and multilevel estimates
# use ggplot2

library(ggplot2)
print(ggplot(rbind(Individual, Multilevel), 
       aes(yearID, Slope, color=Type)) +
  geom_point(size=4) + 
  geom_line(size=1.5) +
  labs(title="Estimating Pythagorean Slopes") +
  theme(plot.title = element_text(size = rel(2))) +
  theme(axis.text = element_text(size = rel(2))) +
  theme(axis.title = element_text(size = rel(2))) +
  theme(legend.text = element_text(size = rel(1.5))) +
  theme(legend.title = element_text(size = rel(1.5))))

##################################################################