# ============================================================================= # HSCI 410 Public Health Assessment and Analysis - Lesson 4: Generalized Linear Models # Answer key for the in-lesson R activities # Data file(s): phaa_survey_clean.csv, phaa_followup.csv (download from the lesson page; save in your R working directory) # Packages: nnet, MASS, brant, AER (install once with install.packages(c("nnet", "MASS", "brant", "AER"))) # Reproduces every code block in the lesson, then answers each activity question. # ============================================================================= # Part 1 (Ordinal and Multinomial Models) uses phaa_survey_clean.csv, the file # saved at the end of Lesson 2. Part 2 (Count and Rate Data) uses # phaa_followup.csv, which adds years of follow-up and GP-visit counts to the # same cohort. Run the script from the folder that holds both files. # ============================================================================= # ==== Part 1: Ordinal and Multinomial Models ==== # ============================================================================= # ==== Section 3: Proportional-Odds Model / Activity: multinomial and proportional-odds models in R ==== library(nnet); library(MASS); library(brant) phaa <- read.csv("phaa_survey_clean.csv", stringsAsFactors = FALSE, na.strings = c("", "NA")) phaa$region <- relevel(factor(phaa$region), ref = "Vancouver") phaa$education <- factor(phaa$education, levels = c("Less than high school", "High school", "Some college", "Bachelor's", "Graduate degree"), ordered = TRUE) # The outcomes we are about to model table(phaa$region) # Vancouver, Burnaby, Other, Richmond, Surrey table(phaa$education) # 5 ordered levels # 1. MULTINOMIAL: nominal outcome (region) mn <- multinom(region ~ age + gender + smoker, data = phaa, trace = FALSE) exp(coef(mn)) # relative-risk ratios z <- summary(mn)$coefficients / summary(mn)$standard.errors round((1 - pnorm(abs(z))) * 2, 3) # p-values # 2. PROPORTIONAL-ODDS: ordered outcome (education) po <- polr(education ~ age + gender + smoker, data = phaa, Hess = TRUE) summary(po) exp(cbind(OR = coef(po), confint(po))) # OR + 95% CI # 3. Brant test: is the proportional-odds assumption defensible? # (brant() warns that a few education x gender cells are empty: only 27 # respondents are Non-binary. The warning is expected; it is not an error.) brant(po) # 4. Multinomial fit of the SAME ordered outcome, the fall-back to compare mn_edu <- multinom(education ~ age + gender + smoker, data = phaa, trace = FALSE) AIC(po, mn_edu) # lower AIC = better-fitting model # ---- Activity questions (Part 1, box rActivity-410-6-1) ---------------------- rrr <- exp(coef(mn)) pmat <- round((1 - pnorm(abs(z))) * 2, 3) cat("\nQ1: RRR for smokerYes, Burnaby vs Vancouver =", round(rrr["Burnaby", "smokerYes"], 2), " (p =", pmat["Burnaby", "smokerYes"], ")\n") cat(" Smokers are about", round(100 * (1 - rrr["Burnaby", "smokerYes"])), "% less likely than non-smokers to live in Burnaby rather than Vancouver,\n", " holding age and gender constant, but the p-value shows this is not\n", " statistically significant. Region was simulated independently of the\n", " predictors, so only age for Burnaby (RRR 0.98, p = 0.038) reaches p < .05,\n", " which is what chance alone produces across 16 tests.\n") # The page's model answer used to cite a 'Fraser Valley' category; region has no # such level (Vancouver, Burnaby, Surrey, Richmond, Other). or_po <- exp(cbind(OR = coef(po), confint(po))) cat("\nQ2: OR for age (proportional odds) =", round(or_po["age", "OR"], 3), " 95% CI", round(or_po["age", "2.5 %"], 3), "to", round(or_po["age", "97.5 %"], 3), "\n") cat(" Each extra year of age multiplies the odds of being at-or-above any given\n", " education level by 0.99 (about 0.6% lower odds per year); the CI includes 1,\n", " so age is not associated with education in this sample. Under proportional\n", " odds the same OR applies at every cut-point.\n") br <- brant(po) cat("\nQ3: Brant omnibus X2 =", round(br["Omnibus", "X2"], 2), "on", br["Omnibus", "df"], "df, p =", round(br["Omnibus", "probability"], 2), "\n") cat(" Predictors with p < .05:", ifelse(any(br[-1, "probability"] < 0.05), paste(rownames(br)[-1][br[-1, "probability"] < 0.05], collapse = ", "), "none (smallest is age, p = 0.14)"), "\n") aic_tab <- AIC(po, mn_edu) cat(" AIC: po =", round(aic_tab["po", "AIC"], 1), " mn_edu =", round(aic_tab["mn_edu", "AIC"], 1), " -> keep the proportional-odds model (lower AIC with 12 fewer parameters).\n") # ============================================================================= # ==== Part 2: Count and Rate Data ==== # ============================================================================= # ==== Section 6: Poisson Regression Model & Interpretation / Activity: Poisson with an offset and a negative-binomial fall-back ==== library(MASS); library(AER) phaa <- read.csv("phaa_followup.csv", stringsAsFactors = FALSE, na.strings = c("", "NA")) phaa$smoker <- factor(phaa$smoker, levels = c("No","Yes")) # NOTE (answer key): one participant (P0533) has fu_years = 0, so log(fu_years) # is -Inf and glm() stops with "NA/NaN/Inf in 'y'". A person with no follow-up # time contributes no person-time, so the row is dropped before modelling. phaa <- phaa[phaa$fu_years > 0, ] # one participant has 0 years of follow-up; log(0) is undefined # 1. A peek at the count outcome summary(phaa$gp_visits) hist(phaa$gp_visits, breaks = 30, main = "GP visits during follow-up", xlab = "Visits") table(phaa$gp_visits == 0) # 34 participants with no visits # 2. Poisson regression with offset to model the RATE per person-year fit_rate <- glm(gp_visits ~ age + smoker + hypertension + offset(log(fu_years)), family = poisson, data = phaa) summary(fit_rate) exp(coef(fit_rate)) # incidence-rate ratios exp(confint(fit_rate)) # 3. Goodness of fit: Pearson chi^2 / df ~ 1 = good sum(residuals(fit_rate, type = "pearson")^2) / fit_rate$df.residual # 4. Formal overdispersion test dispersiontest(fit_rate) # 5. Negative binomial fall-back when Poisson is overdispersed fit_nb <- glm.nb(gp_visits ~ age + smoker + hypertension + offset(log(fu_years)), data = phaa) AIC(fit_rate, fit_nb) cbind(Poisson = exp(coef(fit_rate)), NegBin = exp(coef(fit_nb))) # ---- Activity questions (Part 2, box rActivity-410-7-1) ---------------------- irr <- exp(coef(fit_rate)) irr_ci <- exp(confint(fit_rate)) cat("\nQ1: IRR for smokerYes =", round(irr["smokerYes"], 2), " 95% CI", round(irr_ci["smokerYes", 1], 2), "to", round(irr_ci["smokerYes", 2], 2), "\n") cat(" Smokers have about", round(100 * (irr["smokerYes"] - 1)), "% more GP visits per person-year than non-smokers,\n", " holding age and hypertension constant; the CI excludes 1.\n") pearson_ratio <- sum(residuals(fit_rate, type = "pearson")^2) / fit_rate$df.residual dt <- dispersiontest(fit_rate) cat("\nQ2: Pearson chi-square / df =", round(pearson_ratio, 2), " (a little above 1, so mild overdispersion)\n") cat(" dispersiontest(): estimated dispersion =", round(dt$estimate, 2), ", z =", round(dt$statistic, 2), ", p =", signif(dt$p.value, 2), "\n") cat(" Both point the same way: the variance exceeds the mean by roughly 20-25%,\n", " which is modest but statistically clear, so Poisson SEs are slightly too small.\n") aic2 <- AIC(fit_rate, fit_nb) se_tab <- cbind(Poisson_SE = sqrt(diag(vcov(fit_rate))), NegBin_SE = sqrt(diag(vcov(fit_nb)))) cat("\nQ3: AIC Poisson =", round(aic2["fit_rate", "AIC"], 1), " NegBin =", round(aic2["fit_nb", "AIC"], 1), " (NegBin lower by", round(aic2["fit_rate", "AIC"] - aic2["fit_nb", "AIC"], 1), "points)\n") print(round(cbind(Poisson = exp(coef(fit_rate)), NegBin = exp(coef(fit_nb))), 3)) print(round(se_tab, 4)) cat(" Prefer the negative binomial. The IRRs are almost identical (smokerYes 1.22 vs 1.23);\n", " what changes is the standard errors, about 13-20% larger under NegBin (theta ~ 56),\n", " so the CIs widen a little and the p-values become slightly less extreme.\n") # ---- Stretch (not required): a linear model is the wrong tool for counts ---- lm_naive <- lm(gp_visits ~ age + smoker, data = phaa) range(predict(lm_naive)) # non-integer predictions with no variance-mean link