# ============================================================================= # HSCI 410 Lesson 8: Mediation, Moderation and Path Analysis -- ANSWER KEY # Author: Kiffer G. Card, PhD - Faculty of Health Sciences, SFU # ----------------------------------------------------------------------------- # Page: modules/HSCI_410_Lesson_8_Mediation_Moderation_and_Path_Analysis.html # Data: Canadian Social Connection Survey (CSCS), 2021 wave, loaded from GitHub # with load(url(github_url)), as in the narrated R walkthroughs # r-walkthroughs/HSCI_410_Mediation_and_Moderation.html and # r-walkthroughs/HSCI_410_SEM_Path_Analysis.html # # This script reproduces every code block on the page, in page order, and then # answers each numbered activity question with a printed result. Each activity # loads the data again, as on the page, so that it can be run on its own. # # Packages: install once with # install.packages(c("mediation", "lavaan")) # # Mediation estimates from these cross-sectional data are statistical # decompositions of associations. They do not establish causal pathways. # ============================================================================= options(width = 70) # ============================================================================= # SECTION 1: Load the CSCS data and meet the three variables (Activity 8.1 on the lesson page) # ============================================================================= # Load the CSCS data from GitHub (this takes a few seconds) github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave dim(data) # people (rows) and variables (columns) # support: social support (MSPSS), 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) med_data <- na.omit(med_data) # keep people with no missing values nrow(med_data) # how many people are in the analysis? round(cor(med_data[, c("support", "loneliness", "depression")]), 2) # --- Answers --------------------------------------------------------------- cat("\nQ1. People in the analysis after na.omit():", nrow(med_data), "\n") # Every model must use the same people, so that the total effect equals the # direct effect plus the indirect effect. cat("\nQ2. Correlations:\n") print(round(cor(med_data[, c("support", "loneliness", "depression")]), 2)) # support-loneliness -0.40, support-depression -0.38, loneliness-depression 0.47. # Q3. Correlations are symmetric and the three measures come from one survey # wave, so they cannot show which variable came first. Longitudinal data # (support, then loneliness, then symptoms) would establish temporal order. # ============================================================================= # SECTION 2, ACTIVITY 1: Moderation by age group (Activity 8.2 on the lesson page) # ============================================================================= github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave 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) med_data <- na.omit(med_data) # the same people in every model # The moderator: two age groups, with "Under 50" as the reference 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) # 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) anova(model_main, model_int) # does the interaction improve the model? # Simple slopes: the slope of support in each age group coef(model_int)["support"] # Under 50 coef(model_int)["support"] + coef(model_int)["support:age_group50 and over"] # 50 and over # One fitted line for each age group plot(jitter(med_data$support), jitter(med_data$loneliness), 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) # --- Answers --------------------------------------------------------------- cat("\nQ1. Interaction coefficient, its p-value and 95% CI:\n") print(round(summary(model_int)$coefficients["support:age_group50 and over", ], 6)) print(round(confint(model_int)["support:age_group50 and over", ], 3)) cat("\nQ2. Simple slopes of support by age group:\n") print(round(c(under_50 = unname(coef(model_int)["support"]), age_50_plus = unname(coef(model_int)["support"] + coef(model_int)["support:age_group50 and over"])), 3)) cat("\nQ3. Difference between age groups at support scores of 1, 5 and 7:\n") s <- c(1, 5, 7) print(round(setNames(coef(model_int)["age_group50 and over"] + coef(model_int)["support:age_group50 and over"] * s, paste("support", s)), 2)) # The 1.83 is the difference at support = 0, outside the 1 to 7 scale. # ============================================================================= # SECTION 2, ACTIVITY 2: Centring a continuous moderator and simple slopes (Activity 8.3 on the lesson page) # ============================================================================= github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave 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) med_data <- na.omit(med_data) # the same people in every model # Centre each continuous variable by subtracting its mean med_data$age_c <- med_data$age - mean(med_data$age) med_data$support_c <- med_data$support - mean(med_data$support) round(c(mean_age = mean(med_data$age), mean_support = mean(med_data$support)), 2) model_raw <- lm(loneliness ~ support * age, data = med_data) # uncentred model_cen <- lm(loneliness ~ support_c * age_c, data = med_data) # centred round(summary(model_raw)$coefficients[, 1:3], 4) # estimate, SE, t value round(summary(model_cen)$coefficients[, 1:3], 4) # Simple slopes of support at ages 30, 50 and 70 b <- coef(model_cen) ages <- c(30, 50, 70) slopes <- b["support_c"] + b["support_c:age_c"] * (ages - mean(med_data$age)) names(slopes) <- paste("age", ages) round(slopes, 3) # --- Answers --------------------------------------------------------------- cat("\nQ1. Support coefficient, uncentred (slope at age 0) vs centred (slope at mean age):\n") print(round(c(uncentred = unname(coef(model_raw)["support"]), centred = unname(coef(model_cen)["support_c"])), 4)) cat("\nQ2. Interaction coefficient and t value in both models (identical):\n") print(round(rbind(uncentred = summary(model_raw)$coefficients["support:age", 1:3], centred = summary(model_cen)$coefficients["support_c:age_c", 1:3]), 4)) cat("\nQ3. Simple slopes at ages 30, 50 and 70:\n") print(round(slopes, 3)) # ============================================================================= # SECTION 3, ACTIVITY 1: Mediation, one path at a time (Activity 8.4 on the lesson page) # ============================================================================= github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave 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) med_data <- na.omit(med_data) # the same people in every model nrow(med_data) # Path c, the total effect: support -> depression model_total <- lm(depression ~ support, data = med_data) round(cbind(estimate = coef(model_total), confint(model_total)), 3) # with 95% CIs # Path a: support -> loneliness model_a <- lm(loneliness ~ support, data = med_data) round(cbind(estimate = coef(model_a), confint(model_a)), 3) # with 95% CIs # Path b (loneliness -> depression) and path c' (the direct effect) model_b <- lm(depression ~ support + loneliness, data = med_data) round(cbind(estimate = coef(model_b), confint(model_b)), 3) # with 95% CIs a <- coef(model_a)["support"] b <- coef(model_b)["loneliness"] c_total <- coef(model_total)["support"] c_direct <- coef(model_b)["support"] a * b # indirect effect: product of coefficients c_total - c_direct # indirect effect: difference method # --- Answers --------------------------------------------------------------- cat("\nQ1. Paths with 95% CIs:\n") print(round(rbind(c_total = c(coef(model_total)["support"], confint(model_total)["support", ]), a = c(coef(model_a)["support"], confint(model_a)["support", ]), b = c(coef(model_b)["loneliness"], confint(model_b)["loneliness", ]), c_direct = c(coef(model_b)["support"], confint(model_b)["support", ])), 3)) cat("\nQ2. Product of coefficients and difference method:\n") print(round(c(a_times_b = unname(a * b), c_minus_cprime = unname(c_total - c_direct)), 4)) cat("\nQ3. Share of the total association that is indirect:\n") print(round(unname((a * b) / c_total), 3)) # ============================================================================= # SECTION 3, ACTIVITY 2: Bootstrap mediation and sensitivity analysis (Activity 8.5 on the lesson page) # ============================================================================= # install.packages("mediation") # run once, if not yet installed library(mediation) github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave 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) med_data <- na.omit(med_data) # the same people in every model model_a <- lm(loneliness ~ support, data = med_data) # path a model_b <- lm(depression ~ support + loneliness, data = med_data) # paths b and c-prime 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) plot(med_result, xlim = c(-0.7, 0)) # the three effects with 95% CIs # Sensitivity analysis: how strong would unmeasured confounding of # loneliness and depression have to be to remove the indirect effect? sens <- medsens(med_result, rho.by = 0.1, effect.type = "indirect", sims = 100) summary(sens) # --- Answers --------------------------------------------------------------- cat("\nQ1. ACME (indirect effect) with 95% bootstrap CI:\n") print(round(c(ACME = med_result$d0, lower = med_result$d0.ci[1], upper = med_result$d0.ci[2]), 3)) cat("\nQ2. Proportion mediated with 95% bootstrap CI:\n") print(round(c(prop = med_result$n0, lower = med_result$n0.ci[1], upper = med_result$n0.ci[2]), 2)) cat("\nQ3. Residual correlation (rho) at which the ACME equals zero:\n") print(sens$err.cr.d) # Unmeasured confounding of loneliness and depressive symptoms that induced a # residual correlation of about 0.4 would remove the indirect association. # ============================================================================= # SECTION 4, ACTIVITY 1: The simple mediation model in lavaan (Activity 8.6 on the lesson page) # ============================================================================= # install.packages("lavaan") # run once, if not yet installed library(lavaan) github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave path_data <- data.frame( support = data$PSYCH_zimet_multidimensional_social_support_scale_score, loneliness = data$LONELY_ucla_loneliness_scale_score, depression = data$WELLNESS_phq_score, anxiety = data$WELLNESS_gad_score) # GAD-2 anxiety symptoms, 0 to 6 path_data <- na.omit(path_data) # people with all four variables nrow(path_data) model_med <- ' loneliness ~ a * support depression ~ b * loneliness + c * support indirect := a * b total := c + a * b ' fit_med <- sem(model_med, data = path_data) summary(fit_med) # --- Answers --------------------------------------------------------------- cat("\nQ1 and Q2. Labelled paths and defined parameters:\n") pe <- parameterEstimates(fit_med) print(pe[pe$label != "", c("lhs", "op", "rhs", "label", "est", "se")], digits = 3) cat("\nQ3. Degrees of freedom (0 = saturated model):\n") print(fitMeasures(fit_med, c("chisq", "df"))) # ============================================================================= # SECTION 4, ACTIVITY 2: Two outcomes, model comparison and the chosen model (Activity 8.7 on the lesson page) # ============================================================================= # Model 1, full mediation: no direct paths from support model_full <- ' loneliness ~ support depression ~ loneliness anxiety ~ loneliness depression ~~ anxiety ' fit_full <- sem(model_full, data = path_data) fitMeasures(fit_full, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr")) # Model 2, partial mediation: adds a direct path to each outcome model_partial <- ' loneliness ~ a * support depression ~ b1 * loneliness + c1 * support anxiety ~ b2 * loneliness + c2 * support depression ~~ anxiety indirect_dep := a * b1 indirect_anx := a * b2 ' fit_partial <- sem(model_partial, data = path_data) fitMeasures(fit_partial, c("chisq", "df", "cfi", "rmsea", "srmr")) anova(fit_full, fit_partial) # chi-square difference test, AIC and BIC set.seed(2021) fit_boot <- sem(model_partial, data = path_data, se = "bootstrap", bootstrap = 1000) est <- parameterEstimates(fit_boot, boot.ci.type = "perc") est[est$op %in% c("~", ":="), c("lhs", "op", "rhs", "est", "ci.lower", "ci.upper")] std <- standardizedSolution(fit_boot) std[std$lhs != std$rhs, c("lhs", "op", "rhs", "est.std")] # paths, covariance, indirect lavInspect(fit_boot, "rsquare") # variance explained in each outcome # --- Answers --------------------------------------------------------------- cat("\nQ1. Fit of Model 1 (no direct paths):\n") print(round(fitMeasures(fit_full, c("chisq", "df", "cfi", "tli", "rmsea", "srmr")), 3)) cat("\nQ2. AIC of each model (lower is better):\n") print(round(c(model_1 = AIC(fit_full), model_2 = AIC(fit_partial)), 1)) cat("\nQ3. Indirect effects with 95% bootstrap CIs, and R-squared:\n") print(est[est$op == ":=", c("lhs", "est", "ci.lower", "ci.upper")], digits = 3) print(round(lavInspect(fit_boot, "rsquare"), 3)) # ============================================================================= # SECTION 4, R EXAMPLE: a structural equation model with a latent loneliness factor (Worked code 8.1 on the lesson page) # ============================================================================= library(lavaan) github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) # creates a data frame called data data <- data[data$SURVEY_collection_year == 2021, ] # keep the 2021 wave sem_data <- data.frame( companionship = data$LONELY_ucla_loneliness_scale_companionship_num, left_out = data$LONELY_ucla_loneliness_scale_left_out_num, isolated = data$LONELY_ucla_loneliness_scale_isolated_num, support = data$PSYCH_zimet_multidimensional_social_support_scale_score, depression = data$WELLNESS_phq_score) sem_data <- na.omit(sem_data) nrow(sem_data) model_sem <- ' # measurement model: a latent loneliness factor with three items Lonely =~ companionship + left_out + isolated # structural model: the paths between support, Lonely and depression Lonely ~ a * support depression ~ b * Lonely + c * support indirect := a * b ' fit_sem <- sem(model_sem, data = sem_data) fitMeasures(fit_sem, c("chisq", "df", "cfi", "rmsea", "srmr")) std <- standardizedSolution(fit_sem) std[std$op %in% c("=~", "~", ":="), c("lhs", "op", "rhs", "est.std", "pvalue")]