# =============================================================================
# HSCI 841 Qualitative Research Methods and Analysis in Public Health  -  Lesson 8: Content Analysis
# Answer key for the in-lesson R activities
# Data file(s): term projects/HSCI_841/taguette_export_week8.csv,
#               term projects/HSCI_841/codebook_week8.csv,
#               term projects/HSCI_841/participant_metadata.csv,
#               term projects/HSCI_841/transcripts/P*.txt
#               (download from the lesson page; paths are relative to the course
#                repository root, so set your working directory there first)
# Packages: tidyverse, readtext, quanteda, quanteda.textstats, quanteda.textplots
#           (install once with install.packages(c("tidyverse","readtext","quanteda",
#            "quanteda.textstats","quanteda.textplots")))
# Reproduces every code block in the lesson and reports the numbers the page quotes.
# =============================================================================

suppressPackageStartupMessages({
  library(tidyverse)
  library(readtext)
})

dir.create("figures", showWarnings = FALSE)   # the ggsave() calls below write here

# ==== Section 1: Load Taguette export and merge participant metadata ====

# Read the Taguette export, one row per (passage, code) pair
codings <- read_csv("term projects/HSCI_841/taguette_export_week8.csv", show_col_types = FALSE)
glimpse(codings)
# Columns: id, document, tag (= code), content (= passage text), start, end

# Pull participant_id and pseudonym out of the document filename
codings <- codings |>
  mutate(
    participant_id = str_extract(document, "P[0-9]+"),
    pseudonym      = str_extract(document, "(?<=_)[A-Za-z]+(?=\\.txt)")
  )

# Read the participant metadata table (created from transcript headers)
# Columns: participant_id, pseudonym, age, gender, life_stage, caregiver, immigrant
participants <- read_csv("term projects/HSCI_841/participant_metadata.csv", show_col_types = FALSE)

# NOTE (answer key): the join key is participant_id AND pseudonym. Joining on
# participant_id alone leaves two pseudonym columns (pseudonym.x, pseudonym.y),
# which breaks any later code that refers to `pseudonym`.
codings <- codings |> left_join(participants, by = c("participant_id", "pseudonym"))

glimpse(codings)

# The codebook that defines the 11 codes
codebook <- read_csv("term projects/HSCI_841/codebook_week8.csv", show_col_types = FALSE)
cat("\nCodebook:", nrow(codebook), "codes;",
    n_distinct(codings$id), "highlights;", nrow(codings), "(passage, code) rows\n")

# ==== Section 2: Build the codes x cases matrix ====

# Binary indicator matrix: 1 if code appeared in transcript, 0 otherwise
code_matrix_binary <- codings |>
  distinct(participant_id, tag) |>
  mutate(present = 1) |>
  pivot_wider(
    names_from  = tag,
    values_from = present,
    values_fill = 0
  )

# Count matrix: number of times each code appeared in each transcript
code_matrix_counts <- codings |>
  count(participant_id, tag) |>
  pivot_wider(
    names_from  = tag,
    values_from = n,
    values_fill = 0
  )

# Add the participant-level variables for subgroup analysis
code_matrix_binary <- code_matrix_binary |> left_join(participants, by = "participant_id")

print(code_matrix_binary)

# ==== Section 3: Compute frequency tables: overall and by subgroup ====

# Overall code frequency (n transcripts in which each code appeared)
code_freq_overall <- code_matrix_binary |>
  summarise(across(where(is.numeric) & !c(age), sum)) |>
  pivot_longer(everything(), names_to = "code", values_to = "n_transcripts") |>
  mutate(pct = round(100 * n_transcripts / nrow(code_matrix_binary), 1)) |>
  arrange(desc(n_transcripts))
print(code_freq_overall, n = 20)

# By caregiver status: 2 x 2 contingency for each code
code_freq_by_caregiver <- codings |>
  distinct(participant_id, tag, caregiver) |>
  count(tag, caregiver) |>
  pivot_wider(names_from = caregiver, values_from = n, values_fill = 0)
print(code_freq_by_caregiver, n = 20)

# By age band
code_matrix_binary <- code_matrix_binary |>
  mutate(age_band = case_when(
    age < 30 ~ "18-29",
    age < 50 ~ "30-49",
    age < 65 ~ "50-64",
    TRUE     ~ "65+"
  ))

code_by_age <- code_matrix_binary |>
  group_by(age_band) |>
  summarise(across(starts_with("loneliness-"), sum), n_in_band = n())
print(code_by_age)

# The numbers the page quotes in Section 3.1 and 3.2
cat("\nMost prevalent code:", code_freq_overall$code[1], "in",
    code_freq_overall$n_transcripts[1], "of 20 transcripts\n")
cat("Corpus split:", sum(participants$caregiver == "caregiver"), "caregivers,",
    sum(participants$caregiver == "non-caregiver"), "non-caregivers\n")

# ==== Section 4: Bar chart of code frequencies by life-stage subgroup ====

library(ggplot2)

# Long-format frequency table for plotting
plot_data <- codings |>
  distinct(participant_id, tag, caregiver) |>
  count(tag, caregiver) |>
  left_join(
    participants |> count(caregiver, name = "n_in_group"),
    by = "caregiver"
  ) |>
  mutate(pct = 100 * n / n_in_group)

ggplot(plot_data, aes(x = reorder(tag, pct), y = pct, fill = caregiver)) +
  geom_col(position = "dodge") +
  coord_flip() +
  scale_fill_manual(values = c("caregiver" = "#CC0033", "non-caregiver" = "#0B7B6B")) +
  labs(
    x = NULL,
    y = "% of subgroup with code present",
    fill = "Caregiver status",
    title = "Code prevalence by caregiver status (loneliness corpus, n=20)"
  ) +
  theme_minimal(base_size = 12) +
  theme(panel.grid.major.y = element_blank())

ggsave("figures/code_prevalence_by_caregiver.png", width = 8, height = 5, dpi = 300)

# ==== Section 5: Chi-squared and Fisher's exact tests on codes x subgroup ====

# NOTE (answer key): the page's earlier draft used a code called
# shame-prevents-disclosure. Every code in the week 8 codebook starts with
# "loneliness-", and the code that actually separates the two subgroups is
# loneliness-as-sole-responsibility, so that is the worked example here.

# Build the 2 x 2 contingency table
resp_table <- table(
  caregiver      = code_matrix_binary$caregiver,
  responsibility = code_matrix_binary$`loneliness-as-sole-responsibility`
)
print(resp_table)

# Chi-squared with Yates' continuity correction (default for 2 x 2)
print(suppressWarnings(chisq.test(resp_table)))

# Fisher's exact, the more defensible choice for small expected counts
print(fisher.test(resp_table))

# Loop across all codes to produce a table of p-values
code_names <- setdiff(colnames(code_matrix_binary),
                       c("participant_id", "age", "gender", "life_stage",
                         "caregiver", "immigrant", "pseudonym", "age_band"))

test_results <- map_dfr(code_names, function(cd) {
  tbl <- table(code_matrix_binary$caregiver, code_matrix_binary[[cd]])
  ft  <- fisher.test(tbl)
  tibble(
    code     = cd,
    # NOTE (answer key): a code present in every caregiver (or in none) produces a
    # table without a "0" or a "1" column, so read the counts defensively.
    n_caregiver_pos    = if ("1" %in% colnames(tbl)) tbl["caregiver", "1"] else 0L,
    n_noncaregiver_pos = if ("1" %in% colnames(tbl)) tbl["non-caregiver", "1"] else 0L,
    odds_ratio = unname(ft$estimate),
    p_value    = ft$p.value
  )
}) |> arrange(p_value)

print(test_results, n = 20)

cat("\nOnly", sum(test_results$p_value < 0.05),
    "code differs by caregiver status at p < 0.05:",
    test_results$code[1], "(Fisher p =", signif(test_results$p_value[1], 3), ")\n")

# ==== Section 6: Keyness analysis: which words distinguish caregiver from non-caregiver transcripts? ====

library(quanteda)
library(quanteda.textstats)
library(quanteda.textplots)

# Build the quanteda corpus from the transcripts (as earlier in the course)
loneliness_rt <- readtext("term projects/HSCI_841/transcripts/P*.txt",
                          docvarsfrom = "filenames",
                          docvarnames = c("participant_id", "pseudonym"),
                          dvsep = "_")
loneliness_corpus <- corpus(loneliness_rt)

# Attach caregiver status as a document variable
docvars(loneliness_corpus, "caregiver") <- participants$caregiver[
  match(docvars(loneliness_corpus, "participant_id"),
        participants$participant_id)
]

# Tokenise, lowercase, remove stopwords
loneliness_tokens <- tokens(loneliness_corpus,
                            remove_punct   = TRUE,
                            remove_numbers = TRUE) |>
  tokens_tolower() |>
  tokens_remove(stopwords("en"))

# Build the document-feature matrix and group by caregiver status
loneliness_dfm <- dfm(loneliness_tokens) |>
  dfm_group(groups = docvars(loneliness_corpus, "caregiver"))

# Keyness test: caregiver subcorpus as target, non-caregiver as reference
keyness <- textstat_keyness(loneliness_dfm,
                            target = "caregiver",
                            measure = "lr")  # log-likelihood ratio
print(head(keyness, 30))
cat("\nWords distinctive of the non-caregiver transcripts:\n")
print(tail(keyness, 10))

# Plot top 20 keywords on each side
textplot_keyness(keyness, n = 20)
ggsave("figures/keyness_caregiver_vs_noncaregiver.png", width = 8, height = 6, dpi = 300)

cat("\nDone: Lesson 8 answer key completed.\n")