# =============================================================================
# HSCI 410 Lesson 5: Modelling Dependent Data -- ANSWER KEY
# Author: Kiffer G. Card, PhD - Faculty of Health Sciences, SFU
# -----------------------------------------------------------------------------
# Page:  modules/HSCI_410_Lesson_5_Modelling_Dependent_Data.html
# Data:  r-activities/data/phaa_clinics.csv   (sections 1 to 3)
#        r-activities/data/phaa_repeated.csv  (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: install once with
#   install.packages(c("lme4", "lmerTest", "geepack"))
#
# Both CSV files are simulated teaching data. The previous version of this
# lesson and its answer key are kept as
# HSCI_410_Lesson_5_Modelling_Dependent_Data_old_version.{html,R}.
# =============================================================================


# =============================================================================
# SECTION 1: Detecting clustering in the clinic data (Activity 5.1 on the lesson page)
# =============================================================================
# install.packages(c("lme4", "lmerTest", "geepack"))   # run once, if not yet installed
library(lmerTest)    # loads lme4 (for lmer()) and adds p-values
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
nlevels(clinics$clinic_id)                        # how many clinics?
summary(as.vector(table(clinics$clinic_id)))      # patients per clinic

m0 <- lmer(sbp ~ 1 + (1 | clinic_id), data = clinics)   # a model with clinics only
vc <- as.data.frame(VarCorr(m0))
vc[, c("grp", "vcov")]                    # variance between clinics and within clinics
icc <- vc$vcov[1] / sum(vc$vcov)          # intraclass correlation
icc

m_bar <- mean(table(clinics$clinic_id))   # average clinic size
deff  <- 1 + (m_bar - 1) * icc             # design effect
deff
nrow(clinics) / deff                       # effective sample size

naive <- lm(sbp ~ clinic_urban, data = clinics)                     # ignores clinics
mixed <- lmer(sbp ~ clinic_urban + (1 | clinic_id), data = clinics)  # allows for clinics
summary(naive)$coefficients
summary(mixed)$coefficients

# --- Answers ---------------------------------------------------------------
cat("\nQ1. ICC from the model with clinics only:\n")
print(round(icc, 3))
# About 17.5% of the variation in systolic blood pressure lies between clinics,
# so two patients from the same clinic have blood pressures correlated at about
# 0.17 and tend to be more alike than two patients from different clinics.
cat("\nQ2. Design effect and effective sample size:\n")
print(round(c(design_effect = deff, effective_n = nrow(clinics) / deff), 2))
# For clinic characteristics, the 966 patients carry about as much information
# as 150 independent people.
cat("\nQ3. Urban clinic: estimate, SE and p-value, ordinary vs mixed model:\n")
print(round(rbind(naive = summary(naive)$coefficients["clinic_urbanurban", c(1, 2, 4)],
                  mixed = summary(mixed)$coefficients["clinic_urbanurban", c(1, 2, 5)]), 3))
# Report the mixed model (1.67 mmHg, SE 1.96, p = 0.40): clinic_urban has only
# 30 independent values, and the ordinary regression (SE 0.74, p = 0.044)
# treats it as if it had 966.

# =============================================================================
# SECTION 2: A random-intercept model for blood pressure (Activity 5.2 on the lesson page)
# =============================================================================
library(lmerTest)   # loads lme4 and adds p-values to lmer() output
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
lmm <- lmer(sbp ~ age + female + smoker + clinic_urban + (1 | clinic_id), data = clinics)
summary(lmm)

confint(lmm, parm = "beta_", method = "Wald")   # 95% CIs for the fixed effects
vc <- as.data.frame(VarCorr(lmm))
vc$vcov[1] / sum(vc$vcov)                        # ICC after adjusting for the predictors

par(mfrow = c(1, 3))
plot(fitted(lmm), resid(lmm), main = "Residuals vs fitted"); abline(h = 0, lty = 2)
qqnorm(resid(lmm), main = "Q-Q: residuals"); qqline(resid(lmm))
qqnorm(ranef(lmm)$clinic_id[, 1], main = "Q-Q: clinic effects"); qqline(ranef(lmm)$clinic_id[, 1])
par(mfrow = c(1, 1))
isSingular(lmm)        # FALSE means the clinic variance was estimated without problems

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Fixed effect for smoking (mmHg) and 95% CI:\n")
print(round(c(estimate = fixef(lmm)[["smokerYes"]],
              confint(lmm, parm = "beta_", method = "Wald")["smokerYes", ]), 2))
# Smokers average 5.0 mmHg higher than non-smokers of the same age and gender
# attending the same clinic (95% CI 3.63 to 6.37).
cat("\nQ2. Variance components and adjusted ICC:\n")
print(round(c(clinic = vc$vcov[1], residual = vc$vcov[2], icc = vc$vcov[1] / sum(vc$vcov)), 3))
# The predictors explain much of the variation within clinics (residual
# variance falls from about 104 to 69) and little of the variation between
# clinics (about 22), so the clinic share rises from 0.175 to 0.24.
cat("\nQ3. Clinics with the largest estimated effects (mmHg):\n")
re <- ranef(lmm)$clinic_id[, 1]; names(re) <- rownames(ranef(lmm)$clinic_id)
print(round(sort(re, decreasing = TRUE)[1:3], 1))
cat("
Standard errors with and without allowing for clinics (Section 1, slide 10):
")
print(round(cbind(lm = summary(lm(sbp ~ age + female + smoker + clinic_urban, data = clinics))$coefficients[, 2],
                  lmer = summary(lmm)$coefficients[, 2]), 2))
# The residual plots support linearity, equal variance and normal residuals;
# the Q-Q plot of the clinic effects follows the line apart from these three
# clinics; isSingular() is FALSE; 30 clinics is enough. The model is acceptable.

# =============================================================================
# SECTION 3: A GLMM and GEE for specialist referral (Activity 5.3 on the lesson page)
# =============================================================================
library(lme4); library(geepack)   # glmer() for the GLMM, geeglm() for GEE
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
table(clinics$referred)                  # 1 = referred to a specialist
naive <- glm(referred ~ age + smoker + clinic_urban, family = binomial, data = clinics)
glmm  <- glmer(referred ~ age + smoker + clinic_urban + (1 | clinic_id),
               family = binomial, data = clinics)
summary(glmm)

clinics <- clinics[order(clinics$clinic_id), ]   # GEE needs each clinic's rows together
gee <- geeglm(referred ~ age + smoker + clinic_urban, id = clinic_id,
              family = binomial, corstr = "exchangeable", data = clinics)
summary(gee)
options(digits = 7)   # summary() of a GEE model lowers the printed digits; this restores the default

# Odds ratios and standard errors for urban clinics from the three models
se <- function(fit) summary(fit)$coefficients["clinic_urbanurban", 2]
b  <- c(naive = coef(naive)[["clinic_urbanurban"]], glmm = fixef(glmm)[["clinic_urbanurban"]],
        gee = coef(gee)[["clinic_urbanurban"]])
round(cbind(OR = exp(b), SE = c(se(naive), se(glmm), se(gee))), 3)

# --- Answers ---------------------------------------------------------------
# Q1. Urban clinic OR: GLMM 2.31 (clinic-specific: same clinic baseline, same
#     age and smoking status); GEE 2.04 (population-averaged: all urban-clinic
#     patients compared with all rural-clinic patients).
cat("\nQ2. Standard errors for urban clinic and for smoking in the three models:\n")
se2 <- function(fit, term) summary(fit)$coefficients[term, 2]
print(round(rbind(
  urban   = c(naive = se2(naive, "clinic_urbanurban"), glmm = se2(glmm, "clinic_urbanurban"), gee = se2(gee, "clinic_urbanurban")),
  smoking = c(naive = se2(naive, "smokerYes"),         glmm = se2(glmm, "smokerYes"),         gee = se2(gee, "smokerYes"))), 3))
# Urban clinic is a clinic-level predictor with information from 30 clinics,
# so its SE rises by about 70%; smoking varies within every clinic, so its SE
# barely changes.
cat("\nQ3. GEE odds ratio for urban clinics with 95% CI and p-value:\n")
print(round(c(exp(c(OR = coef(gee)[["clinic_urbanurban"]], confint.default(gee)["clinic_urbanurban", ])),
              p = summary(gee)$coefficients["clinic_urbanurban", 4]), 3))
cat("Number of clusters:", length(unique(clinics$clinic_id)), "\n")
# Report the GEE result for the provincial question (OR 2.04, p = 0.014) and
# note that 30 clinics is at the lower limit for GEE standard errors.

# =============================================================================
# SECTION 4, ACTIVITY 1: A mixed model for blood pressure over time (Activity 5.4 on the lesson page)
# =============================================================================
library(lmerTest)   # loads lme4 and adds p-values
visits <- read.csv("phaa_repeated.csv")   # file must be in the working directory
visits$id  <- factor(visits$id)
visits$arm <- factor(visits$arm, levels = c("control", "intervention"))
head(visits, 8)                                      # long format: one row per visit
table(visit = visits$visit, missing = is.na(visits$sbp_mmhg))
round(tapply(visits$sbp_mmhg, list(visits$arm, visits$visit), mean, na.rm = TRUE), 1)

lmm_t <- lmer(sbp_mmhg ~ arm * visit + (1 | id), data = visits)
summary(lmm_t)

vc <- as.data.frame(VarCorr(lmm_t))
vc$vcov[1] / sum(vc$vcov)                            # how alike one person's visits are
18 * fixef(lmm_t)[["armintervention:visit"]]          # extra change over 18 months

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Arm-by-visit interaction: estimate, SE and p-value:\n")
print(round(summary(lmm_t)$coefficients["armintervention:visit", c(1, 2, 5)], 4))
# Blood pressure fell 0.149 mmHg per month faster in the intervention arm, an
# extra fall of about 2.7 mmHg over 18 months (p = 0.014).
cat("\nQ2. Measurements used, people, and people with all four visits:\n")
complete <- tapply(!is.na(visits$sbp_mmhg), visits$id, all)
print(c(measurements = nobs(lmm_t), people = nlevels(visits$id), all_four_visits = sum(complete)))
# The mixed model uses every attended visit from all 200 people; it is valid if
# missed visits are missing at random (missingness depends only on observed
# information such as arm and earlier readings).
cat("\nQ3. ICC for repeated visits:\n")
print(round(vc$vcov[1] / sum(vc$vcov), 2))
# 0.56, compared with 0.175 for patients in clinics: one person's visits are
# much more alike than the patients of one clinic.

# =============================================================================
# SECTION 4, ACTIVITY 2: GEE for adherence over time (Activity 5.5 on the lesson page)
# =============================================================================
library(geepack)
adh <- visits[!is.na(visits$adherent), ]             # visits with adherence recorded
adh <- adh[order(adh$id, adh$visit), ]               # each person's rows together
round(tapply(adh$adherent, list(adh$arm, adh$visit), mean), 2)   # share adherent
gee_a <- geeglm(adherent ~ arm * visit, id = id, family = binomial,
                corstr = "exchangeable", data = adh)
summary(gee_a)
options(digits = 7)   # restore the default number of printed digits

round(exp(cbind(OR = coef(gee_a), confint.default(gee_a))), 3)   # odds ratios, 95% CIs

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Interaction odds ratio per month, and over six months:\n")
or_int <- exp(coef(gee_a)[["armintervention:visit"]])
print(round(c(per_month = or_int, per_six_months = or_int^6), 3))
cat("\nQ2. Number of clusters (people) in the GEE model:\n")
print(length(unique(adh$id)))
# geeglm() treats consecutive rows with the same id as one cluster, so the data
# are sorted by id (and visit) before fitting.
cat("\nQ3. Estimated exchangeable correlation (alpha):\n")
print(round(gee_a$geese$alpha, 3))
# Even a small correlation breaks independence; GEE allows for it at little
# cost, so the GEE result is the one to report.