# =============================================================================
# HSCI 410 Public Health Assessment and Analysis  -  Lesson 6: Modelling Dependent Data
# Answer key for the in-lesson R activities
# Data file(s): phaa_clinics.csv, phaa_repeated.csv  (download from the lesson page; save in your R working directory)
# Packages: lme4, lmerTest, performance, sandwich, lmtest, geepack, nlme
#           (install once with install.packages(c("lme4", "lmerTest", "performance", "sandwich", "lmtest", "geepack", "nlme")))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================
# Parts 1 to 3 use phaa_clinics.csv: 966 patients nested in 30 primary-care
# clinics (18 to 45 patients each), with a continuous outcome (sbp) and a binary
# outcome (referred). Part 4 uses phaa_repeated.csv: 200 trial participants
# measured at months 0, 6, 12 and 18 (sbp_mmhg and adherent), with about 12%
# of visits missing.

library(lme4);  library(lmerTest);  library(performance)
library(sandwich);  library(lmtest);  library(geepack);  library(nlme)
# NOTE (answer key): recent lme4 releases print variance components with the
# `%||%` operator, which base R only gained in version 4.4. On R 4.3 that makes
# summary() of any lmer()/glmer() fit stop with "could not find function
# '%||%'". Defining it once here is harmless on newer R and fixes older R.
if (!exists("%||%")) `%||%` <- function(x, y) if (is.null(x)) y else x

# =============================================================================
# ==== Part 1: Introduction to Clustered Data ====
# =============================================================================

# ==== Section 2: Effects of Clustering on Statistical Analysis / Activity: ICC and design effect from a clustered dataset ====
clinics <- read.csv("phaa_clinics.csv", stringsAsFactors = FALSE,
                    na.strings = c("", "NA"))
clinics$clinic_id <- factor(clinics$clinic_id)

table(clinics$clinic_id)                      # 18 to 45 patients per clinic
clinic_means <- aggregate(sbp ~ clinic_id, data = clinics, FUN = mean)
range(clinic_means$sbp)                       # clinic means spread from 102 to 126 mmHg

# 1. ICC from a one-way ANOVA
aov_fit <- aov(sbp ~ clinic_id, data = clinics)
ms     <- summary(aov_fit)[[1]][, "Mean Sq"]
n_per  <- mean(table(clinics$clinic_id))
icc_h  <- (ms[1] - ms[2]) / (ms[1] + (n_per - 1) * ms[2])
icc_h

# 2. Same ICC, conveniently, from a null mixed model
m_null <- lmer(sbp ~ 1 + (1 | clinic_id), data = clinics)
icc(m_null)

# 3. Design effect
deff <- 1 + (n_per - 1) * icc_h;  deff
nrow(clinics) / deff                          # effective sample size

# 4. Cluster-robust SEs as a quick fix for naive OLS
naive <- lm(sbp ~ age + smoker + bmi + clinic_urban, data = clinics)
summary(naive)$coef[, "Std. Error"]           # naive (independence) SEs
coeftest(naive, vcov. = vcovCL, cluster = ~ clinic_id)

# ---- Activity questions (Part 1, box rActivity-410-9-1) ----------------------
icc_null <- icc(m_null)$ICC_adjusted
cat("\nQ1: icc_h (one-way ANOVA) =", round(icc_h, 3),
    "  icc(m_null) =", round(icc_null, 3), "\n")
cat("    Essentially the same (the ANOVA formula is the method-of-moments version\n",
    "   of the REML estimate): about 17% of the total variance in SBP lies between\n",
    "   clinics and 83% between patients within clinics.\n")

cat("\nQ2: average cluster size =", round(n_per, 1), " deff = 1 + (", round(n_per, 1),
    "- 1) x", round(icc_h, 3), "=", round(deff, 2), "\n")
cat("    effective sample size = 966 /", round(deff, 2), "=", round(nrow(clinics) / deff), "\n")
cat("    Each clustered observation carries the information of about 1/", round(deff, 1),
    " = ", round(1 / deff, 2), " of an independent observation, so the 966 patients\n",
    "    are worth roughly ", round(nrow(clinics) / deff), " independent ones.\n", sep = "")

se_naive  <- summary(naive)$coef[, "Std. Error"]
se_robust <- coeftest(naive, vcov. = vcovCL, cluster = ~ clinic_id)[, "Std. Error"]
cat("\nQ3: naive vs cluster-robust SEs\n")
print(round(cbind(naive = se_naive, robust = se_robust, ratio = se_robust / se_naive), 3))
cat("    The clinic-level predictor clinic_urban sees by far the biggest inflation\n",
    "   (about 3.4 times); bmi and smoker grow by 20-25%, and the age SE barely moves.\n",
    "   Naive OLS treats the 966 rows as independent, but patients in the same clinic\n",
    "   share a clinic effect, so each row carries less new information than OLS\n",
    "   assumes; for a variable that is constant within a clinic there are really\n",
    "   only 30 independent values, which is why its SE is understated the most.\n")

# =============================================================================
# ==== Part 2: Mixed Models for Continuous Data ====
# =============================================================================

# ==== Section 5: Introduction & The Linear Mixed Model / Activity: Random-intercept linear mixed model with lme4::lmer() ====
clinics <- read.csv("phaa_clinics.csv", stringsAsFactors = FALSE,
                    na.strings = c("", "NA"))
clinics$clinic_id    <- factor(clinics$clinic_id)
clinics$clinic_urban <- factor(clinics$clinic_urban,
                                  levels = c("rural","urban"))
clinics$smoker       <- factor(clinics$smoker, levels = c("No","Yes"))

# 1. Null model: variance partitioning
m0 <- lmer(sbp ~ 1 + (1 | clinic_id), data = clinics)
summary(m0)                   # random-effects table: between-clinic and residual variance
icc(m0)

# 2. Add patient-level fixed effects
m1 <- lmer(sbp ~ age + smoker + bmi + female + (1 | clinic_id),
           data = clinics)
summary(m1)
confint(m1, method = "Wald")

# 3. Add cluster-level fixed effects (clinic_urban, clinic_size)
m2 <- lmer(sbp ~ age + smoker + bmi + female
                  + clinic_urban + scale(clinic_size)
                  + (1 | clinic_id),
           data = clinics)
summary(m2)
icc(m2)

# 4. Is the random intercept needed?
ranova(m2)

# 5. Random slope: does the smoker effect vary by clinic?
m3 <- lmer(sbp ~ age + smoker + bmi + female + clinic_urban + scale(clinic_size)
                  + (1 + smoker | clinic_id), data = clinics)
anova(m2, m3)        # LRT: random-intercept-only vs random-slope

# Contextual effects (mentioned in the lesson note): split age into the clinic
# mean and the within-clinic deviation
clinics$age_clinic_mean <- ave(clinics$age, clinics$clinic_id, FUN = mean)
clinics$age_within      <- clinics$age - clinics$age_clinic_mean
m4 <- lmer(sbp ~ age_within + age_clinic_mean + smoker + bmi + female
                  + clinic_urban + (1 | clinic_id), data = clinics)
round(fixef(m4)[c("age_within", "age_clinic_mean")], 3)

# ---- Activity questions (Part 2, box rActivity-410-10-1) ---------------------
vc0 <- as.data.frame(VarCorr(m0))
s2u <- vc0$vcov[vc0$grp == "clinic_id"];  s2e <- vc0$vcov[vc0$grp == "Residual"]
cat("\nQ1: m0 between-clinic variance sigma^2_u =", round(s2u, 1),
    " residual variance sigma^2 =", round(s2e, 1), "\n")
cat("    ICC by hand =", round(s2u, 1), "/ (", round(s2u, 1), "+", round(s2e, 1), ") =",
    round(s2u / (s2u + s2e), 3), "  icc(m0) =", round(icc(m0)$ICC_adjusted, 3), "\n")
cat("    About 17% of the total SBP variance lies between clinics; the other 83% is\n",
    "   between patients within a clinic. Note the scale: SD 4.7 mmHg between clinics\n",
    "   against 10.2 mmHg within.\n")

vc2 <- as.data.frame(VarCorr(m2))
cat("\nQ2: icc(m0) adjusted =", round(icc(m0)$ICC_adjusted, 3),
    "  icc(m2) adjusted =", round(icc(m2)$ICC_adjusted, 3),
    " unadjusted =", round(icc(m2)$ICC_unadjusted, 3), "\n")
cat("    between-clinic variance: m0 =", round(s2u, 1), " m2 =",
    round(vc2$vcov[vc2$grp == "clinic_id"], 1),
    "  residual variance: m0 =", round(s2e, 1), " m2 =",
    round(vc2$vcov[vc2$grp == "Residual"], 1), "\n")
cat("    The between-clinic variance hardly shrank (22.1 to 21.6), so clinic_urban and\n",
    "   clinic_size explain almost none of the clustering (neither is significant).\n",
    "   The patient-level covariates removed a large share of the WITHIN-clinic\n",
    "   variance, which is why the adjusted ICC rises to 0.25: the clinic share of\n",
    "   what remains unexplained is larger. Clinics still differ for reasons the\n",
    "   model does not measure.\n")

lrt <- anova(m2, m3)
cat("\nQ3: random slope on smoker: chi-square =", round(lrt$Chisq[2], 2),
    "on", lrt$Df[2], "df, p =", round(lrt$`Pr(>Chisq)`[2], 2), "\n")
cat("    Not significant, so keep the random-intercept-only model m2. There is no\n",
    "   evidence that the smoking effect (about +5 mmHg) differs between clinics,\n",
    "   and the two extra parameters (slope variance and its correlation with the\n",
    "   intercept) do not improve fit. The simpler model is more stable, easier to\n",
    "   report, and less prone to fitting clinic-level noise.\n")

# =============================================================================
# ==== Part 3: Mixed Models for Discrete Data ====
# =============================================================================

# ==== Section 10: GLMMs for Count, Binary & Categorical Data / Activity: Logistic GLMM and GEE on the same clustered data ====
clinics <- read.csv("phaa_clinics.csv", stringsAsFactors = FALSE,
                    na.strings = c("", "NA"))
clinics$clinic_id    <- factor(clinics$clinic_id)
clinics$clinic_urban <- factor(clinics$clinic_urban,
                                  levels = c("rural","urban"))
clinics$smoker       <- factor(clinics$smoker, levels = c("No","Yes"))

# 1. Logistic GLMM (subject-specific)
m_glmm <- glmer(referred ~ age + female + smoker + bmi + clinic_urban
                            + (1 | clinic_id),
                data    = clinics,
                family  = binomial,
                control = glmerControl(optimizer = "bobyqa"))
summary(m_glmm)
exp(fixef(m_glmm))                     # subject-specific ORs
exp(confint(m_glmm, method = "Wald"))  # Wald 95% CIs for the ORs
icc(m_glmm)                            # latent-scale ICC

# 2. GEE (population-averaged)
m_gee <- geeglm(referred ~ age + female + smoker + bmi + clinic_urban,
                id     = clinic_id, data = clinics,
                family = binomial, corstr = "exchangeable")
summary(m_gee)
options(digits = 7)   # geepack's summary() print method leaves options(digits) at 4; restore it
exp(coef(m_gee))                       # population-averaged ORs

# 3. Compare side-by-side
cbind(GLMM_OR = exp(fixef(m_glmm)),
      GEE_OR  = exp(coef(m_gee)))

# ---- Activity questions (Part 3, box rActivity-410-11-1) ---------------------
or_glmm <- exp(fixef(m_glmm))
ci_glmm <- exp(confint(m_glmm, method = "Wald"))
cat("\nQ1: subject-specific OR for smokerYes =", round(or_glmm["smokerYes"], 2),
    " Wald 95% CI", round(ci_glmm["smokerYes", 1], 2), "to", round(ci_glmm["smokerYes", 2], 2), "\n")
cat("    Within the same clinic, a patient who smokes has about 1.7 times the odds of\n",
    "   a specialist referral of a non-smoking patient of the same age, sex and BMI.\n",
    "   The comparison holds the clinic's random intercept fixed, which is what makes\n",
    "   it conditional (subject-specific).\n")

s2u_bin <- as.data.frame(VarCorr(m_glmm))$vcov[1]
cat("\nQ2: latent-scale ICC from icc(m_glmm) =", round(icc(m_glmm)$ICC_adjusted, 3),
    " (by hand:", round(s2u_bin, 3), "/ (", round(s2u_bin, 3), "+ pi^2/3) =",
    round(s2u_bin / (s2u_bin + pi^2 / 3), 3), ")\n")
cat("    A binary outcome has no free residual variance: on the logistic latent\n",
    "   scale the within-cluster variance is fixed at pi^2/3 = 3.29 by the link, so the\n",
    "   ICC is sigma^2_u / (sigma^2_u + 3.29) and depends only on the between-clinic\n",
    "   variance. About 7% of the latent variation in referral is between clinics.\n")

cmp <- cbind(GLMM_OR = exp(fixef(m_glmm)), GEE_OR = exp(coef(m_gee)))
cat("\nQ3: GLMM vs GEE odds ratios\n")
print(round(cmp, 3))
cat("    The GLMM (conditional) ORs are the larger ones for every predictor, for\n",
    "   example smokerYes 1.71 vs 1.68 and clinic_urban 2.36 vs 2.09. With a logit\n",
    "   link, averaging subject-specific probabilities over the clinic random effects\n",
    "   pulls the marginal OR toward 1 by a factor of roughly 1/sqrt(1 + 0.35 x\n",
    "   sigma^2_u); here sigma^2_u is only 0.26, so the attenuation is small.\n",
    "   A public-health audience prefers the GEE (population-averaged) OR when the\n",
    "   question is what would change across the whole population of clinics; the\n",
    "   GLMM OR answers the within-clinic, patient-level question.\n")

# =============================================================================
# ==== Part 4: Repeated Measures ====
# =============================================================================

# ==== Section 14: Univariate & Multivariate Approaches / Activity: Longitudinal mixed model with autoregressive errors ====
dat <- read.csv("phaa_repeated.csv", stringsAsFactors = FALSE,
                na.strings = c("", "NA"))
dat$id  <- factor(dat$id)
dat$arm <- factor(dat$arm, levels = c("control","intervention"))

table(table(dat$id))              # 4 rows per person
sum(is.na(dat$sbp_mmhg))          # 96 missing visits (12%)

# 1. Random-intercept mixed model on the continuous outcome
m_lmm <- lmer(sbp_mmhg ~ visit * arm + age + female + (1 | id),
              data = dat)
summary(m_lmm)

# 2. Add an AR(1) within-subject correlation structure with nlme
m_lme <- lme(sbp_mmhg ~ visit * arm + age + female,
             random      = ~ 1 | id,
             correlation = corAR1(form = ~ visit | id),
             data        = dat, na.action = na.omit)
summary(m_lme)

# 3. Compare correlation structures with AIC
m_cs <- update(m_lme, correlation = corCompSymm(form = ~ visit | id))
m_un <- update(m_lme, correlation = corSymm(form    = ~ 1     | id))
AIC(m_lme, m_cs, m_un)

# 4. Same trial, binary outcome (adherence) -- GLMM and GEE
m_bin <- glmer(adherent ~ visit * arm + age + female + (1 | id),
               data    = dat, family = binomial,
               control = glmerControl(optimizer = "bobyqa"))
exp(fixef(m_bin))               # subject-specific ORs

m_gee <- geeglm(adherent ~ visit * arm + age + female,
                id = id, data = dat, family = binomial,
                corstr = "exchangeable")
summary(m_gee)
options(digits = 7)             # restore the digits option after geepack's summary()
exp(coef(m_gee))                # population-averaged ORs

# ---- Activity questions (Part 4, box rActivity-410-12-1) ---------------------
cf <- summary(m_lmm)$coefficients
cat("\nQ1: visit:armintervention =", round(cf["visit:armintervention", "Estimate"], 3),
    " mmHg per month, SE =", round(cf["visit:armintervention", "Std. Error"], 3),
    ", t =", round(cf["visit:armintervention", "t value"], 2),
    ", p =", round(cf["visit:armintervention", "Pr(>|t|)"], 3), "\n")
cat("    Each month, SBP in the intervention arm changes by about 0.15 mmHg more\n",
    "   (more negative) than in the control arm, on top of the control drift of\n",
    "   -0.16 mmHg per month. Over the 18-month trial that is roughly 2.6 mmHg\n",
    "   more reduction than control, a modest but statistically significant effect.\n")

aics <- AIC(m_lme, m_cs, m_un)
cat("\nQ2: AIC  AR(1) =", round(aics["m_lme", "AIC"], 1), "  CS =", round(aics["m_cs", "AIC"], 1),
    "  UN =", round(aics["m_un", "AIC"], 1), "\n")
cat("    AR(1) phi =", round(as.numeric(coef(m_lme$modelStruct$corStruct, unconstrained = FALSE)), 3),
    "  CS rho =", round(as.numeric(coef(m_cs$modelStruct$corStruct, unconstrained = FALSE)), 3), "\n")
cat("    AR(1) and compound symmetry tie exactly, because both estimate a residual\n",
    "   correlation of 0: the random intercept already captures all the within-person\n",
    "   correlation. Unstructured is 8.4 units worse with 5 extra parameters. Report\n",
    "   the simplest model, the random intercept alone (equivalently CS), and mention\n",
    "   that AR(1) and UN change nothing.\n")

or_bin <- exp(fixef(m_bin));  or_gee <- exp(coef(m_gee))
cat("\nQ3: OR for visit:armintervention (per month): GLMM =", round(or_bin["visit:armintervention"], 3),
    "  GEE =", round(or_gee["visit:armintervention"], 3), "\n")
cat("    The subject-specific (GLMM) OR is the larger one, 1.100 vs 1.095, but only\n",
    "   just: the between-person variance on the logit scale is small (0.21), so\n",
    "   averaging over people barely attenuates the effect. With a logit link the\n",
    "   population-averaged OR is always closer to 1 than the conditional OR, and the\n",
    "   gap grows with the random-effect variance. Within a person, each month in the\n",
    "   intervention arm raises the odds of adherence by about 10% relative to control.\n")