# =============================================================================
# HSCI 410 Lesson 7: Measurement and Psychometrics -- ANSWER KEY
# Author: Kiffer G. Card, PhD - Faculty of Health Sciences, SFU
# -----------------------------------------------------------------------------
# Page:  modules/HSCI_410_Lesson_7_Measurement_and_Psychometrics.html
# Data:  Canadian Social Connection Survey (CSCS), 2021 wave, loaded from GitHub
#        (an internet connection is needed). The six De Jong Gierveld items are
#        prepared exactly as in the walkthrough
#        r-walkthroughs/HSCI_410_Factor_Analysis_and_Scale_Scoring.html,
#        so every number matches it.
#
# This script reproduces every code block on the page, in page order, and then
# answers each numbered activity question with a printed result.
#
# Packages: install.packages(c("psych", "GPArotation", "lavaan"))
#           install.packages("mirt")   # only for the optional last section
# =============================================================================

options(width = 80)


# =============================================================================
# SECTION 1, ACTIVITY 1: Load the CSCS data and prepare the six items (Activity 7.1 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))     # loads a data frame called data
dim(data)                 # rows (survey responses) and columns

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

# 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)

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

# 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

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

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

# Reverse-code the three positively worded items: for these
# items "No" (1) is the lonely answer. Subtracting from 4 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
table(dj$rely_n, dj$rely_r)      # check: 1 -> 3, 2 -> 2, 3 -> 1

# 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)

# --- Answers ---------------------------------------------------------------
# Q1. The rely, trust and close items are positively worded, so "No" (code 1) is
#     the lonely answer; without reversal they would count in the opposite
#     direction to the emotional items. The table shows 1 -> 3 (546 people),
#     2 -> 2 (1,526) and 3 -> 1 (1,409), with no other cells filled.
cat("\nQ2. People in the 2021 wave, and people who answered all six items:\n")
print(c(wave_2021 = nrow(data), complete_six_items = nrow(items)))
# "Presented but no response" was not given a number, so it stayed NA and
# na.omit() removed those people.
# Q3. Each item is ordinal (ordered answers, unknown spacing). Adding or averaging
#     the codes treats them as equally spaced (approximately interval).

# =============================================================================
# SECTION 1, ACTIVITY 2: Score the scale with its published rule (Activity 7.2 on the lesson page)
# =============================================================================
# The published scoring rule: an item scores 1 for "More or
# less" or the lonely answer (a code of 2 or 3), and 0 otherwise.
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)

# 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")

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

# 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")

# --- Answers ---------------------------------------------------------------
# Q1. For "There are enough people I feel close to", "No" and "More or less"
#     score 1 and "Yes" scores 0; after reversal these are close_r = 3 and 2,
#     which is why the code uses dj$close_r >= 2.
cat("\nQ2. People classified as lonely (total 2 to 6) and the percentage:\n")
n_lonely <- sum(dj$total_score >= 2, na.rm = TRUE)
print(c(lonely = n_lonely, percent = round(100 * n_lonely / sum(!is.na(dj$total_score)), 1)))
# Q3. Agreement with the supplied score checks the recoding, the reversal and
#     the scoring rule at once; all 3,415 totals match.

# =============================================================================
# SECTION 2: Cronbach's alpha, omega and a missed reversal (Activity 7.3 on the lesson page)
# =============================================================================
# install.packages("psych")   # run once, if not yet installed
library(psych)    # alpha() and omega()
emotional_items <- items[, c("emptiness", "miss", "rejected")]
social_items    <- items[, c("rely", "trust", "close")]
alpha(social_items)

round(alpha(emotional_items)$total, 2)               # overall alpha only
round(alpha(emotional_items)$alpha.drop[, 1:2], 2)   # alpha if each item is dropped
round(alpha(emotional_items)$item.stats[, "r.drop", drop = FALSE], 2)   # item-total

# Alpha for all six items, correctly reverse-coded
round(alpha(items)$total[, 1:2], 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)
round(alpha(na.omit(not_reversed))$total[, 1:2], 2)

# McDonald's omega from a one-factor model of each subscale
f_social <- fa(social_items, nfactors = 1)   # one-factor model
lambda <- f_social$loadings[, 1]             # the three loadings
round(lambda, 2)
sum(lambda)^2 / (sum(lambda)^2 + sum(f_social$uniquenesses))   # omega

f_emo <- fa(emotional_items, nfactors = 1)
lambda <- f_emo$loadings[, 1]
round(lambda, 2)
sum(lambda)^2 / (sum(lambda)^2 + sum(f_emo$uniquenesses))

# An index: the Steptoe social isolation index (CSCS supplies 0/1 codes,
# where 1 = isolated on that component)
iso <- data.frame(unmarried = data$LONELY_steptoe_isolation_index_unmarried_num,
                  clubs     = data$LONELY_steptoe_isolation_index_clubs_num,
                  friends   = data$LONELY_steptoe_isolation_index_friends_num,
                  children  = data$LONELY_steptoe_isolation_index_kids_num,
                  family    = data$LONELY_steptoe_isolation_index_other_fam_num)
iso <- na.omit(iso)
round(cor(iso), 2)
round(alpha(iso)$total[, 1:2], 2)   # shown only to explain why it is not reported

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Alpha (with Feldt 95% CI) for the social subscale, and alpha for the emotional subscale:\n")
f <- alpha(social_items)$feldt
print(round(c(lower = f$lower.ci[[1]], alpha = f$alpha[[1]], upper = f$upper.ci[[1]]), 2))
print(round(alpha(emotional_items)$total$raw_alpha, 2))
cat("\nQ2. Emotional subscale: alpha if each item is dropped, and corrected item-total r:\n")
a_emo <- alpha(emotional_items)
print(round(data.frame(alpha_if_dropped = a_emo$alpha.drop$raw_alpha,
                       r_drop = a_emo$item.stats$r.drop,
                       row.names = rownames(a_emo$alpha.drop)), 2))
# The weakest item is miss; keep the published subscale as the main measure and
# run a sensitivity analysis with the two-item version.
cat("\nAlpha for the two-item emotional subscale (sensitivity analysis):\n")
print(round(alpha(emotional_items[, c("emptiness", "rejected")])$total$raw_alpha, 2))
# Q3. Omega equals alpha for the social items because their loadings are similar
#     (0.73, 0.73, 0.67); for the emotional items the loadings differ (0.83,
#     0.20, 0.58), so alpha (0.51) underestimates reliability and omega is 0.57.

# =============================================================================
# SECTION 3: Convergent, discriminant and known-groups evidence (Activity 7.4 on the lesson page)
# =============================================================================
# Other measures in the CSCS file (already scored by the survey team)
val <- data.frame(total     = dj$total_score,
                  social    = dj$social_score,
                  emotional = dj$emotional_score,
                  ucla      = data$LONELY_ucla_loneliness_scale_score,
                  support   = data$PSYCH_zimet_multidimensional_social_support_scale_score,
                  anxiety   = data$WELLNESS_gad_score)
val <- na.omit(val)     # people with every score
nrow(val)
round(cor(val), 2)      # correlations between every pair of measures

# Known groups: people with no close friends should score higher
friends <- data$CONNECTION_social_num_close_friends_grouped
friends[friends == "Presented but no response"] <- NA
friends <- factor(friends, levels = c("None", "1–2", "3–4", "5 or more"))
round(tapply(dj$social_score, friends, mean, na.rm = TRUE), 2)
round(tapply(dj$emotional_score, friends, mean, na.rm = TRUE), 2)

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Correlation of each subscale with perceived social support:\n")
print(round(cor(val$social, val$support), 2))
print(round(cor(val$emotional, val$support), 2))
# The social subscale's association is stronger, as hypothesized.
# Q2. Convergent: UCLA 0.45 and social support -0.49; discriminant: anxiety 0.33.
cat("\nQ3. Drop in mean subscale score from no close friends to five or more:\n")
soc <- tapply(dj$social_score, friends, mean, na.rm = TRUE)
emo <- tapply(dj$emotional_score, friends, mean, na.rm = TRUE)
print(round(c(social = soc[["None"]] - soc[["5 or more"]], emotional = emo[["None"]] - emo[["5 or more"]]), 2))

# =============================================================================
# SECTION 4, ACTIVITY 1: Exploratory factor analysis (Activity 7.5 on the lesson page)
# =============================================================================
# install.packages(c("psych", "GPArotation"))   # run once
library(psych)        # KMO(), cortest.bartlett(), fa.parallel(), fa()
library(GPArotation)  # the rotation methods that fa() uses
round(cor(items), 2)  # correlations between every pair of items

# Explore in one half, test in the other half
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)

KMO(efa_half)    # overall MSA above 0.6 is usually adequate
cortest.bartlett(cor(efa_half), n = nrow(efa_half))

fa.parallel(efa_half, fa = "fa")   # how many factors?

efa_2 <- fa(efa_half, nfactors = 2, rotate = "oblimin")
print(efa_2$loadings, cutoff = 0.3)  # hide loadings below 0.3
round(efa_2$Phi, 2)                  # correlation between the factors

# --- Answers ---------------------------------------------------------------
# Q1. Overall KMO = 0.68 (adequate); parallel analysis suggests 2 factors. It
#     keeps the factors whose eigenvalues exceed those from random data.
cat("\nQ2. Loadings (all values shown, rounded):\n")
print(round(unclass(efa_2$loadings), 2))
# Social items load on MR1, emptiness and rejected on MR2; miss is below 0.3 on both.
# Q3. Oblimin lets the factors correlate; the factor correlation is 0.29.

# =============================================================================
# SECTION 4, ACTIVITY 2: Confirmatory factor analysis in lavaan (Activity 7.6 on the lesson page)
# =============================================================================
# install.packages("lavaan")   # run once
library(lavaan)    # cfa() for confirmatory factor analysis
# "=~" 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)

fitMeasures(cfa_2f, c("cfi", "tli", "rmsea", "srmr"))

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)    # chi-square difference test

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Fit indices, two-factor model:\n")
print(round(fitMeasures(cfa_2f, c("cfi", "tli", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr")), 3))
cat("\nQ2. Standardized loadings and the share of variance explained (squared):\n")
std <- standardizedSolution(cfa_2f)
std <- std[std$op == "=~", c("lhs", "rhs", "est.std")]
std$explained <- std$est.std^2
print(std, digits = 2)
# Q3. The chi-square difference (386.95 on 1 df, p < 2.2e-16) and the lower AIC
#     favour two factors; wording direction may contribute to the separation.

# =============================================================================
# OPTIONAL (not assessed): a graded response model with mirt
# =============================================================================
if (requireNamespace("mirt", quietly = TRUE)) {
  # OPTIONAL: install.packages("mirt")   # run once
  library(mirt)
  # A graded response model for the three social loneliness items
  grm <- mirt(social_items, model = 1, itemtype = "graded", verbose = FALSE)
  print(round(coef(grm, IRTpars = TRUE, simplify = TRUE)$items, 2))

  print(plot(grm, type = "trace"))    # category curves for each item
  print(plot(grm, type = "info"))     # test information curve
} else {
  message("Install the mirt package to run this optional section.")
}