library(haven) # import Stata data and labels
library(dplyr) # data manipulation
library(tidyr) # reshape data: pivot_longer()
library(ggplot2) # graphs
library(MASS) # negative binomial: glm.nb()
library(pscl) # zero-inflated models: zeroinfl()
library(cleanplots) # graph scheme
cleanplots_defaults()
source("https://raw.githubusercontent.com/tdmize/Rfunctions/main/Count_helper_functions.R")
gss <- read_dta("https://tdmize.github.io/data/data/cda_gss.dta")
gss <- gss |> filter(year >= 2002 & year <= 2014)
# Data management: Christian (Protestant, Catholic, Christian) vs not religious
gss <- gss |>
mutate(notrelig = case_when(relig %in% c(1, 2, 11) ~ 0, relig == 4 ~ 1))
# Keep model variables; drop any missing data
gss <- gss[c("mntlhlth", "physhlth", "woman", "race3", "year", "college", "parent",
"faminc", "metro", "notrelig", "married", "age")]
gss <- gss[complete.cases(gss), ]
# Treat binary/nominal variables as factors
gss <- gss |>
mutate(across(c(woman, married, race3, college, parent), as_factor),
notrelig = factor(notrelig, labels = c("Christian", "Not Religious")),
year = factor(year))
# 6.12.c - compare predictions to choose model
count_f <- mntlhlth ~ woman + married + poly(age, 2) + poly(faminc, 2) + race3 + college +
parent + notrelig + year
zero_f <- mntlhlth ~ woman + married + poly(age, 2) + poly(faminc, 2) + race3 + college +
parent + notrelig + year |
woman + married + poly(age, 2) + faminc + race3 + college + parent + notrelig + year
pmod <- glm(count_f, family = "poisson", data = gss)
nbmod <- glm.nb(count_f, data = gss)
zipmod <- zeroinfl(zero_f, data = gss, dist = "poisson")
zinbmod <- zeroinfl(zero_f, data = gss, dist = "negbin")
# Observed proportions and average predicted probabilities of 0 to 9 days
meanpred <- data.frame(k = 0:9) |>
mutate(Observed = sapply(k, function(j) mean(gss$mntlhlth == j)),
Poisson = sapply(k, function(j) mean(dpois(j, fitted(pmod)))),
"Neg. Binomial" = sapply(k, function(j)
mean(dnbinom(j, size = nbmod$theta, mu = fitted(nbmod)))),
"Z.I. Poisson" = colMeans(predict(zipmod, type = "prob", at = 0:9)),
"Z.I.N.B." = colMeans(predict(zinbmod, type = "prob", at = 0:9))) |>
pivot_longer(-k, names_to = "model", values_to = "pr") |>
mutate(model = factor(model, levels = unique(model)))
ggplot(meanpred, aes(x = k, y = pr, color = model)) +
geom_line() + geom_point() +
scale_x_continuous(breaks = 0:9) +
scale_y_continuous(breaks = seq(0, .6, .1), limits = c(0, .6)) +
labs(title = "[6.12.c] Comparison of Mean Predictions from Various Count Models",
x = "Number of Bad Mental Health Days", y = "Pr(y = k)", color = NULL)
![[6.12.c] Observed proportions of 0 to 9 bad mental health days with the average predicted probabilities from four count models: about 60% observed at 0, which the negative binomial, zero-inflated Poisson, and zero-inflated negative binomial all match closely while the Poisson predicts only about 6%; the Poisson also over-predicts 1 to 6 days, and the zero-inflated Poisson under-predicts 1 to 3 days.](fig/cda-06_12_c-mntlhlth.png)
![[6.12.c] Observed proportions of 0 to 9 bad mental health days with the average predicted probabilities from four count models: about 60% observed at 0, which the negative binomial, zero-inflated Poisson, and zero-inflated negative binomial all match closely while the Poisson predicts only about 6%; the Poisson also over-predicts 1 to 6 days, and the zero-inflated Poisson under-predicts 1 to 3 days.](6-12-model-selection_files/figure-html/unnamed-chunk-1-1.png)
![[6.12.a] Observed minus average predicted probability of 0 to 9 bad mental health days for the Poisson, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial models: the Poisson misses badly, under-predicting zeros by about 0.54 and over-predicting 1 to 4 days by about 0.1, while the other three models stay close to 0.](fig/cda-06_12_a-countfit.png)
![[6.12.a] Observed minus average predicted probability of 0 to 9 bad mental health days for the Poisson, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial models: the Poisson misses badly, under-predicting zeros by about 0.54 and over-predicting 1 to 4 days by about 0.1, while the other three models stay close to 0.](6-12-model-selection_files/figure-html/unnamed-chunk-2-1.png)
![[6.12.b] Observed minus average predicted probability of 0 to 9 bad mental health days without the Poisson model: the negative binomial over-predicts 1 day by about 0.06, the zero-inflated Poisson under-predicts 1 to 3 days by up to about 0.06 and over-predicts 6 to 9 days, and the zero-inflated negative binomial stays within about 0.025 of 0.](fig/cda-06_12_b-countfit.png)
![[6.12.b] Observed minus average predicted probability of 0 to 9 bad mental health days without the Poisson model: the negative binomial over-predicts 1 day by about 0.06, the zero-inflated Poisson under-predicts 1 to 3 days by up to about 0.06 and over-predicts 6 to 9 days, and the zero-inflated negative binomial stays within about 0.025 of 0.](6-12-model-selection_files/figure-html/unnamed-chunk-3-1.png)