# =============================================================================
# HSCI 410 Lesson 6: Exploratory Data Analysis and Visualization -- ANSWER KEY
# Author: Kiffer G. Card, PhD - Faculty of Health Sciences, SFU
# -----------------------------------------------------------------------------
# Page:  modules/HSCI_410_Lesson_6_Exploratory_Data_Analysis_and_Visualization.html
# Data:  Canadian Social Connection Survey (CSCS), 2021 wave, loaded from the
#        course GitHub repository (an internet connection is needed). Section 1
#        uses the anscombe data frame that comes with R.
#
# This script reproduces every code block on the page, in page order, and then
# answers each numbered activity question with a printed result. It runs from
# top to bottom in one R session. Plots appear in the Plots pane, and three
# files are written to the working directory: briefing_base_r.png,
# loneliness_age_gender.png and briefing_figure.png.
#
# Packages: ggplot2, ggExtra and patchwork. Install them once with
#   install.packages(c("ggplot2", "ggExtra", "patchwork"))
# Data preparation uses base R only (bracket recoding and na.omit()).
# =============================================================================

options(width = 70)


# =============================================================================
# SECTION 1 ACTIVITY: Anscombe's quartet in R (Activity 6.1 on the lesson page)
# =============================================================================
anscombe                          # four small datasets that come with R
round(colMeans(anscombe), 2)      # the mean of each column
round(sapply(anscombe, sd), 2)    # the standard deviation of each column

round(cor(anscombe$x1, anscombe$y1), 3)   # correlation in set 1
round(cor(anscombe$x2, anscombe$y2), 3)   # set 2
round(cor(anscombe$x3, anscombe$y3), 3)   # set 3
round(cor(anscombe$x4, anscombe$y4), 3)   # set 4
coef(lm(y1 ~ x1, data = anscombe))        # intercept and slope of the line, set 1
coef(lm(y4 ~ x4, data = anscombe))        # the same for set 4

par(mfrow = c(2, 2))    # four plots in one window: two rows, two columns
plot(anscombe$x1, anscombe$y1, main = "Set 1", xlab = "x", ylab = "y")
plot(anscombe$x2, anscombe$y2, main = "Set 2", xlab = "x", ylab = "y")
plot(anscombe$x3, anscombe$y3, main = "Set 3", xlab = "x", ylab = "y")
plot(anscombe$x4, anscombe$y4, main = "Set 4", xlab = "x", ylab = "y")
par(mfrow = c(1, 1))    # back to one plot per window

# --- Answers ---------------------------------------------------------------
# Q1. All four sets share: mean x = 9.0, mean y = 7.5, SD x = 3.32,
#     SD y = 2.03, r = 0.816 (0.817 in set 4) and the line y = 3.00 + 0.500x.
cat("\nQ1. Intercept and slope in each set:\n")
fits <- sapply(1:4, function(i) coef(lm(anscombe[[paste0("y", i)]] ~ anscombe[[paste0("x", i)]])))
colnames(fits) <- paste("Set", 1:4); rownames(fits) <- c("intercept", "slope")
print(round(fits, 3))
# Q2. Set 1: a linear cloud (the line fits). Set 2: a smooth curve. Set 3: a
#     near-perfect line with one outlier. Set 4: ten points at x = 8 and one
#     point at x = 19 that creates the whole relationship.
cat("\nQ2. Set 2 is exactly quadratic; a squared term fits it almost perfectly:\n")
print(round(summary(lm(y2 ~ x2 + I(x2^2), data = anscombe))$r.squared, 4))
# Q3. Identical summaries can hide a curve, an outlier or a relationship
#     created by one point, so the briefing shows the data before any model.

# =============================================================================
# SECTION 2, ACTIVITY 1: preparing the briefing file and one variable (Activity 6.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))
data <- data[data$SURVEY_collection_year == 2021, ]   # the 2021 wave: one row per person
nrow(data)                                             # how many people?

# Short names for the briefing variables
data$loneliness <- data$LONELY_ucla_loneliness_scale_score   # 3 to 9, higher = lonelier
data$support    <- data$PSYCH_zimet_multidimensional_social_support_scale_score   # 1 to 7
data$depression <- data$WELLNESS_phq_score                   # 0 to 6, higher = more symptoms
data$anxiety    <- data$WELLNESS_gad_score                   # 0 to 6, higher = more symptoms
data$life_sat   <- data$WELLNESS_life_satisfaction_num       # 1 to 10, higher = more satisfied
data$age        <- data$DEMO_age                             # years

# Gender: the answer "Presented but no response" becomes missing (NA)
data$gender <- as.character(data$DEMO_gender)
data$gender[data$gender == "Presented but no response"] <- NA
data$gender <- factor(data$gender, levels = c("Woman", "Man", "Non-binary"))

# Self-rated mental health, in order from Poor to Excellent
data$mental_health <- as.character(data$WELLNESS_self_rated_mental_health)
data$mental_health[data$mental_health == "Presented but no response"] <- NA
data$mental_health <- factor(data$mental_health,
  levels = c("Poor", "Fair", "Good", "Very good", "Excellent"))

# Four age groups
data$age_group[data$age < 30] <- "16 to 29"
data$age_group[data$age >= 30 & data$age < 45] <- "30 to 44"
data$age_group[data$age >= 45 & data$age < 65] <- "45 to 64"
data$age_group[data$age >= 65] <- "65 and older"
data$age_group <- factor(data$age_group,
  levels = c("16 to 29", "30 to 44", "45 to 64", "65 and older"))

# The briefing file: the briefing variables for people with no missing values
brief <- na.omit(data[, c("loneliness", "support", "depression", "anxiety",
                          "life_sat", "age", "age_group", "gender", "mental_health")])
nrow(brief)                                            # people in the briefing file

summary(brief$loneliness)              # a numeric summary first
hist(brief$loneliness,
     breaks = seq(2.5, 9.5, by = 1),       # one bar for each score from 3 to 9
     main = "Loneliness (UCLA 3-item scale)",
     xlab = "Loneliness score (3 to 9)", ylab = "Number of people")

table(brief$mental_health)             # the counts that the bars will show
barplot(table(brief$mental_health),
        main = "Self-rated mental health",
        xlab = "Rating", ylab = "Number of people")

summary(data$GEO_housing_household_size)   # people each respondent lives with
hist(data$GEO_housing_household_size,
     breaks = seq(-0.5, 20.5, by = 1),             # one bar for each value from 0 to 20
     main = "Household size",
     xlab = "Number of people the respondent lives with", ylab = "Number of people")
sum(data$GEO_housing_household_size == 20, na.rm = TRUE)   # how many at the top value?

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Counts and percentages for each loneliness score:\n")
print(table(brief$loneliness))
print(round(100 * prop.table(table(brief$loneliness)), 1))
# Most common score 6 (755), then 5 (689); 270 people (8.8%) at the maximum of 9,
# more than at 8 (179): a ceiling pile-up. Mean 5.57 > median 5.
cat("\nQ2. Percentage rating mental health as poor or fair:\n")
poor_fair <- sum(brief$mental_health %in% c("Poor", "Fair"))
print(c(people = poor_fair, percent = round(100 * poor_fair / nrow(brief), 1)))
cat("\nQ3. Recode the top-coded household value 20 to missing:\n")
data$household <- data$GEO_housing_household_size
data$household[data$household == 20] <- NA
print(summary(data$household))

# --- Checks behind the reading ----------------------------------------------
cat("\nHousehold size is the sum of eight 'live with' questions, capped at 20:\n")
live_with <- c("GEO_housing_live_with_partner", "GEO_housing_live_with_children",
               "GEO_housing_live_with_grandkids", "GEO_housing_live_with_parent",
               "GEO_housing_live_with_in_laws", "GEO_housing_live_with_siblings",
               "GEO_housing_live_with_roommate", "GEO_housing_live_with_other")
total <- rowSums(data[, live_with], na.rm = TRUE)
hh <- data$GEO_housing_household_size
print(table(capped_sum_matches = pmin(total, 20)[!is.na(hh)] == hh[!is.na(hh)]))
print(c(at_20 = sum(hh == 20, na.rm = TRUE),
        total_above_20 = sum(total > 20 & hh == 20, na.rm = TRUE),
        largest_total = max(total[!is.na(hh) & hh == 20]),
        living_alone = sum(hh == 0, na.rm = TRUE)))
cat("\nHours worked per week: heaping and implausible values:\n")
hours <- data$WORK_hours_per_week
print(c(answered = sum(!is.na(hours)), at_39 = sum(hours == 39, na.rm = TRUE),
        at_40 = sum(hours == 40, na.rm = TRUE), at_41 = sum(hours == 41, na.rm = TRUE),
        at_35 = sum(hours == 35, na.rm = TRUE), above_112 = sum(hours > 112, na.rm = TRUE),
        maximum = max(hours, na.rm = TRUE)))
hist(hours, breaks = seq(0, 140, by = 2), main = "Hours worked per week",
     xlab = "Hours", ylab = "Number of people")

# =============================================================================
# SECTION 2, ACTIVITY 2: two variables, layouts and saving a plot (Activity 6.3 on the lesson page)
# =============================================================================
boxplot(loneliness ~ age_group, data = brief,
        main = "Loneliness by age group",
        xlab = "Age group", ylab = "Loneliness score (3 to 9)")
tapply(brief$loneliness, brief$age_group, median)   # the thick line in each box

plot(brief$support, brief$loneliness,              # every person as one point
     xlab = "Social support (1 to 7)", ylab = "Loneliness score (3 to 9)")
set.seed(2021)                                       # the same jitter every time
plot(jitter(brief$support), jitter(brief$loneliness),   # nudge each point a little
     pch = 16, col = rgb(0, 0, 0, 0.15),             # small grey see-through dots
     xlab = "Social support (1 to 7)", ylab = "Loneliness score (3 to 9)")

pairs(brief[, c("loneliness", "support", "depression", "age")],
      pch = 16, col = rgb(0, 0, 0, 0.1))   # every pair of variables in one grid

png("briefing_base_r.png", width = 1600, height = 800, res = 200)   # open a file
par(mfrow = c(1, 2))                       # one row, two plots
hist(brief$loneliness, breaks = seq(2.5, 9.5, by = 1),
     main = "A. Loneliness", xlab = "Loneliness score", ylab = "Number of people")
boxplot(loneliness ~ age_group, data = brief,
        main = "B. Loneliness by age group", xlab = "Age group", ylab = "Loneliness score")
par(mfrow = c(1, 1))
dev.off()                                  # close the file, which saves it

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Boxplot statistics (whisker, quartile, median, quartile, whisker):\n")
bp <- boxplot(loneliness ~ age_group, data = brief, plot = FALSE)
stats <- bp$stats; colnames(stats) <- bp$names
print(stats)
print(c(points_beyond_whiskers_16_to_29 = sum(bp$group == 1)))
# Medians 5, 5, 6, 6; the 45 to 64 box is widest (4 to 8); 16 to 29 is narrow (5 to 6).
# Q2. Jitter and transparency show most people between support 4 and 6 and
#     loneliness 5 or 6, and a downward trend.
# Q3. png() sent the plots to briefing_base_r.png in the working directory:
cat("\nQ3. The working directory:\n")
print(getwd())

# =============================================================================
# SECTION 3, ACTIVITY 3: first plots in ggplot2 (Activity 6.4 on the lesson page)
# =============================================================================
# install.packages("ggplot2")   # run once, if not yet installed
library(ggplot2)
ggplot(brief, aes(x = loneliness)) +          # data and the aesthetic mapping
  geom_histogram(binwidth = 1,                # the geom: one bar per score
                 fill = "#0B7B6B", colour = "white") +
  labs(title = "Loneliness scores",           # labels
       x = "Loneliness score (3 to 9)", y = "Number of people")

ggplot(brief, aes(x = mental_health)) +
  geom_bar() +                                # geom_bar() counts the rows in each category
  labs(x = "Self-rated mental health", y = "Number of people")

means <- aggregate(loneliness ~ age_group, data = brief, FUN = mean)
means                                         # one row per age group
ggplot(means, aes(x = age_group, y = loneliness)) +
  geom_col() +                                # geom_col() draws the values it is given
  labs(x = "Age group", y = "Mean loneliness score")

set.seed(2021)                                # the same jitter every time
ggplot(brief, aes(x = age_group, y = loneliness)) +
  geom_boxplot(outlier.shape = NA) +          # the box summarizes each group
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.1) +   # every person as a faint dot
  labs(x = "Age group", y = "Loneliness score (3 to 9)")

# --- Answers ---------------------------------------------------------------
# Q1. Data: brief. Mapping: aes(x = loneliness). Geom: geom_histogram().
#     Labels: labs(). The fill and outline colours are set outside aes().
cat("\nQ2. Mean loneliness by age group, two decimals:\n")
print(round(tapply(brief$loneliness, brief$age_group, mean), 2))
cat("\nQ3. Number of people in each age group:\n")
print(table(brief$age_group))

# =============================================================================
# SECTION 3, ACTIVITY 4: smoothers, facets, colour and saving (Activity 6.5 on the lesson page)
# =============================================================================
set.seed(2021)
ggplot(brief, aes(x = age, y = loneliness)) +
  geom_jitter(height = 0.2, alpha = 0.1) +
  geom_smooth(method = "loess") +             # a smooth curve through the average score
  labs(x = "Age (years)", y = "Loneliness score (3 to 9)")

ggplot(brief, aes(x = loneliness)) +
  geom_histogram(binwidth = 1, fill = "#0B7B6B", colour = "white") +
  facet_wrap(~ age_group) +                   # one panel for each age group
  labs(x = "Loneliness score (3 to 9)", y = "Number of people")

set.seed(2021)
ggplot(brief, aes(x = support, y = loneliness)) +
  geom_jitter(height = 0.2, alpha = 0.2) +
  geom_smooth(method = "lm") +                # a straight line in each panel
  facet_grid(gender ~ age_group) +            # rows: gender; columns: age group
  labs(x = "Social support (1 to 7)", y = "Loneliness score (3 to 9)")

table(brief$gender, brief$age_group)         # people in each panel

p <- ggplot(brief, aes(x = age_group, y = loneliness, fill = gender)) +
  geom_boxplot() +
  scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73")) +  # Okabe-Ito colours
  labs(x = "Age group", y = "Loneliness score (3 to 9)", fill = "Gender") +
  theme_minimal(base_size = 12)               # a plain theme with larger text
p                                             # show the plot
ggsave("loneliness_age_gender.png", p, width = 7, height = 4.5, dpi = 300)

# --- Answers ---------------------------------------------------------------
cat("\nQ1. The loess curve at selected ages, and people at the two ends:\n")
lo <- loess(loneliness ~ age, data = brief)   # the same smoother as geom_smooth()
ages <- c(16, 30, 55, 70, 80)
print(round(setNames(predict(lo, data.frame(age = ages)), paste("age", ages)), 2))
print(c(under_20 = sum(brief$age < 20), aged_80_plus = sum(brief$age >= 80)))
cat("\nQ2. Slope of loneliness on support within each age group:\n")
slopes <- sapply(levels(brief$age_group), function(g)
  coef(lm(loneliness ~ support, data = brief[brief$age_group == g, ]))[["support"]])
print(round(slopes, 2))
cat("\nLoneliness scores within each age group (row percentages):\n")
print(round(100 * prop.table(table(brief$age_group, brief$loneliness), 1), 1))
# Q3. Okabe-Ito colours stay distinct with colour-vision deficiency; width and
#     height are in inches and dpi = 300 gives 2,100 by 1,350 pixels.

# =============================================================================
# SECTION 4, ACTIVITY 5: a correlation matrix and a heatmap (Activity 6.6 on the lesson page)
# =============================================================================
vars <- c("loneliness", "support", "depression", "anxiety", "life_sat", "age")
cor_mat <- round(cor(brief[, vars]), 2)     # Pearson correlations, rounded
cor_mat

cor_long <- as.data.frame(as.table(cor_mat))   # one row for each pair
names(cor_long) <- c("var1", "var2", "r")
head(cor_long)

ggplot(cor_long, aes(x = var1, y = var2, fill = r)) +
  geom_tile(colour = "white") +                     # one coloured square per pair
  geom_text(aes(label = r), size = 3.5) +           # print the value on each square
  scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       limits = c(-1, 1)) +         # blue negative, red positive
  labs(x = NULL, y = NULL, fill = "r") +
  theme_minimal()

# --- Answers ---------------------------------------------------------------
cat("\nQ1. Correlations with loneliness, strongest first:\n")
r_lonely <- cor_mat["loneliness", vars != "loneliness"]
print(r_lonely[order(-abs(r_lonely))])
# Q2. r = 0.07 for age, yet the loess curve (Activity 6.5) bends: r measures only
#     straight-line association.
cat("\nSpearman correlation of loneliness and support (reading):\n")
print(round(cor(brief$loneliness, brief$support, method = "spearman"), 2))
# Q3. A negative correlation (-0.40) from one survey wave cannot show cause.

# =============================================================================
# SECTION 4, ACTIVITY 6: marginal histograms and the combined figure (Activity 6.7 on the lesson page)
# =============================================================================
# install.packages(c("ggExtra", "patchwork"))   # run once, if not yet installed
library(ggExtra)     # ggMarginal() adds histograms to the margins of a plot
library(patchwork)   # + and / combine ggplot2 plots into one figure
set.seed(2021)
p_scatter <- ggplot(brief, aes(x = support, y = loneliness)) +
  geom_jitter(height = 0.2, width = 0, alpha = 0.15) +
  geom_smooth(method = "lm", colour = "#CC0033") +
  labs(x = "Social support (1 to 7)", y = "Loneliness score (3 to 9)")
ggMarginal(p_scatter, type = "histogram")    # a histogram on each axis

p_hist <- ggplot(brief, aes(x = loneliness)) +
  geom_histogram(binwidth = 1, fill = "#0B7B6B", colour = "white") +
  labs(x = "Loneliness score (3 to 9)", y = "Number of people")
p_box <- ggplot(brief, aes(x = age_group, y = loneliness)) +
  geom_boxplot(fill = "#E6F3F0") +
  labs(x = "Age group", y = "Loneliness score (3 to 9)")
p_heat <- ggplot(cor_long, aes(x = var1, y = var2, fill = r)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = r), size = 3) +
  scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       limits = c(-1, 1)) +
  labs(x = NULL, y = NULL, fill = "r")

figure <- (p_hist + p_box) / (p_scatter + p_heat) +   # + side by side, / stacked
  plot_annotation(tag_levels = "A",                     # label the panels A to D
    title = "Loneliness and social connection, CSCS 2021 (n = 3,083)")
figure
ggsave("briefing_figure.png", figure, width = 11, height = 8.5, dpi = 300)

# --- Answers ---------------------------------------------------------------
cat("\nQ1. The most common social support scores (spikes at whole numbers):\n")
print(head(sort(table(round(brief$support, 2)), decreasing = TRUE), 5))
# Q2. + places plots side by side, / stacks the rows; tag_levels = "A" labels
#     the panels A to D.
cat("\nQ3. Plain labels for the heatmap, one theme for every panel:\n")
nice <- c(loneliness = "Loneliness", support = "Social support",
          depression = "Depression", anxiety = "Anxiety",
          life_sat = "Life satisfaction", age = "Age")
levels(cor_long$var1) <- nice[levels(cor_long$var1)]
levels(cor_long$var2) <- nice[levels(cor_long$var2)]
print(levels(cor_long$var1))
p_heat <- ggplot(cor_long, aes(x = var1, y = var2, fill = r)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = r), size = 3) +
  scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       limits = c(-1, 1)) +
  labs(x = NULL, y = NULL, fill = "Pearson r")
figure_final <- ((p_hist + p_box) / (p_scatter + p_heat) +
  plot_annotation(tag_levels = "A")) & theme_minimal()
figure_final