# ========================================================
# FACTOR ANALYSIS, RELIABILITY AND SCALE SCORES IN R
# Canadian Social Connection Survey (CSCS), 2021 wave
# ========================================================
# Question: do the six items of the De Jong Gierveld
# Loneliness Scale measure the two kinds of loneliness the
# scale was designed to measure (emotional and social), and
# how are the answers turned into a score?
# ========================================================


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

# install.packages(c("psych", "GPArotation", "lavaan"))   # run once
library(psych)        # fa.parallel(), fa(), KMO() and alpha()
library(GPArotation)  # the rotation methods that fa() uses
library(lavaan)       # cfa() for confirmatory factor analysis

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

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


# ========================================================
# SECTION 2: THE SIX ITEMS
# ========================================================
# Each statement is answered "Yes", "More or less" or "No".
# Emotional loneliness (negatively worded: "Yes" = lonely)
#   emptiness: I experience a general sense of emptiness
#   miss:      I miss having people around me
#   rejected:  I often feel rejected
# Social loneliness (positively worded: "No" = lonely)
#   rely:      There are plenty of people I can rely on
#              when I have problems
#   trust:     There are many people I can trust completely
#   close:     There are enough people I feel close to

# A small data frame with short names for the six items
dj <- data.frame(
  emptiness = data$LONELY_dejong_emotional_social_loneliness_scale_emptiness,
  miss      = data$LONELY_dejong_emotional_social_loneliness_scale_miss,
  rejected  = data$LONELY_dejong_emotional_social_loneliness_scale_rejected,
  rely      = data$LONELY_dejong_emotional_social_loneliness_scale_rely,
  trust     = data$LONELY_dejong_emotional_social_loneliness_scale_trust,
  close     = data$LONELY_dejong_emotional_social_loneliness_scale_close)
head(dj)

# The answers are stored as text categories
table(dj$emptiness, useNA = "ifany")
table(dj$rely, useNA = "ifany")


# ========================================================
# SECTION 3: RECODE THE ANSWERS
# ========================================================

# --- 3.1 Text to numbers: No = 1, More or less = 2, Yes = 3
# Anything else ("Presented but no response" or NA) stays NA
dj$emptiness_n <- NA
dj$emptiness_n[dj$emptiness == "No"] <- 1
dj$emptiness_n[dj$emptiness == "More or less"] <- 2
dj$emptiness_n[dj$emptiness == "Yes"] <- 3

# Check: each answer should map to one number
table(dj$emptiness, dj$emptiness_n, useNA = "ifany")

dj$miss_n <- NA
dj$miss_n[dj$miss == "No"] <- 1
dj$miss_n[dj$miss == "More or less"] <- 2
dj$miss_n[dj$miss == "Yes"] <- 3

dj$rejected_n <- NA
dj$rejected_n[dj$rejected == "No"] <- 1
dj$rejected_n[dj$rejected == "More or less"] <- 2
dj$rejected_n[dj$rejected == "Yes"] <- 3

dj$rely_n <- NA
dj$rely_n[dj$rely == "No"] <- 1
dj$rely_n[dj$rely == "More or less"] <- 2
dj$rely_n[dj$rely == "Yes"] <- 3

dj$trust_n <- NA
dj$trust_n[dj$trust == "No"] <- 1
dj$trust_n[dj$trust == "More or less"] <- 2
dj$trust_n[dj$trust == "Yes"] <- 3

dj$close_n <- NA
dj$close_n[dj$close == "No"] <- 1
dj$close_n[dj$close == "More or less"] <- 2
dj$close_n[dj$close == "Yes"] <- 3

# --- 3.2 Reverse-code the three positively worded items
# For these items "No" (1) is the lonely answer. Subtracting
# from 4 (the highest code plus the lowest code) flips the
# scale: 1 becomes 3, 2 stays 2, and 3 becomes 1.
dj$rely_r  <- 4 - dj$rely_n
dj$trust_r <- 4 - dj$trust_n
dj$close_r <- 4 - dj$close_n

# Check: 1 -> 3, 2 -> 2, 3 -> 1
table(dj$rely_n, dj$rely_r)

# --- 3.3 One data frame of six items, higher = lonelier
items <- data.frame(emptiness = dj$emptiness_n,
                    miss      = dj$miss_n,
                    rejected  = dj$rejected_n,
                    rely      = dj$rely_r,
                    trust     = dj$trust_r,
                    close     = dj$close_r)
items <- na.omit(items)   # people who answered all six items
nrow(items)


# ========================================================
# SECTION 4: LOOK AT HOW THE ITEMS RELATE
# ========================================================

# Correlations between every pair of items
round(cor(items), 2)


# ========================================================
# SECTION 5: SPLIT THE SAMPLE IN HALF
# ========================================================
# Explore the structure in one half, then test it in the
# other half, so the test uses people the exploration did
# not see.
set.seed(2021)    # makes the random split the same every time
efa_rows <- sample(nrow(items), size = round(nrow(items) / 2))
efa_half <- items[efa_rows, ]     # half 1: exploratory
cfa_half <- items[-efa_rows, ]    # half 2: confirmatory
nrow(efa_half)
nrow(cfa_half)


# ========================================================
# SECTION 6: EXPLORATORY FACTOR ANALYSIS (EFA)
# ========================================================

# --- 6.1 Are the items suitable for factor analysis?
# KMO above 0.6 is usually considered adequate
KMO(efa_half)

# --- 6.2 How many factors? Parallel analysis compares the
# data with random data of the same size
fa.parallel(efa_half, fa = "fa")

# --- 6.3 Fit a two-factor model
# "oblimin" rotation lets the two factors be correlated
efa_2 <- fa(efa_half, nfactors = 2, rotate = "oblimin")

# Loadings below 0.3 are hidden to make the pattern easy to see
print(efa_2$loadings, cutoff = 0.3)

# The correlation between the two factors
round(efa_2$Phi, 2)


# ========================================================
# SECTION 7: CONFIRMATORY FACTOR ANALYSIS (CFA)
# ========================================================

# --- 7.1 Write the two-factor model in lavaan syntax
# "=~" means "is measured by"
model_2f <- '
  emotional =~ emptiness + miss + rejected
  social    =~ rely + trust + close
'
cfa_2f <- cfa(model_2f, data = cfa_half)
summary(cfa_2f, fit.measures = TRUE, standardized = TRUE)

# --- 7.2 The four fit indices most often reported
fitMeasures(cfa_2f, c("cfi", "tli", "rmsea", "srmr"))

# --- 7.3 Compare with a one-factor model
model_1f <- '
  loneliness =~ emptiness + miss + rejected + rely + trust + close
'
cfa_1f <- cfa(model_1f, data = cfa_half)
fitMeasures(cfa_1f, c("cfi", "tli", "rmsea", "srmr"))
anova(cfa_1f, cfa_2f)


# ========================================================
# SECTION 8: RELIABILITY: CRONBACH'S ALPHA
# ========================================================

# --- 8.1 Alpha for each subscale (all 3,415 people)
emotional_items <- items[, c("emptiness", "miss", "rejected")]
social_items    <- items[, c("rely", "trust", "close")]
alpha(emotional_items)
alpha(social_items)

# --- 8.2 What happens if we forget to reverse-code?
not_reversed <- data.frame(emptiness = dj$emptiness_n,
                           miss      = dj$miss_n,
                           rejected  = dj$rejected_n,
                           rely      = dj$rely_n,
                           trust     = dj$trust_n,
                           close     = dj$close_n)
alpha(na.omit(not_reversed))$total


# ========================================================
# SECTION 9: SCORE THE SCALE
# ========================================================

# --- 9.1 The published scoring rule: an item scores 1 if the
# answer is "More or less" or the lonely answer (a code of 2
# or 3 on our recoded items), and 0 otherwise.
# A comparison gives TRUE or FALSE; as.numeric() turns
# TRUE into 1 and FALSE into 0.
dj$emptiness_01 <- as.numeric(dj$emptiness_n >= 2)
dj$miss_01      <- as.numeric(dj$miss_n >= 2)
dj$rejected_01  <- as.numeric(dj$rejected_n >= 2)
dj$rely_01      <- as.numeric(dj$rely_r >= 2)
dj$trust_01     <- as.numeric(dj$trust_r >= 2)
dj$close_01     <- as.numeric(dj$close_r >= 2)

# --- 9.2 Add the items: two subscales (0 to 3) and a total (0 to 6)
dj$emotional_score <- dj$emptiness_01 + dj$miss_01 + dj$rejected_01
dj$social_score    <- dj$rely_01 + dj$trust_01 + dj$close_01
dj$total_score     <- dj$emotional_score + dj$social_score
table(dj$total_score, useNA = "ifany")

# --- 9.3 Check against the score supplied with the CSCS data
table(dj$total_score == data$LONELY_dejong_emotional_social_loneliness_scale_score,
      useNA = "ifany")

# --- 9.4 Plot the total score
barplot(table(dj$total_score),
        main = "De Jong Gierveld Loneliness Scale scores",
        xlab = "Total score (0 = not lonely, 6 = most lonely)",
        ylab = "Number of participants",
        col = "grey80")