library(ggplot2) library(dplyr) library(lme4) # might need to install # baseball df <- read.csv("http://roualdes.sfo3.digitaloceanspaces.com/data/bball1970.csv") # no pooling fit_np <- glm(cbind(Hits, AB - Hits) ~ Player, data = df, family = "binomial") phat_np <- plogis(predict(fit_np)) # partial pooling fit_pp <- glmer(cbind(Hits, AB - Hits) ~ (1 | Player), data = df, family = "binomial") phat_pp <- plogis(predict(fit_pp)) # complete pooling phat_cp <- rep(plogis(fixef(fit_pp)), length(df$Player)) ## will be the same as ## fit_cp <- glm(cbind(Hits, AB - Hits) ~ 1, data = df, family = "binomial") ## phat_cp <- plogis(predict(fit_cp)) ## since there are the same number of observations in each level of Player pdf <- data.frame(phat = c(phat_cp, phat_np, phat_pp), player = rep(df$Player, 3), pooling = rep(c("complete", "no", "partial"), each = length(df$Player))) ggplot(pdf, aes(player, phat, color = pooling, shape = pooling)) + geom_point(size = 1.8) + labs(x = "Player", y = "Predicted batting average") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) # radon: random intercepts df <- read.csv("http://roualdes.sfo3.digitaloceanspaces.com/data/radon.csv") # no pooling fit_np <- lm(log_radon ~ county, data = df) lr_np <- predict(fit_np) # partial pooling fit_pp <- lmer(log_radon ~ (1 | county), data = df) lr_pp <- predict(fit_pp) # complete pooling ## but when there are not the same number of observations in each level, ## fit_cp <- lm(log_radon ~ 1, data = df) ## lr_cp <- predict(fit_cp) ## does not produce the correct global mean ## use this instead lr_cp <- rep(fixef(fit_pp), length(df$county)) pdf <- data.frame(yhat = c(lr_cp, lr_np, lr_pp), county = rep(df$county, 3), pooling = rep(c("complete", "no", "partial"), each = length(df$county))) ggplot(pdf, aes(yhat, county, color = pooling, shape = pooling)) + geom_point(size = 1.8) + labs(y = "County", x = "Predicted log(radon)") + theme_minimal() # carnivora: random slopes df <- read.csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/carnivora.csv") %>% mutate(brain_weight = log10(SB), body_weight = log10(SW)) %>% select(Family, body_weight, brain_weight) %>% na.omit # no pooling fit_np <- lm(brain_weight ~ Family + Family : body_weight, data = df) # lm(brain_weight ~ Family * body_weight) families <- sort(unique(df$Family)) K <- length(families) beta_np <- coef(fit_np)[(K+1):(2*K)] # partial pooling fit_pp <- lmer(brain_weight ~ Family + (body_weight | Family), data = df, control = lmerControl(optimizer = "bobyqa")) beta_pp <- coef(fit_pp)$Family$body_weight ## complete pooling fit_cp <- lm(brain_weight ~ body_weight, data = df) beta_cp <- rep(coef(fit_np)[2], K) pdf <- data.frame(beta = c(beta_cp, beta_np, beta_pp), Family = rep(families, 3), pooling = rep(c("complete", "no", "partial"), each = length(families))) ggplot(pdf, aes(beta, Family, color = pooling, shape = pooling)) + geom_point(size = 1.8) + labs(y = "Family", x = "Estimated coefficient") + theme_minimal() # references on hierarchical/mixed effects models # OER: Open Educational Resources are both free and publicly accessible # OER 1. https://www.tjmahr.com/plotting-partial-pooling-in-mixed-effects-models/ # OER 2. https://m-clark.github.io/mixed-models-with-R/introduction.html # $93 3. https://www.routledge.com/Statistical-Rethinking-A-Bayesian-Course-with-Examples-in-R-and-STAN/McElreath/p/book/9780367139919