arm_ex_12.2_3.R

library(ggplot2)
# (c)  Investigate the change in partial pooling from (a)
# to (b) both graphically and numerically.

# Change in standard errors

# First and second model intercepts
df1 <- coef(mod1)$newpid[,1 , drop = F]
df2 <- coef(mod2)$newpid[,1 , drop = F]
names(df1) <- c("int")
names(df2) <- c("int")

# Confidence intervals for each intercept for both moels
df1$ci_bottom <- df1$int + (-2 * se.ranef(mod1)$newpid[,1])
df1$ci_upper <- df1$int + (2 * se.ranef(mod1)$newpid[,1])

df2$ci_bottom <- df2$int + (-2 * se.ranef(mod2)$newpid[,1])
df2$ci_upper <- df2$int + (2 * se.ranef(mod2)$newpid[,1])

# Now we need to compare whether the CI's shrunk from
# the first to the second model

# Calculate difference
df1$diff <- df1$ci_upper - df1$ci_bottom
df2$dff <- df2$ci_upper - df2$ci_bottom

# Create a df with both differences
df3 <- data.frame(cbind(df$diff, df2$dff))

# Create a difference out of that
df3$diff <- df3$X1 - df3$X2

# Graph it
ggplot(df3, aes(diff)) + geom_histogram(bins = 100) +
  xlim(0, 0.2)

# It looks like the difference is always higher than zero which
# means that in the second model the difference between
# the upper and lower CI is smaller than in the first model.
# This suggests we have greater certainty of our estimation
# by including the two predictors in the model.

# Numerically, the between-child variance in the first
# model was:
display(mod1)
11.65 / (11.65 + 7.31)

# 0.614%

# For the second model
display(mod2)
11.45 / (11.45 + 7.32)

# 0.6100%

# The between variance went down JUST a little, in line
# with the tiny reduction in the standard errors of the intercept.