# =============================================================================
# HSCI 841 Qualitative Research Methods and Analysis in Public Health  -  Lesson 12: Computational Text and LLM Analysis
# Answer key for the in-lesson R activities
# Data file(s): term projects/HSCI_841/transcripts/P*.txt,
#               term projects/HSCI_841/freelist_what_helps.csv,
#               term projects/HSCI_841/pilesort_relationships.csv,
#               term projects/HSCI_841/consensus_loneliness.csv,
#               term projects/HSCI_841/outputs/wk12_hand_codings.csv,
#               term projects/HSCI_841/outputs/wk12_llm_codings.csv
#               (download from the lesson page; paths are relative to the course
#                repository root, so set your working directory there first)
# Packages: tidyverse, quanteda, quanteda.textstats, quanteda.textplots, readtext,
#           tidytext, igraph, tidygraph, ggraph, irr
#           (install once with install.packages(c("tidyverse","quanteda",
#            "quanteda.textstats","quanteda.textplots","readtext","tidytext",
#            "igraph","tidygraph","ggraph","irr")))
# Reproduces every code block in the lesson and writes the outputs/wk12_*.csv files.
# =============================================================================

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

# every write_csv() below lands here
dir.create("term projects/HSCI_841/outputs", showWarnings = FALSE, recursive = TRUE)
out <- function(f) file.path("term projects/HSCI_841/outputs", f)

# ==== Section 1: KWIC analysis of the loneliness corpus ====

loneliness_rt <- readtext(
  "term projects/HSCI_841/transcripts/P*.txt",
  docvarsfrom = "filenames",
  docvarnames = c("participant_id", "pseudonym"),
  dvsep = "_"
)
loneliness_corpus <- corpus(loneliness_rt)
loneliness_tokens <- tokens(loneliness_corpus, remove_punct = TRUE) |>
  tokens_tolower()

# KWIC for "chair" across the 20 transcripts, 6 words of context on either side
chair_kwic <- kwic(loneliness_tokens, pattern = phrase("chair*"), window = 6)
print(chair_kwic, max_ndoc = 20)

# Convert to a tibble so you can write it out for close reading
chair_df <- as_tibble(chair_kwic)
write_csv(chair_df, out("wk12_kwic_chair.csv"))

# Contrastive KWIC: "alone" vs "lonely" -- different words, different meanings
alone_kwic  <- kwic(loneliness_tokens, pattern = "alone",  window = 5)
lonely_kwic <- kwic(loneliness_tokens, pattern = "lonely", window = 5)

cat("alone:  ", nrow(as_tibble(alone_kwic)),  "occurrences\n")
cat("lonely: ", nrow(as_tibble(lonely_kwic)), "occurrences\n")

# KWIC for multiple terms at once
candidate_terms <- c("chair", "empty", "hollow", "fading", "tired", "wahda")
multi_kwic <- kwic(loneliness_tokens, pattern = candidate_terms, window = 5)
print(head(multi_kwic, 20))

# ==== Section 2: Word frequencies, type-token ratios, lexical diversity ====

library(quanteda.textstats)

loneliness_tokens_nostop <- tokens_remove(loneliness_tokens, stopwords("en"))
loneliness_dfm <- dfm(loneliness_tokens_nostop)

top50 <- topfeatures(loneliness_dfm, n = 50)
print(top50)

freq_tbl <- textstat_frequency(loneliness_dfm, n = 100)
print(head(freq_tbl, 20))
write_csv(freq_tbl, out("wk12_word_frequencies.csv"))

ttr_basic <- textstat_lexdiv(loneliness_dfm, measure = "TTR")
print(ttr_basic)

ttr_mattr <- textstat_lexdiv(
  loneliness_tokens_nostop,
  measure = c("TTR", "MATTR", "MSTTR"),
  MATTR_window = 100,
  MSTTR_segment = 100
)
print(ttr_mattr)

ttr_mattr |>
  as_tibble() |>
  ggplot(aes(reorder(document, MATTR), MATTR)) +
  geom_col(fill = "#0B7B6B") +
  coord_flip() +
  labs(x = "", y = "MATTR (length-corrected lexical diversity)",
       title = "Lexical diversity across the 20 loneliness transcripts") +
  theme_minimal()

# ==== Section 3: TF-IDF for case-distinctive vocabulary ====

loneliness_tfidf <- dfm_tfidf(loneliness_dfm)

# Top 10 most-distinctive words per transcript
# NOTE (answer key): apply() over a dfm needs a plain matrix, so convert first.
top_tfidf_per_doc <- apply(as.matrix(loneliness_tfidf), 1, function(x) {
  names(sort(x, decreasing = TRUE))[1:10]
})
print(top_tfidf_per_doc)

library(tidytext)
loneliness_long <- as_tibble(convert(loneliness_dfm, to = "data.frame")) |>
  pivot_longer(-doc_id, names_to = "word", values_to = "n") |>
  filter(n > 0)

loneliness_tfidf_tidy <- loneliness_long |>
  bind_tf_idf(word, doc_id, n)

loneliness_tfidf_tidy |>
  group_by(doc_id) |>
  slice_max(tf_idf, n = 5, with_ties = FALSE) |>
  arrange(doc_id, desc(tf_idf)) |>
  print(n = 100)

# ==== Section 4: Keyness analysis: older vs younger participants ====

library(quanteda.textplots)

older   <- c("P05", "P11", "P17", "P20")   # 70s and 80s
younger <- c("P01", "P10", "P12", "P19")   # 20s

docvars(loneliness_dfm, "age_group") <- case_when(
  docvars(loneliness_dfm, "participant_id") %in% older   ~ "older",
  docvars(loneliness_dfm, "participant_id") %in% younger ~ "younger",
  TRUE                                                   ~ NA_character_
)

keyness_dfm <- dfm_subset(loneliness_dfm, !is.na(age_group))

keyness_older <- textstat_keyness(
  keyness_dfm,
  target  = docvars(keyness_dfm, "age_group") == "older",
  measure = "lr"   # log-likelihood ratio; "chi2" is also available
)

print(head(keyness_older, 30))  # top 30 words distinctive of OLDER participants
print(tail(keyness_older, 30))  # top 30 words distinctive of YOUNGER participants

textplot_keyness(keyness_older, n = 20,
  color = c("#0B7B6B", "#CC0033"),
  labelsize = 3.5
)

write_csv(as_tibble(keyness_older), out("wk12_keyness_older_v_younger.csv"))

# ==== Section 5: Collocations: empty chairs, quiet apartments, group chats ====

collocs <- textstat_collocations(
  loneliness_tokens_nostop,
  size = 2:3,
  min_count = 4
)

print(collocs |>
  as_tibble() |>
  arrange(desc(z)) |>
  head(40))

alone_neighbors <- tokens_select(loneliness_tokens_nostop,
  pattern = "alone",
  window = 4,
  selection = "keep"
) |>
  dfm() |>
  topfeatures(n = 30)
print(alone_neighbors)

lonely_neighbors <- tokens_select(loneliness_tokens_nostop,
  pattern = "lonely",
  window = 4,
  selection = "keep"
) |>
  dfm() |>
  topfeatures(n = 30)
print(lonely_neighbors)

print(setdiff(names(alone_neighbors),  names(lonely_neighbors)))
print(setdiff(names(lonely_neighbors), names(alone_neighbors)))

# ==== Section 6: Smith's salience from a free-listing dataset ====

# NOTE (answer key): AnthroTools is not on CRAN for R 4.3, so Smith's salience is
# computed directly. For one mention at rank R in a list of length L the salience
# is (L - R + 1) / L; the item's Smith's S is the mean of that over all
# respondents, counting 0 for the respondents who did not list it.
free_lists <- read_csv("term projects/HSCI_841/freelist_what_helps.csv",
                       show_col_types = FALSE)
n_respondents <- n_distinct(free_lists$Subj)

salience <- free_lists |>
  group_by(Subj) |>
  mutate(list_length = n(),
         Smith.Salience = (list_length - Order + 1) / list_length) |>
  ungroup()
print(head(salience))

salience_summary <- salience |>
  group_by(CODE) |>
  summarise(SmithsS   = sum(Smith.Salience) / n_respondents,
            n_lists   = n(),
            mean_rank = round(mean(Order), 2),
            .groups   = "drop") |>
  arrange(desc(SmithsS))
print(head(salience_summary, 20))

cat("\n", n_respondents, "free lists,", n_distinct(free_lists$CODE), "distinct items.",
    "The most salient item is", salience_summary$CODE[1],
    "(Smith's S =", round(salience_summary$SmithsS[1], 3), ").\n")
write_csv(salience_summary, out("wk12_freelist_salience.csv"))

# ==== Section 7: MDS on aggregate pile-sort data ====

# Pile-sort data: 10 kinds of relationship, 20 sorters, each row one partition
piles <- read_csv("term projects/HSCI_841/pilesort_relationships.csv",
                  show_col_types = FALSE)

items <- setdiff(names(piles), "participant")
M <- matrix(0, length(items), length(items),
            dimnames = list(items, items))
for (p in 1:nrow(piles)) {
  for (i in items) for (j in items) {
    if (piles[[i]][p] == piles[[j]][p]) M[i, j] <- M[i, j] + 1
  }
}
print(M)

# Convert similarity to dissimilarity (max - similarity)
D <- max(M) - M
diag(D) <- 0

mds_fit <- cmdscale(as.dist(D), k = 2)

print(tibble(item = rownames(mds_fit), x = mds_fit[, 1], y = mds_fit[, 2]) |>
  ggplot(aes(x, y, label = item)) +
  geom_point(size = 4, color = "#0B7B6B") +
  geom_text(aes(label = item), vjust = -1.2) +
  labs(title = "MDS of pile-sort co-membership: kinds of social relationship",
       x = "Dimension 1", y = "Dimension 2") +
  theme_minimal())

# ==== Section 8: Consensus analysis from the agreement matrix ====

# NOTE (answer key): the AnthroTools ConsensusPipeline() is replaced by the two
# quantities it reports, computed in base R: the ratio of the first to the second
# eigenvalue of the respondent-by-respondent agreement matrix (>= 3 supports a
# single shared culture) and the first eigenvector, which is the competence score.
consensus_raw <- read_csv("term projects/HSCI_841/consensus_loneliness.csv",
                          show_col_types = FALSE)
answers <- as.matrix(consensus_raw[, -1])       # items x participants, 1 = yes
rownames(answers) <- consensus_raw$item
print(dim(answers))

A <- t(answers)                                  # participants x items
n_items <- ncol(A)

# proportion of items on which each pair of participants gives the same answer
agreement_matrix <- (A %*% t(A) + (1 - A) %*% t(1 - A)) / n_items
print(round(agreement_matrix[1:6, 1:6], 2))

ev <- eigen(agreement_matrix, symmetric = TRUE)
eigen_ratio <- ev$values[1] / ev$values[2]

competence <- ev$vectors[, 1]
if (mean(competence) < 0) competence <- -competence      # sign is arbitrary
competence <- competence / max(competence)               # scale to a 0-1 range
names(competence) <- rownames(A)

# competence-weighted vote gives the model's "cultural answer" for each item
cultural_answers <- as.integer(colSums(A * competence) / sum(competence) > 0.5)
names(cultural_answers) <- colnames(A)

cat("Eigenvalue ratio:", round(eigen_ratio, 2), "\n")
cat("Mean competence:", round(mean(competence), 2), "\n")
cat("Cultural answers:\n"); print(cultural_answers)
cat(ifelse(eigen_ratio >= 3,
           "The ratio is above 3, so a single shared cultural model fits these answers.\n",
           "The ratio is below 3, so there is no single shared cultural model here.\n"))

write_csv(tibble(participant = names(competence), competence = round(competence, 3)),
          out("wk12_consensus_competence.csv"))

# ==== Section 9: Building a feature co-occurrence matrix with quanteda ====

top50_words <- names(topfeatures(loneliness_dfm, n = 50))

top_tokens <- tokens_select(loneliness_tokens_nostop,
  pattern = top50_words,
  selection = "keep"
)

fcm_window <- fcm(top_tokens, context = "window", window = 5, tri = FALSE)
print(dim(fcm_window))          # 50 x 50
# NOTE (answer key): topfeatures() on an fcm is defunct in quanteda 4, so the
# total co-occurrence weight per word is summed directly.
print(sort(colSums(fcm_window), decreasing = TRUE)[1:10])

fcm_doc <- fcm(top_tokens, context = "document", tri = FALSE)

set.seed(841)
textplot_network(fcm_window,
  min_freq = 0.5,
  edge_color = "#0B7B6B",
  edge_alpha = 0.6,
  vertex_labelfont = "sans",
  vertex_labelsize = 3.5
)

# ==== Section 10: Convert FCM to igraph and compute centrality ====

library(igraph)

# NOTE (answer key): quanteda 4 no longer supplies an as.igraph() method for an
# fcm, so the graph is built from the co-occurrence matrix itself. mode = "max"
# makes it undirected, and the isolated vertices are dropped afterwards.
fcm_mat <- as.matrix(fcm_window)
diag(fcm_mat) <- 0
g <- graph_from_adjacency_matrix(fcm_mat, mode = "max", weighted = TRUE, diag = FALSE)
g <- delete_vertices(g, degree(g) == 0)
print(summary(g))
print(vcount(g))   # number of vertices (words)
print(ecount(g))   # number of edges

deg <- degree(g)
print(sort(deg, decreasing = TRUE)[1:15])

str <- strength(g)
print(sort(str, decreasing = TRUE)[1:15])

btw <- betweenness(g, weights = 1 / E(g)$weight)
print(sort(btw, decreasing = TRUE)[1:15])

cls <- closeness(g, weights = 1 / E(g)$weight)
print(sort(cls, decreasing = TRUE)[1:15])

eig <- eigen_centrality(g)$vector
print(sort(eig, decreasing = TRUE)[1:15])

centrality_tbl <- tibble(
  word        = V(g)$name,
  degree      = deg,
  strength    = str,
  betweenness = btw,
  closeness   = cls,
  eigenvector = eig
) |> arrange(desc(eigenvector))
print(centrality_tbl, n = 20)
write_csv(centrality_tbl, out("wk12_semantic_centrality.csv"))

# ==== Section 11: Community detection with Louvain ====

set.seed(841)  # for reproducibility
comm_louvain <- cluster_louvain(g, weights = E(g)$weight)

print(membership(comm_louvain))
print(sizes(comm_louvain))
print(modularity(comm_louvain))
print(length(comm_louvain))

community_tbl <- tibble(
  word      = V(g)$name,
  community = as.factor(membership(comm_louvain))
) |>
  arrange(community, word)
print(community_tbl, n = 50)

comm_walktrap <- cluster_walktrap(g, weights = E(g)$weight, steps = 4)
print(modularity(comm_walktrap))
print(length(comm_walktrap))

print(compare(membership(comm_louvain), membership(comm_walktrap), method = "rand"))

# ==== Section 12: Force-directed visualization with ggraph ====

library(ggraph)
library(tidygraph)

# NOTE (answer key): inside activate(nodes) the edge attribute `weight` is not in
# scope, so the betweenness weights are taken from the edge table with .E()$weight.
g_tidy <- as_tbl_graph(g) |>
  activate(nodes) |>
  mutate(
    community = as.factor(membership(comm_louvain)),
    degree    = centrality_degree(),
    btwn      = centrality_betweenness(weights = 1 / .E()$weight)
  )

set.seed(841)
print(ggraph(g_tidy, layout = "fr") +
  geom_edge_link(aes(width = weight, alpha = weight),
                 edge_colour = "#888", show.legend = FALSE) +
  scale_edge_width(range = c(0.1, 1.5)) +
  scale_edge_alpha(range = c(0.1, 0.7)) +
  geom_node_point(aes(size = degree, colour = community)) +
  geom_node_text(aes(label = name), repel = TRUE, size = 3.2) +
  scale_colour_brewer(palette = "Set2") +
  labs(title = "Semantic network of the loneliness corpus (top 50 content words)",
       subtitle = "Edges = sliding-window co-occurrence (window = 5); communities by Louvain",
       colour = "Community", size = "Degree") +
  # NOTE (answer key): theme_graph() defaults to the Arial Narrow font, which is
  # not installed on every machine; base_family = "sans" always works.
  theme_graph(base_family = "sans"))

# ==== Section 13: Network-level statistics ====

print(edge_density(g))
print(mean_distance(g, weights = NA))
print(transitivity(g, type = "global"))
print(diameter(g, weights = NA))

set.seed(841)
g_random <- sample_gnm(vcount(g), ecount(g), directed = FALSE)
print(transitivity(g_random, type = "global"))
print(mean_distance(g_random))

# A real semantic network should have substantially higher clustering than the
# random graph and a similar or shorter mean distance: the small-world property.

# ==== Section 14: Compute human-LLM agreement with Krippendorff's alpha ====

library(irr)

# Read both coding sets (long format: doc, passage_id, code)
hand <- read_csv(out("wk12_hand_codings.csv"), show_col_types = FALSE)
llm  <- read_csv(out("wk12_llm_codings.csv"),  show_col_types = FALSE)
cat("hand:", nrow(hand), "rows;  llm:", nrow(llm), "rows\n")

# Join by passage_id (you have constructed this in step 3 above)
agreement <- hand |>
  rename(code_hand = code) |>
  left_join(rename(llm, code_llm = code), by = c("doc", "passage_id"))

# NOTE (answer key): both files carry a `passage` column, so the join renames them
# passage.x and passage.y. Keep the hand coder's text and drop the duplicate.
agreement <- agreement |>
  select(doc, passage_id, code_hand, code_llm, passage = passage.x)

# Encode codes as integers (irr requires numeric)
all_codes <- unique(c(agreement$code_hand, agreement$code_llm))
agreement <- agreement |>
  mutate(code_hand_int = match(code_hand, all_codes),
         code_llm_int  = match(code_llm,  all_codes))

# Build the 2-row matrix that irr expects (rows = coders, columns = passages)
codings_matrix <- rbind(
  hand = agreement$code_hand_int,
  llm  = agreement$code_llm_int
)

# Krippendorff's alpha (nominal codes)
print(kripp.alpha(codings_matrix, method = "nominal"))

# Cohen's kappa for direct comparison
# NOTE (answer key): kappa2() needs complete pairs, so drop the passages the LLM
# did not return before computing it.
complete_pairs <- codings_matrix[, complete.cases(t(codings_matrix)), drop = FALSE]
print(kappa2(t(complete_pairs), weight = "unweighted"))

# Confusion matrix: which codes does the LLM systematically disagree on?
print(agreement |>
  count(code_hand, code_llm) |>
  pivot_wider(names_from = code_llm, values_from = n, values_fill = 0),
  width = 200)

cat("\nRaw agreement:",
    round(mean(agreement$code_hand == agreement$code_llm, na.rm = TRUE), 3),
    " matched passages:", sum(!is.na(agreement$code_llm)), "of", nrow(agreement), "\n")

# ==== Section 15: Hallucination audit: verbatim-quote verification ====

library(stringr)

# Sample 30 random LLM-coded passages
set.seed(841)
audit_sample <- llm |> slice_sample(n = 30)

# For each row, load the source transcript and check whether the passage is verbatim
audit_sample <- audit_sample |>
  rowwise() |>
  mutate(
    transcript_text = paste(readLines(file.path(
      "term projects/HSCI_841/transcripts", paste0(doc, ".txt"))),
      collapse = "\n"),
    is_verbatim = str_detect(transcript_text,
                              fixed(str_squish(passage)))
  ) |>
  ungroup()

# Summary: hallucination rate
cat("proportion verbatim:", mean(audit_sample$is_verbatim), "(target: >= 0.95)\n")
print(audit_sample |>
  filter(!is_verbatim) |>
  select(doc, code, passage))

# Save the audit log for the appendix
write_csv(audit_sample |> select(-transcript_text), out("wk12_hallucination_audit.csv"))

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