# ========================================================
# MEDIATION AND MODERATION IN R
# Canadian Social Connection Survey (CSCS), 2021 wave
# ========================================================
# Question 1 (mediation): does loneliness help explain why
# people with less social support report more depressive
# symptoms?
# Question 2 (moderation): is the link between social
# support and loneliness the same for younger and older
# adults?
# ========================================================


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

# install.packages("mediation")   # run once
library(mediation)   # mediate() tests the indirect effect

github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData"
load(url(github_url))

# Keep the 2021 wave, which gives one row per person
data <- data[data$SURVEY_collection_year == 2021, ]
dim(data)


# ========================================================
# SECTION 2: PREPARE THE VARIABLES
# ========================================================
# support:    Multidimensional Scale of Perceived Social
#             Support, 1 to 7 (higher = more support)
# loneliness: UCLA 3-item Loneliness Scale, 3 to 9
#             (higher = lonelier)
# depression: PHQ-2 depressive symptoms, 0 to 6
#             (higher = more symptoms)
med_data <- data.frame(
  support    = data$PSYCH_zimet_multidimensional_social_support_scale_score,
  loneliness = data$LONELY_ucla_loneliness_scale_score,
  depression = data$WELLNESS_phq_score,
  age        = data$DEMO_age)
summary(med_data)

# Every model must use the same people, so keep the people
# with no missing values
med_data <- na.omit(med_data)
nrow(med_data)

# Correlations between the three main variables
round(cor(med_data[, c("support", "loneliness", "depression")]), 2)


# ========================================================
# SECTION 3: MEDIATION, ONE PATH AT A TIME
# ========================================================
# X = support, M = loneliness (the mediator), Y = depression

# --- 3.1 Path c, the total effect: X -> Y
model_total <- lm(depression ~ support, data = med_data)
summary(model_total)

# --- 3.2 Path a: X -> M
model_a <- lm(loneliness ~ support, data = med_data)
summary(model_a)

# --- 3.3 Path b (M -> Y) and path c' (the direct effect),
# estimated together in one model
model_b <- lm(depression ~ support + loneliness, data = med_data)
summary(model_b)

# --- 3.4 The indirect effect is path a times path b
a <- coef(model_a)["support"]
b <- coef(model_b)["loneliness"]
a * b

# Check: total effect = direct effect + indirect effect
coef(model_b)["support"] + a * b
coef(model_total)["support"]


# ========================================================
# SECTION 4: TEST THE INDIRECT EFFECT
# ========================================================
# mediate() takes the two models (path a, then paths b and c')
# and uses the bootstrap: it resamples the data 1,000 times
# to build a confidence interval for the indirect effect.
set.seed(2021)    # makes the bootstrap results repeatable
med_result <- mediate(model_a, model_b,
                      treat = "support", mediator = "loneliness",
                      boot = TRUE, sims = 1000)
summary(med_result)

# A plot of the three effects with their 95% confidence intervals
plot(med_result, xlim = c(-0.7, 0),
     main = "Indirect (ACME), direct (ADE) and total effects")


# ========================================================
# SECTION 5: MODERATION WITH AN INTERACTION TERM
# ========================================================

# --- 5.1 The moderator: two age groups
med_data$age_group <- NA
med_data$age_group[med_data$age < 50] <- "Under 50"
med_data$age_group[med_data$age >= 50] <- "50 and over"
med_data$age_group <- factor(med_data$age_group,
                             levels = c("Under 50", "50 and over"))
table(med_data$age_group)

# --- 5.2 Without and with the interaction term
# support * age_group is shorthand for
# support + age_group + support:age_group
model_main <- lm(loneliness ~ support + age_group, data = med_data)
model_int  <- lm(loneliness ~ support * age_group, data = med_data)
summary(model_int)

# Does adding the interaction improve the model?
anova(model_main, model_int)

# --- 5.3 The slope of support in each age group
# Under 50 (the reference group): the support coefficient
coef(model_int)["support"]
# 50 and over: the support coefficient plus the interaction
coef(model_int)["support"] + coef(model_int)["support:age_group50 and over"]

# --- 5.4 Plot one line for each age group
plot(jitter(med_data$support), jitter(med_data$loneliness),
     main = "Social support and loneliness, by age group",
     xlab = "Social support (1 to 7)",
     ylab = "UCLA loneliness score (3 to 9)",
     pch = 16, col = rgb(0, 0, 0, 0.10))
abline(lm(loneliness ~ support,
          data = med_data[med_data$age_group == "Under 50", ]),
       col = "steelblue", lwd = 3)
abline(lm(loneliness ~ support,
          data = med_data[med_data$age_group == "50 and over", ]),
       col = "firebrick", lwd = 3)
legend("topright", legend = c("Under 50", "50 and over"),
       col = c("steelblue", "firebrick"), lwd = 3)