# ========================================================
# MIXED MODELS AND GEE IN R
# Patients in 30 clinics, and a trial with repeated visits
# ========================================================
# Questions: once patients in the same clinic are allowed
# to be alike, do patients of urban and rural clinics
# differ in blood pressure and in specialist referral?
# In a wellness trial with four visits, did blood pressure
# and adherence change at different rates in the two arms?
# ========================================================


# ========================================================
# SECTION 1: PACKAGES AND THE CLINIC DATA
# ========================================================

# install.packages(c("lme4", "lmerTest", "geepack"))   # run once
library(lmerTest)   # loads lme4 (lmer, glmer) and adds p-values
library(geepack)    # geeglm() fits GEE models

# The data files are stored on the course website
clinics_url <- "https://www.sfu-epi.ca/r-activities/data/phaa_clinics.csv"
clinics <- read.csv(clinics_url)
dim(clinics)
head(clinics)

# Factors: the first level is the reference group
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"))

# How many clinics, and how many patients in each?
nlevels(clinics$clinic_id)
summary(as.vector(table(clinics$clinic_id)))

# Blood pressure in each clinic, sorted by the clinic average
clinics$clinic_sorted <- reorder(clinics$clinic_id,
                                 clinics$sbp, FUN = mean)
boxplot(sbp ~ clinic_sorted, data = clinics, col = "grey85",
        las = 2, cex.axis = 0.7,
        xlab = "Clinic, sorted by its average",
        ylab = "Systolic blood pressure (mmHg)",
        main = "Blood pressure in 30 clinics")
abline(h = mean(clinics$sbp), lty = 2)   # overall average


# ========================================================
# SECTION 2: HOW STRONG IS THE CLUSTERING?
# ========================================================

# --- 2.1 A mixed model with clinics only, no predictors
m0 <- lmer(sbp ~ 1 + (1 | clinic_id), data = clinics)

# Variance between clinics and variance within clinics
vc <- as.data.frame(VarCorr(m0))
vc[, c("grp", "vcov")]

# The intraclass correlation: the share between clinics
icc <- vc$vcov[1] / sum(vc$vcov)
icc

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


# ========================================================
# SECTION 3: A CLINIC-LEVEL PREDICTOR, TWO WAYS
# ========================================================

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


# ========================================================
# SECTION 4: A RANDOM-INTERCEPT MODEL FOR BLOOD PRESSURE
# ========================================================

lmm <- lmer(sbp ~ age + female + smoker + clinic_urban +
              (1 | clinic_id), data = clinics)
summary(lmm)

# 95% confidence intervals for the fixed effects
confint(lmm, parm = "beta_", method = "Wald")

# The ICC after adjusting for the predictors
vc <- as.data.frame(VarCorr(lmm))
vc$vcov[1] / sum(vc$vcov)


# ========================================================
# SECTION 5: CHECK THE MIXED MODEL
# ========================================================

# --- 5.1 Residuals against fitted values
plot(fitted(lmm), resid(lmm),
     xlab = "Fitted value (mmHg)", ylab = "Residual (mmHg)",
     main = "Residuals vs fitted")
abline(h = 0, lty = 2)

# --- 5.2 Q-Q plots: residuals, then the 30 clinic effects
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])

# --- 5.3 FALSE means the clinic variance was estimated
isSingular(lmm)

# A random slope would also let the effect of age differ
# between clinics: (age | clinic_id) in place of
# (1 | clinic_id). This script fits random intercepts only.


# ========================================================
# SECTION 6: A BINARY OUTCOME WITH A GLMM
# ========================================================

table(clinics$referred)   # 1 = referred to a specialist

# --- 6.1 Ordinary logistic regression, ignoring clinics
naive <- glm(referred ~ age + smoker + clinic_urban,
             family = binomial, data = clinics)

# --- 6.2 Logistic regression with a random intercept
glmm <- glmer(referred ~ age + smoker + clinic_urban +
                (1 | clinic_id), family = binomial, data = clinics)
summary(glmm)

# Clinic-specific odds ratios with 95% CIs
round(exp(cbind(OR = fixef(glmm),
                confint(glmm, parm = "beta_", method = "Wald"))), 3)


# ========================================================
# SECTION 7: THE SAME OUTCOME WITH GEE
# ========================================================

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

# Population-averaged odds ratios with 95% CIs
round(exp(cbind(OR = coef(gee), confint.default(gee))), 3)

# --- 7.1 Urban clinics in the three models: OR and SE
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)


# ========================================================
# SECTION 8: REPEATED MEASURES IN A WELLNESS TRIAL
# ========================================================

visits_url <- "https://www.sfu-epi.ca/r-activities/data/phaa_repeated.csv"
visits <- read.csv(visits_url)
visits$id <- factor(visits$id)
visits$arm <- factor(visits$arm,
                     levels = c("control", "intervention"))
head(visits, 8)   # long format: one row per visit

# --- 8.1 Missed visits
table(visit = visits$visit, missing = is.na(visits$sbp_mmhg))
complete <- tapply(!is.na(visits$sbp_mmhg), visits$id, all)
sum(complete)     # people who attended all four visits

# --- 8.2 Average blood pressure in each arm at each visit
sbp_means <- round(tapply(visits$sbp_mmhg,
                          list(visits$arm, visits$visit),
                          mean, na.rm = TRUE), 1)
sbp_means

plot(c(0, 6, 12, 18), sbp_means["control", ], type = "b",
     pch = 16, lwd = 2, col = "steelblue", ylim = c(122, 132),
     xlab = "Months since the start of the trial",
     ylab = "Average blood pressure (mmHg)",
     main = "Blood pressure by arm and visit")
lines(c(0, 6, 12, 18), sbp_means["intervention", ], type = "b",
      pch = 16, lwd = 2, col = "firebrick")
legend("topright", legend = c("Control", "Intervention"),
       col = c("steelblue", "firebrick"), lwd = 2, pch = 16)


# ========================================================
# SECTION 9: A MIXED MODEL FOR CHANGE OVER TIME
# ========================================================

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

# How alike are one person's visits?
vc <- as.data.frame(VarCorr(lmm_t))
vc$vcov[1] / sum(vc$vcov)

# The extra change in the intervention arm over 18 months
18 * fixef(lmm_t)[["armintervention:visit"]]


# ========================================================
# SECTION 10: GEE FOR ADHERENCE OVER TIME
# ========================================================

adh <- visits[!is.na(visits$adherent), ]   # adherence recorded
adh <- adh[order(adh$id, adh$visit), ]     # each person together
round(tapply(adh$adherent, list(adh$arm, adh$visit), mean), 2)

gee_a <- geeglm(adherent ~ arm * visit, id = id,
                family = binomial, corstr = "exchangeable",
                data = adh)
summary(gee_a)
options(digits = 7)   # restore the default printed digits

# Odds ratios with 95% CIs, then the interaction over 6 months
round(exp(cbind(OR = coef(gee_a), confint.default(gee_a))), 3)
exp(6 * coef(gee_a)[["armintervention:visit"]])