mixed_beer.R

#
# Generate random beer data
#
set.seed(314159)
ibv <- .5                    # Inner Beer Variation
mu.i <- rnorm(8,18,4)
sodium <- rnorm(6*8,mu.i,ibv)
beer <- rep(1:8,6)
plot(sodium~beer)

#
# Approach 1: Assume independence
#
t.test(sodium)

#
# Approach 2: Assume correlation of 1
#
b.mean <- tapply(sodium,beer,mean)
t.test(b.mean)

#
# Approach 3: Mixed Model Approach
#
install.packages("lme4")
library(lme4)

res.m <- lmer(sodium~1+(1|beer))
summary(res.m)

# Assuming IBV = .5 and seed = 314159
# 95% C.I. From Fixed Effects
18.051 + c(-1,1)*qt(.975,39)*1.227
#Correlation within each beer
12.013/(12.013+.238)

#
#  Try with IBV = 4
#

#
# 48 Unique Beers
#

set.seed(123456)
sodium <- rnorm(6*8,18,4)
beer <- rep(1:8,6)
plot(sodium~beer)
res.m <- lmer(sodium~1+(1|beer))
summary(res.m)
18.50 + c(-1,1)*qt(.975,39)*.5981
t.test(sodium)

#
# Is the random effect important
#
et.seed(314159)
ibv <- 1                    # Inner Beer Variation
mu.i <- rnorm(8,18,4)
sodium <- rnorm(6*8,mu.i,ibv)
beer <- rep(1:8,6)
plot(sodium~beer)
res.m <- lmer(sodium~1+(1|beer))
res <- lm(sodium~1)
summary(res)
summary(res.m)
AIC(res.m)
AIC(res)

# No correlation

set.seed(123456)
sodium <- rnorm(6*8,18,4)
beer <- rep(1:8,6)
plot(sodium~beer)
res.m <- lmer(sodium~1+(1|beer))
res <- lm(sodium~1)
AIC(res.m)
AIC(res)
summary(res.m)

#
# Another Example
#

install.packages("nlme")
library(nlme)
Rail
?Rail
attach(Rail)
plot(travel~Rail)

t.test(travel)
travel.m <- tapply(travel,Rail,mean)
travel.m
t.test(travel.m)

res.m <- lmer(travel~1+(1|Rail))
summary(res.m)
615.311/(615.311+16.167)

66.50 + c(-1,1)*qt(.975,11)*10.17