# ========================================================
# PATH ANALYSIS WITH STRUCTURAL EQUATION MODELLING IN R
# Canadian Social Connection Survey (CSCS), 2021 wave
# ========================================================
# Question: does loneliness account for the association between
# social support and two mental health outcomes, depressive
# symptoms and anxiety symptoms?
# ========================================================


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

# install.packages("lavaan")   # run once
library(lavaan)   # sem() fits path models and other SEMs

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:    social support, 1 to 7 (higher = more support)
# loneliness: UCLA 3-item Loneliness Scale, 3 to 9
# depression: PHQ-2 depressive symptoms, 0 to 6
# anxiety:    GAD-2 anxiety symptoms, 0 to 6
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)
path_data <- na.omit(path_data)   # people with all four variables
nrow(path_data)
summary(path_data)

# Correlations between the four variables
round(cor(path_data), 2)


# ========================================================
# SECTION 3: A FIRST PATH MODEL: SIMPLE MEDIATION
# ========================================================
# lavaan syntax:
#   y ~ x     y is regressed on x (an arrow from x to y)
#   a * x     gives the path from x the label "a"
#   :=        defines a new quantity from labelled paths
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)


# ========================================================
# SECTION 4: A PATH MODEL WITH TWO OUTCOMES
# ========================================================

# --- 4.1 Model 1, full mediation: support affects both
# outcomes only through loneliness (no direct paths)
#   y1 ~~ y2   lets the two outcomes stay correlated
model_full <- '
  loneliness ~ support
  depression ~ loneliness
  anxiety    ~ loneliness
  depression ~~ anxiety
'
fit_full <- sem(model_full, data = path_data)
summary(fit_full, fit.measures = TRUE, standardized = TRUE)

# The fit indices most often reported
fitMeasures(fit_full, c("chisq", "df", "pvalue", "cfi", "tli",
                        "rmsea", "srmr"))

# --- 4.2 Model 2, partial mediation: adds direct paths from
# support to each outcome, and labels every path
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"))

# --- 4.3 Compare the two models
anova(fit_full, fit_partial)


# ========================================================
# SECTION 5: RESULTS FROM THE CHOSEN MODEL
# ========================================================

# --- 5.1 Bootstrap confidence intervals for the indirect
# effects (lavaan refits the model 1,000 times)
set.seed(2021)
fit_boot <- sem(model_partial, data = path_data,
                se = "bootstrap", bootstrap = 1000)
parameterEstimates(fit_boot, boot.ci.type = "perc")

# --- 5.2 Standardized paths and variance explained (R-squared)
standardizedSolution(fit_boot)
lavInspect(fit_boot, "rsquare")