# =============================================================================
# HSCI 410 Lesson 4: Generalized Linear Models -- ANSWER KEY
# Author: Kiffer G. Card, PhD - Faculty of Health Sciences, SFU
# -----------------------------------------------------------------------------
# Page:  modules/HSCI_410_Lesson_4_Generalized_Linear_Models.html
# Data:  r-activities/data/phaa_outcomes.csv  (sections 1 to 4)
#        r-activities/data/phaa_followup.csv  (offset example in section 4)
#
# This script reproduces every code block on the page, in page order, and then
# answers each numbered activity question with a printed result. Run it from a
# folder that holds both CSV files (for example r-activities/data/).
#
# Packages: MASS and nnet are installed with R; brant may need
#   install.packages("brant")
#
# phaa_outcomes.csv is simulated teaching data (see
# r-activities/generators/00_generate_phaa_outcomes.R). The previous version of
# this lesson and its answer key are kept as
# HSCI_410_Lesson_4_Generalized_Linear_Models_old_version.{html,R}.
# =============================================================================


# =============================================================================
# SECTION 1, ACTIVITY 1: Linear regression and its checks (Activity 4.1 on the lesson page)
# =============================================================================
# install.packages(c("MASS", "nnet", "brant"))   # run once, if not yet installed
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
nrow(phaa)                   # how many people are in the analysis?
summary(phaa$systolic_bp)    # a first look at the outcome

lin <- lm(systolic_bp ~ age + smoker, data = phaa)   # fit the linear model
summary(lin)                 # coefficients, standard errors and p-values
confint(lin)                 # 95% confidence intervals

par(mfrow = c(2, 2))         # show four plots on one screen
plot(lin)                    # the four standard checks of a linear model
par(mfrow = c(1, 1))         # back to one plot per screen

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Smoking coefficient (mmHg) and 95% CI:\n")
print(round(c(estimate = coef(lin)[["smokerYes"]], confint(lin)["smokerYes", ]), 2))
# Smokers average about 6.06 mmHg higher than non-smokers of the same age
# (95% CI 4.33 to 7.79).
# Q2. Residuals vs Fitted checks linearity (flat red line) and equal variance
#     (even band); Q-Q Residuals checks normality (points on the diagonal).
#     Both look acceptable for this model.
# Q3. Independence cannot be checked from the plots; it is judged from the study
#     design (each PHAA participant was surveyed once, not in clusters).

# =============================================================================
# SECTION 1, ACTIVITY 2: Logistic regression and its checks (Activity 4.2 on the lesson page)
# =============================================================================
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$hypertension <- factor(phaa$hypertension, levels = c("No", "Yes"))   # "Yes" is the event
table(phaa$hypertension)     # check: how many people had the event?

logit <- glm(hypertension ~ age + smoker, data = phaa, family = binomial)
summary(logit)
exp(cbind(OR = coef(logit), confint(logit)))   # odds ratios with 95% CIs

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Odds ratio for smoking and 95% CI:\n")
print(round(exp(cbind(OR = coef(logit), confint(logit)))["smokerYes", ], 2))
cat("\nQ2. Events and events per predictor if five predictors were used:\n")
n_events <- sum(phaa$hypertension == "Yes")
print(c(events = n_events, per_predictor_with_2 = n_events / 2, per_predictor_with_5 = n_events / 5))
# 23 events support about two predictors; five predictors would give fewer than
# 5 events each, which is too few.
# Q3. A Q-Q plot checks normal residuals, an assumption of linear regression
#     only. A logistic model has a binomial outcome and its residuals are not
#     expected to be normal.

# =============================================================================
# SECTION 2: Ordinal logistic regression for self-rated health (Activity 4.3 on the lesson page)
# =============================================================================
library(MASS)    # polr() fits the ordinal logistic model
library(brant)   # brant() checks the proportional-odds assumption
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$self_rated_health <- factor(phaa$self_rated_health,
  levels = c("Poor", "Fair", "Good", "Very good", "Excellent"), ordered = TRUE)
table(phaa$smoker, phaa$self_rated_health)   # check: people in each category, by smoking

ord <- polr(self_rated_health ~ age + smoker, data = phaa, Hess = TRUE)
summary(ord)

exp(cbind(OR = coef(ord), confint(ord)))   # odds ratios with 95% CIs

brant(ord)       # check the proportional-odds assumption

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Odds ratio for smoking (higher category) and 95% CI:\n")
print(round(exp(cbind(OR = coef(ord), confint(ord)))["smokerYes", ], 2))
cat("\nQ2. Odds ratio for a ten-year difference in age:\n")
print(round(exp(10 * coef(ord)[["age"]]), 2))
# Q3. Brant omnibus p = 0.35, age p = 0.23, smoking p = 0.53: all above 0.05,
#     so the proportional-odds assumption is reasonable.
cat("\nCrude odds ratios for smoking at each split (worked example):\n")
tab <- table(phaa$smoker, phaa$self_rated_health)
crude <- sapply(2:5, function(j) {
  above <- rowSums(tab[, j:5, drop = FALSE]); below <- rowSums(tab[, 1:(j - 1), drop = FALSE])
  (above["Yes"] / below["Yes"]) / (above["No"] / below["No"])
})
names(crude) <- paste("split", 1:4); print(round(crude, 2))

# =============================================================================
# SECTION 3: Multinomial logistic regression for usual place of care (Activity 4.4 on the lesson page)
# =============================================================================
library(nnet)    # multinom() fits the multinomial logistic model
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$usual_care <- relevel(factor(phaa$usual_care), ref = "Family doctor")   # reference category
table(phaa$usual_care, phaa$smoker)   # check: people in each category, by smoking

mn <- multinom(usual_care ~ age + smoker, data = phaa, trace = FALSE)
summary(mn)

round(exp(coef(mn)), 2)      # relative risk ratios (RRR)
round(exp(confint(mn)), 2)   # 95% confidence intervals for the RRRs

new_people <- data.frame(age = 30, smoker = c("No", "Yes"))   # two 30-year-olds
round(predict(mn, newdata = new_people, type = "probs"), 2)  # predicted probabilities

# --- Answers ---------------------------------------------------------------
cat("\nQ1. RRR for smoking, No usual place vs Family doctor, with 95% CI:\n")
print(round(c(RRR = exp(coef(mn)["No usual place", "smokerYes"]),
              exp(confint(mn))["smokerYes", , "No usual place"]), 2))
cat("\nQ2. RRR for ten years of age, by comparison:\n")
print(round(exp(10 * coef(mn)[, "age"]), 2))
# Q3. Age 30: non-smoker 0.51 family doctor, 0.28 walk-in, 0.06 ED, 0.15 none;
#     smoker 0.32, 0.30, 0.12, 0.26.
cat("\nCrude RRR for smoking, Emergency department vs Family doctor (worked example):\n")
t3 <- table(phaa$usual_care, phaa$smoker)
print(round((t3["Emergency department", "Yes"] / t3["Family doctor", "Yes"]) /
            (t3["Emergency department", "No"] / t3["Family doctor", "No"]), 2))

# =============================================================================
# SECTION 4, ACTIVITY 1: Poisson regression and the overdispersion check (Activity 4.5 on the lesson page)
# =============================================================================
library(MASS)    # glm.nb() fits the negative binomial model (next activity)
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
table(phaa$urgent_visits)    # how many people had 0, 1, 2, ... visits
mean(phaa$urgent_visits)     # the average count
var(phaa$urgent_visits)      # the variance (spread) of the counts

pois <- glm(urgent_visits ~ age + smoker, data = phaa, family = poisson)
summary(pois)
exp(cbind(IRR = coef(pois), confint(pois)))   # incidence rate ratios with 95% CIs

# Check 1: the dispersion ratio is close to 1 when the Poisson model fits
sum(residuals(pois, type = "pearson")^2) / df.residual(pois)
# Check 2: compare the zeros we observed with the zeros the model expects
sum(phaa$urgent_visits == 0)
sum(dpois(0, fitted(pois)))

# --- Answers ---------------------------------------------------------------
cat("\nQ1. IRR for smoking (Poisson) and 95% CI:\n")
print(round(exp(cbind(IRR = coef(pois), confint(pois)))["smokerYes", ], 2))
cat("\nQ2. Dispersion ratio, observed zeros, expected zeros:\n")
print(round(c(dispersion = sum(residuals(pois, type = "pearson")^2) / df.residual(pois),
              observed_zeros = sum(phaa$urgent_visits == 0),
              expected_zeros = sum(dpois(0, fitted(pois)))), 2))
# Q3. The point estimate may be reasonable, but the Poisson standard errors are
#     too small, so its intervals are too narrow and its p-values too small.
cat("\nCrude ratio of mean visits, smokers / non-smokers (worked example):\n")
m <- tapply(phaa$urgent_visits, phaa$smoker, mean); print(round(m, 3)); print(round(m["Yes"] / m["No"], 2))

# =============================================================================
# SECTION 4, ACTIVITY 2: Negative binomial regression and model comparison (Activity 4.6 on the lesson page)
# =============================================================================
nb <- glm.nb(urgent_visits ~ age + smoker, data = phaa)
summary(nb)
exp(cbind(IRR = coef(nb), confint(nb)))       # incidence rate ratios with 95% CIs

AIC(pois, nb)                                               # lower AIC = better fit
sum(residuals(nb, type = "pearson")^2) / df.residual(nb)    # dispersion ratio for NB
sum(dnbinom(0, mu = fitted(nb), size = nb$theta))           # zeros expected by NB

# --- Answers ---------------------------------------------------------------
cat("\nQ1. IRR for smoking: Poisson vs negative binomial:\n")
print(round(rbind(Poisson = exp(cbind(coef(pois), confint(pois)))["smokerYes", ],
                  NegBin  = exp(cbind(coef(nb),   confint(nb)))["smokerYes", ]), 2))
print(round(c(p_Poisson = summary(pois)$coefficients["smokerYes", 4],
              p_NegBin  = summary(nb)$coefficients["smokerYes", 4]), 4))
# Q2. AIC 1954.3 (NB) vs 2154.9 (Poisson); NB dispersion ratio 0.99; NB expects
#     about 463 zeros vs 464 observed. Report the negative binomial model.
# Q3. "In a negative binomial regression adjusted for age, smokers made 1.44
#     times as many emergency department or urgent care visits in the past year
#     as non-smokers (IRR 1.44, 95% CI 1.08 to 1.91)."

# =============================================================================
# SECTION 4, R EXAMPLE: an offset for GP visits over follow-up (Worked code 4.1 on the lesson page)
# =============================================================================
fu <- read.csv("phaa_followup.csv")      # GP visits over 0 to 10 years of follow-up
fu$smoker <- factor(fu$smoker, levels = c("No", "Yes"))
gp <- glm(gp_visits ~ age + smoker + offset(log(fu_years)), data = fu, family = poisson)
round(exp(cbind(IRR = coef(gp), confint(gp))), 3)