# ========================================================
# LATENT CLASS ANALYSIS IN R
# Canadian Social Connection Survey (CSCS), 2021 wave
# ========================================================
# Question: when people feel lonely, are there distinct
# groups (latent classes) that differ in what they miss?
# Survey question: "When you have felt lonely, what did
# you miss most? (Check all that apply)"
# ========================================================


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

# install.packages("poLCA")   # run once
library(poLCA)   # poLCA() fits latent class models

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: RECODE THE SEVEN CHECKBOX ITEMS
# ========================================================
# poLCA() needs each item coded 1, 2, 3, ... (no zeros), so:
#   1 = not selected, 2 = selected

# --- 2.1 Look at one item first
table(data$LONELY_miss_fun, useNA = "ifany")

# --- 2.2 Recode all seven items
miss <- data.frame(loneliness = data$LONELY_ucla_loneliness_scale_score,
                   age        = data$DEMO_age)

miss$fun <- NA
miss$fun[data$LONELY_miss_fun == "Not selected"] <- 1
miss$fun[data$LONELY_miss_fun == "Fun and laughter / Leisure"] <- 2

miss$conversation <- NA
miss$conversation[data$LONELY_miss_meaningful_conversation == "Not selected"] <- 1
miss$conversation[data$LONELY_miss_meaningful_conversation == "Meaningful conversation"] <- 2

miss$touch <- NA
miss$touch[data$LONELY_miss_touch == "Not selected"] <- 1
miss$touch[data$LONELY_miss_touch == "Physical touch/hug, affection"] <- 2

miss$hanging_out <- NA
miss$hanging_out[data$LONELY_miss_hanging_out == "Not selected"] <- 1
miss$hanging_out[data$LONELY_miss_hanging_out == "Being with other people/ hanging out"] <- 2

miss$understand <- NA
miss$understand[data$LONELY_miss_understand == "Not selected"] <- 1
miss$understand[data$LONELY_miss_understand == "Someone to understand me"] <- 2

miss$doing <- NA
miss$doing[data$LONELY_miss_doing == "Not selected"] <- 1
miss$doing[data$LONELY_miss_doing == "Doing something with other people"] <- 2

miss$mattering <- NA
miss$mattering[data$LONELY_miss_mattering == "Not selected"] <- 1
miss$mattering[data$LONELY_miss_mattering ==
                 "Mattering to someone/being able to help someone"] <- 2

# Check one recode
table(data$LONELY_miss_fun, miss$fun, useNA = "ifany")

# --- 2.3 Keep the people with no missing values
miss <- na.omit(miss)
nrow(miss)


# ========================================================
# SECTION 3: LOOK AT THE ITEMS
# ========================================================

# Percent of people who selected each item
round(100 * colMeans(miss[, 3:9] == 2), 1)

# How many items did each person select?
table(rowSums(miss[, 3:9] == 2))


# ========================================================
# SECTION 4: FIT MODELS WITH 1 TO 5 CLASSES
# ========================================================

# --- 4.1 The formula: the items go inside cbind(), and "~ 1"
# means no covariates
lca_formula <- cbind(fun, conversation, touch, hanging_out,
                     understand, doing, mattering) ~ 1

# --- 4.2 Fit each model. nrep = 10 starts each model from 10
# different random starting points and keeps the best one.
set.seed(2021)
lca_1 <- poLCA(lca_formula, data = miss, nclass = 1, verbose = FALSE)
lca_2 <- poLCA(lca_formula, data = miss, nclass = 2, nrep = 10, verbose = FALSE)
lca_3 <- poLCA(lca_formula, data = miss, nclass = 3, nrep = 10, verbose = FALSE)
lca_4 <- poLCA(lca_formula, data = miss, nclass = 4, nrep = 10, verbose = FALSE)
lca_5 <- poLCA(lca_formula, data = miss, nclass = 5, nrep = 10, verbose = FALSE)

# --- 4.3 Compare the models: lower AIC and BIC are better
fit_table <- data.frame(
  classes = 1:5,
  AIC = c(lca_1$aic, lca_2$aic, lca_3$aic, lca_4$aic, lca_5$aic),
  BIC = c(lca_1$bic, lca_2$bic, lca_3$bic, lca_4$bic, lca_5$bic))
round(fit_table, 1)

plot(fit_table$classes, fit_table$BIC, type = "b", pch = 16,
     main = "BIC for models with 1 to 5 classes",
     xlab = "Number of classes",
     ylab = "BIC (lower is better)")

# --- 4.4 What do the extra classes add? Class sizes in the
# 4- and 5-class models, and two items from the 4-class model
round(lca_4$P, 3)
round(lca_5$P, 3)
round(lca_4$probs$touch[, 2], 2)   # probability of selecting "touch"
round(lca_4$probs$fun[, 2], 2)     # probability of selecting "fun"


# ========================================================
# SECTION 5: THE THREE-CLASS MODEL
# ========================================================

# --- 5.1 The full output: item probabilities in each class
lca_3

# --- 5.2 Probability of SELECTING each item (the second
# column, code 2), one column per class. The row names are
# short labels for the plot.
yes_prob <- rbind(Fun        = lca_3$probs$fun[, 2],
                  Talk       = lca_3$probs$conversation[, 2],
                  Touch      = lca_3$probs$touch[, 2],
                  Company    = lca_3$probs$hanging_out[, 2],
                  Understood = lca_3$probs$understand[, 2],
                  Doing      = lca_3$probs$doing[, 2],
                  Mattering  = lca_3$probs$mattering[, 2])
round(yes_prob, 2)

barplot(t(yes_prob), beside = TRUE, ylim = c(0, 1.2),
        main = "What each class misses",
        ylab = "Probability of selecting the item",
        col = c("grey20", "grey55", "grey85"),
        legend.text = c("Class 1", "Class 2", "Class 3"),
        args.legend = list(x = "top", horiz = TRUE, bty = "n"),
        cex.names = 0.7)


# ========================================================
# SECTION 6: ASSIGN PEOPLE TO CLASSES
# ========================================================

# --- 6.1 Each person's probability of belonging to each class
round(head(lca_3$posterior), 3)

# --- 6.2 Assign each person to their most likely class
miss$class <- lca_3$predclass
table(miss$class)


# ========================================================
# SECTION 7: DESCRIBE THE CLASSES
# ========================================================

# Average loneliness score and age in each class
tapply(miss$loneliness, miss$class, mean)
tapply(miss$age, miss$class, mean)

# Average number of boxes ticked in each class
tapply(rowSums(miss[, 3:9] == 2), miss$class, mean)

boxplot(loneliness ~ class, data = miss,
        main = "Loneliness score by latent class",
        xlab = "Latent class",
        ylab = "UCLA loneliness score (3 to 9)",
        col = "grey80")