1 Introduction

1.1 Purpose of the project

This project applies text mining techniques to analyze mental health discussions on Reddit, specifically examining how online mental health communities evolved during the COVID-19 pandemic. The analysis uses selected files from the Reddit Mental Health Dataset (Low et al., 2020) published on Zenodo: https://zenodo.org/records/3941387, which contains posts from 28 subreddits including 15 mental health support groups.

Selected data include:

  • ADHD subreddits from three timeframes:
    • 2018: Jan 1 to April 20, 2018 (control for seasonal fluctuations)
    • 2019: Jan 1 to April 20, 2019 (control for seasonal fluctuations)
    • post: Jan 1 to April 20, 2020 (mid-pandemic period)

1.2 Objectives

  1. Exploratory Data Analysis - understand the distribution and characteristics of posts
  2. Clustering - group similar terms to identify common themes and patterns
  3. Topic Modeling - discover latent topics in ADHD discussions using LDA
  4. Sentiment Analysis - examine emotional tone and its evolution over time periods
  5. Text Networks - visualize word co-occurrence patterns using bigrams

1.3 Assumptions

  1. Posts from the same timeframe in 2018 and 2019 serve as appropriate control periods, accounting for seasonal variations in mental health discussions
  2. The mid-pandemic period (Jan-April 2020) represents a distinct change point in online mental health discourse
  3. Text preprocessing (stemming, stopword removal) preserves meaningful content while reducing noise
  4. Bigram networks effectively capture contextual relationships between mental health terms
  5. Dictionary-based sentiment analysis (Bing, AFINN, NRC lexicons) is valid for mental health discourse despite being for general purpose, not domain-specific.

2 Data

2.1 Files

Files used in this analysis:

  • adhd_2018_features_tfidf_256.csv
  • adhd_2019_features_tfidf_256.csv
  • adhd_post_features_tfidf_256.csv

2.2 Required packages

# Required packages
library(tm)           
library(SnowballC)    
library(lsa)          
library(tidyverse)    # zawiera dplyr, tidyr, ggplot2
library(tidytext)     
library(topicmodels)  
library(cluster)      
library(factoextra)   
library(fpc)          
library(igraph)       
library(ggraph)       
library(RColorBrewer) 
library(wordcloud)    
library(DT)           
library(reshape2)     
library(forcats)
library(gridExtra)

2.3 Loading data

# Load ADHD datasets
adhd_2018 <- read.csv('adhd_2018_features_tfidf_256.csv')
adhd_2019 <- read.csv('adhd_2019_features_tfidf_256.csv')
adhd_post <- read.csv('adhd_post_features_tfidf_256.csv')

2.4 Data structure

cat("Structure of adhd_2018 dataset:\n")
## Structure of adhd_2018 dataset:
cat("Dimensions:", dim(adhd_2018),"\n")
## Dimensions: 5110 350
cat("Column names (first 10):\n")
## Column names (first 10):
print(names(adhd_2018)[1:min(10, ncol(adhd_2018))])
##  [1] "subreddit"                   "author"                     
##  [3] "date"                        "post"                       
##  [5] "automated_readability_index" "coleman_liau_index"         
##  [7] "flesch_kincaid_grade_level"  "flesch_reading_ease"        
##  [9] "gulpease_index"              "gunning_fog_index"
cat("\nStructure of adhd_2019 dataset:\n")
## 
## Structure of adhd_2019 dataset:
cat("Dimensions:", dim(adhd_2019),"\n")
## Dimensions: 10223 350
cat("Column names (first 10):\n")
## Column names (first 10):
print(names(adhd_2019)[1:min(10, ncol(adhd_2019))])
##  [1] "subreddit"                   "author"                     
##  [3] "date"                        "post"                       
##  [5] "automated_readability_index" "coleman_liau_index"         
##  [7] "flesch_kincaid_grade_level"  "flesch_reading_ease"        
##  [9] "gulpease_index"              "gunning_fog_index"
cat("\nStructure of adhd_post dataset:\n")
## 
## Structure of adhd_post dataset:
cat("Dimensions:", dim(adhd_post),"\n")
## Dimensions: 13602 350
cat("Column names (first 10):\n")
## Column names (first 10):
print(names(adhd_post)[1:min(10, ncol(adhd_post))])
##  [1] "subreddit"                   "author"                     
##  [3] "date"                        "post"                       
##  [5] "automated_readability_index" "coleman_liau_index"         
##  [7] "flesch_kincaid_grade_level"  "flesch_reading_ease"        
##  [9] "gulpease_index"              "gunning_fog_index"

2.5 Combine datasets

adhd_2018_clean <- adhd_2018 %>%
  mutate(
    period = "control_2018",
    date = as.POSIXct(date, origin = "1970-01-01")
  ) %>%
  select(subreddit, period, author, date, post)

adhd_2019_clean <- adhd_2019 %>%
  mutate(
    period = "control_2019",
    date = as.POSIXct(date, origin = "1970-01-01")
  ) %>%
  select(subreddit, period, author, date, post)

adhd_post_clean <- adhd_post %>%
  mutate(
    period = "mid_pandemic",
    date = as.POSIXct(date, origin = "1970-01-01")
  ) %>%
  select(subreddit, period, author, date, post)

# combine periods
adhd_combined <- rbind(adhd_2018_clean, adhd_2019_clean, adhd_post_clean)

cat("Structure of adhd_combined dataset:\n")
## Structure of adhd_combined dataset:
cat("Dimensions:", dim(adhd_combined),"\n")
## Dimensions: 28935 5
str(adhd_combined)
## 'data.frame':    28935 obs. of  5 variables:
##  $ subreddit: chr  "adhd" "adhd" "adhd" "adhd" ...
##  $ period   : chr  "control_2018" "control_2018" "control_2018" "control_2018" ...
##  $ author   : chr  "Seftari" "LisaLoves2" "featherflutter" "vanillasky0" ...
##  $ date     : POSIXct, format: "2018-01-01" "2018-01-01" ...
##  $ post     : chr  "Lethargic/Depressed when off meds First I'll give some background to my medical history. I have struggled with "| __truncated__ "Concerta not working on the first day?! Update!: Thank you all for your insightful and kind responses, if not f"| __truncated__ "Comorbid anxiety and ADHD-PI Medication Question Hi, I was recently diagnosed with ADHD-Pi. I am an adult, and "| __truncated__ "Fist Day on Concerta 18mg UPDATE! Update!: Thank you all for your insightful and kind responses, if not for the"| __truncated__ ...
cat("\nBreakdown by period:\n")
## 
## Breakdown by period:
print(table(adhd_combined$period))
## 
## control_2018 control_2019 mid_pandemic 
##         5110        10223        13602
cat("\nFirst 3 rows of combined data:\n")
## 
## First 3 rows of combined data:
head(adhd_combined, 3)

3 ADHD dataset analysis

3.1 Preprocessing

# create corpus from post text
corpus_adhd <- VCorpus(VectorSource(adhd_combined$post))

# check size
print(corpus_adhd) 
## <<VCorpus>>
## Metadata:  corpus specific: 0, document level (indexed): 0
## Content:  documents: 28935
cat("\nExample document content:\n")
## 
## Example document content:
inspect(corpus_adhd[1])
## <<VCorpus>>
## Metadata:  corpus specific: 0, document level (indexed): 0
## Content:  documents: 1
## 
## [[1]]
## <<PlainTextDocument>>
## Metadata:  7
## Content:  chars: 761
# ensure encoding across the entire corpus
corpus_adhd <- tm_map(corpus_adhd, content_transformer(function(x) iconv(x, to = "UTF-8", sub = "byte")))

# function to replace characters with space
toSpace <- content_transformer(function (x, pattern) gsub(pattern, " ", x))

# Remove URLs and web remnants
corpus_adhd <- tm_map(corpus_adhd, toSpace, "@")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "@\\w+")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "\\|")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "[ \t]{2,}")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "(s?)(f|ht)tp(s?)://\\S+\\b")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "http\\w*")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "/")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "www")
corpus_adhd <- tm_map(corpus_adhd, toSpace, "~")

# text cleaning
corpus_adhd <- tm_map(corpus_adhd, content_transformer(tolower))
corpus_adhd <- tm_map(corpus_adhd, removePunctuation)
corpus_adhd <- tm_map(corpus_adhd, removeNumbers) 
corpus_adhd <- tm_map(corpus_adhd, removeWords, stopwords("english"))

# remove custom stopwords
corpus_adhd <- tm_map(corpus_adhd, removeWords, c("adhd", "im", "dont", "can", "get", "just", "like", "ive", "ill", "really"))

corpus_adhd <- tm_map(corpus_adhd, stripWhitespace)

# check
cat("\nCleaned document content:\n")
## 
## Cleaned document content:
corpus_adhd[[1]][[1]]
## [1] "lethargic depressed meds first give background medical history struggled depression years nowim currently bouts crying absolute consumption life every week months ago diagnosed add started taking mg xr adderall morning noon montg starting noticed longer bouts depression overall happier mom super anti meds blames everything bad adderall will take away randomly every time feel withdrawals bouts depression come back feel want lie die anyone else feel way "
# keep a copy of the corpus to use as dictionary for stem completion
corpus_adhd_copy <- corpus_adhd

# perform stemming
corpus_adhd_stemmed <- tm_map(corpus_adhd, stemDocument)

# check stemmed content
cat("Original:\n")
## Original:
corpus_adhd[[1]][[1]]
## [1] "lethargic depressed meds first give background medical history struggled depression years nowim currently bouts crying absolute consumption life every week months ago diagnosed add started taking mg xr adderall morning noon montg starting noticed longer bouts depression overall happier mom super anti meds blames everything bad adderall will take away randomly every time feel withdrawals bouts depression come back feel want lie die anyone else feel way "
cat("\n\nStemmed:\n")
## 
## 
## Stemmed:
corpus_adhd_stemmed[[1]][[1]]
## [1] "letharg depress med first give background medic histori struggl depress year nowim current bout cri absolut consumpt life everi week month ago diagnos add start take mg xr adderal morn noon montg start notic longer bout depress overal happier mom super anti med blame everyth bad adderal will take away random everi time feel withdraw bout depress come back feel want lie die anyon els feel way"

3.2 Tokenization and Document-Term Matrix

dtm_adhd <- DocumentTermMatrix(corpus_adhd_stemmed)

cat("DTM dimensions:\n")
## DTM dimensions:
print(dtm_adhd)
## <<DocumentTermMatrix (documents: 28935, terms: 34524)>>
## Non-/sparse entries: 1718786/997233154
## Sparsity           : 100%
## Maximal term length: 477
## Weighting          : term frequency (tf)
# reduce matrix sparsity - keep terms that appear in at least 1% of documents
dtm_sparse_adhd <- removeSparseTerms(dtm_adhd, sparse = 0.99)

cat("Reduced DTM dimensions:\n")
## Reduced DTM dimensions:
print(dtm_sparse_adhd)
## <<DocumentTermMatrix (documents: 28935, terms: 1019)>>
## Non-/sparse entries: 1378978/28105787
## Sparsity           : 95%
## Maximal term length: 12
## Weighting          : term frequency (tf)
# convert to matrix
dtm_adhd_m <- as.matrix(dtm_sparse_adhd)

3.3 Exploratory Data Analysis

3.3.1 Word frequency analysis

# count word frequencies
word_freq <- sort(colSums(dtm_adhd_m), decreasing = TRUE)
freq_df_adhd <- data.frame(word = names(word_freq), freq = word_freq, row.names = NULL)

cat("Top 20 most frequent terms:\n")
## Top 20 most frequent terms:
head(freq_df_adhd, 20)
# count word frequencies by period
freq_by_period <- adhd_combined %>%
  unnest_tokens(word, post) %>%
  mutate(word = str_replace_all(word, "[’']", "")) %>%
  mutate(word = str_trim(word)) %>%
  anti_join(stop_words, by = "word") %>%
  filter(!word %in% c("adhd", "im", "dont", "can", "just", "like", "get", "amp", "ve", "ive", "cant", "wont", "didnt", "x200b")) %>%
  filter(!str_detect(word, "^[0-9]*$")) %>%
  count(period, word, sort = TRUE) %>%
  ungroup()

cat("Top frequencies calculated by period.\n")
## Top frequencies calculated by period.
head(freq_by_period)

3.3.2 Visualization of Frequent Terms

# bar plot of top 20 terms
top_20_adhd <- freq_df_adhd %>%
  slice_max(freq, n = 20) %>%
  mutate(word = fct_reorder(word, freq))

ggplot(top_20_adhd, aes(x = freq, y = word)) +
  geom_col(fill = "#2a9d8f", width = 0.7) +
  geom_text(aes(label = freq), hjust = -0.2, size = 3.5) +
  labs(
    title = "Top 20 most frequent terms in ADHD posts on Reddit",
    subtitle = "Jan 1 to April 20 of: 2018, 2019, and 2020",
    x = "Frequency",
    y = "Term"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.text.y = element_text(size = 11)
  )

# Top 20 terms for each period
top_20_by_period <- freq_by_period %>%
  group_by(period) %>%
  slice_max(n, n = 20) %>%
  ungroup() %>%
  mutate(word = reorder_within(word, n, period))

# plot
ggplot(top_20_by_period, aes(x = n, y = word, fill = period)) +
  geom_col(show.legend = FALSE, width = 0.7) +
  geom_text(aes(label = n), hjust = -0.1, size = 3) +

  facet_wrap(~ period, scales = "free_y") +
  
  scale_y_reordered() +
  scale_fill_manual(values = c("#2a9d8f", "#e9c46a", "#f4a261")) +
  
  labs(
    title = "Top 20 most frequent terms by period",
    x = "Frequency",
    y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold", size = 12)
  )

3.3.3 Word Cloud

# Word cloud

wordcloud(words = freq_df_adhd$word, 
          freq = freq_df_adhd$freq, 
          min.freq = 50,
          max.words = 100,
          random.order = FALSE,
          colors = brewer.pal(8, "Dark2"),
          scale = c(4, 0.5))

title("ADHD posts - most frequent terms")

# set up plotting area
par(mfrow = c(1, 3))
periods <- c("control_2018", "control_2019", "mid_pandemic")

for(p in periods) {
  subset_data <- freq_by_period %>% 
    filter(period == p)
  wordcloud(words = subset_data$word, 
            freq = subset_data$n, 
            min.freq = 50,
            max.words = 80,    
            random.order = FALSE,
            colors = brewer.pal(8, "Dark2"),
            scale = c(3, 0.5)) 
  title(p, line = -2)
}

par(mfrow = c(1, 1))

3.3.4 TF-IDF Matrix

# Term-frequency matrix with TF-IDF
dtm_tfidf_adhd <- DocumentTermMatrix(corpus_adhd_stemmed,
                                     control = list(weighting = function(x) weightTfIdf(x, normalize = FALSE)))

cat("TF-IDF DTM dimensions:\n")
## TF-IDF DTM dimensions:
print(dtm_tfidf_adhd)
## <<DocumentTermMatrix (documents: 28935, terms: 34524)>>
## Non-/sparse entries: 1718786/997233154
## Sparsity           : 100%
## Maximal term length: 477
## Weighting          : term frequency - inverse document frequency (tf-idf)
# reduce sparsity
dtm_tfidf_sparse_adhd <- removeSparseTerms(dtm_tfidf_adhd, sparse = 0.99)

cat("\nReduced TF-IDF DTM dimensions:\n")
## 
## Reduced TF-IDF DTM dimensions:
print(dtm_tfidf_sparse_adhd)
## <<DocumentTermMatrix (documents: 28935, terms: 1019)>>
## Non-/sparse entries: 1378978/28105787
## Sparsity           : 95%
## Maximal term length: 12
## Weighting          : term frequency - inverse document frequency (tf-idf)
dtm_tfidf_adhd_m <- as.matrix(dtm_tfidf_sparse_adhd)
# TF-IDF word frequencies
# count TF-IDF weighted frequencies
tfidf_freq <- sort(colSums(dtm_tfidf_adhd_m), decreasing = TRUE)
tfidf_df_adhd <- data.frame(word = names(tfidf_freq), tfidf = tfidf_freq, row.names = NULL)

cat("Top 20 terms by TF-IDF:\n")
## Top 20 terms by TF-IDF:
print(head(tfidf_df_adhd, 20))
##      word    tfidf
## 1    feel 32683.99
## 2    work 29957.17
## 3    take 28230.61
## 4    time 28131.17
## 5   thing 27890.66
## 6     day 27374.80
## 7   medic 26283.24
## 8    year 24878.83
## 9    help 24817.34
## 10   know 24811.29
## 11   want 24778.27
## 12    tri 24392.35
## 13  start 24197.78
## 14    now 23004.74
## 15   make 22964.58
## 16   even 22924.27
## 17 school 22852.51
## 18  think 22837.56
## 19  peopl 22663.41
## 20    med 22639.03
# Word cloud - TF-IDF

wordcloud(words = tfidf_df_adhd$word, 
          freq = tfidf_df_adhd$tfidf, 
          min.freq = 50,
          max.words = 100,
          random.order = FALSE,
          colors = brewer.pal(8, "Dark2"),
          scale = c(4, 0.5))

junk_words <- c("x200b", "row.two", "timedelta", "self.task", "iâ", "iâ t", "l")

# count words per Period (TF)
period_words <- adhd_combined %>%
  unnest_tokens(word, post) %>%
  anti_join(stop_words, by = "word") %>%
  filter(!word %in% junk_words) %>%
  filter(!str_detect(word, "^www")) %>%
  filter(!str_detect(word, ".com$|.net$")) %>%
  filter(!str_detect(word, "[^a-z']")) %>%
  count(period, word, sort = TRUE) %>%
  ungroup()

# calculate TF-IDF
total_words <- period_words %>% 
  group_by(period) %>% 
  summarize(total = sum(n))

period_words <- left_join(period_words, total_words)

period_tfidf <- period_words %>%
  bind_tf_idf(word, period, n)

# the most distinctive words for 2020 (mid-pandemic), which appear in 2020 but not in 2018 or 2019
cat("Most distinctive words for mid-pandemic (high TF-IDF):\n")
## Most distinctive words for mid-pandemic (high TF-IDF):
period_tfidf %>%
  filter(period == "mid_pandemic") %>%
  arrange(desc(tf_idf)) %>%
  select(period, word, n, tf_idf) %>%
  head(15) %>%
  print()
##          period        word   n       tf_idf
## 1  mid_pandemic  quarantine 301 4.892618e-04
## 2  mid_pandemic       covid 210 3.413455e-04
## 3  mid_pandemic    pandemic 127 2.064327e-04
## 4  mid_pandemic coronavirus 116 1.885527e-04
## 5  mid_pandemic      corona  87 1.414146e-04
## 6  mid_pandemic    lockdown  71 4.259339e-05
## 7  mid_pandemic        zoom  57 3.419470e-05
## 8  mid_pandemic    outbreak  19 3.088364e-05
## 9  mid_pandemic         wfh  19 3.088364e-05
## 10 mid_pandemic      racket  15 2.438182e-05
## 11 mid_pandemic  telehealth  15 2.438182e-05
## 12 mid_pandemic     discord  39 2.339637e-05
## 13 mid_pandemic    youtuber  14 2.275637e-05
## 14 mid_pandemic        etip  12 1.950546e-05
## 15 mid_pandemic       terry  11 1.788000e-05
# plot the highest TF-IDF words for each period
period_tfidf %>%
  group_by(period) %>%
  slice_max(tf_idf, n = 10) %>%
  ungroup() %>%
  mutate(word = reorder_within(word, tf_idf, period)) %>%
  ggplot(aes(word, tf_idf, fill = period)) +
  geom_col(show.legend = FALSE) +
  labs(title = "Most distinctive words by period (TF-IDF)",
       x = NULL, y = "TF-IDF") +
  facet_wrap(~period, scales = "free") +
  coord_flip() +
  scale_x_reordered() +
  theme_minimal()

3.4 Clustering

3.4.1 Hierarchical Clustering

# Hierarchical Clustering

# calculate distance matrix
dist_matrix_adhd <- dist(t(dtm_tfidf_adhd_m), method = "euclidean")

# Hierarchical Clustering using Ward's method
hc_adhd <- hclust(dist_matrix_adhd, method = "ward.D2")

# plot dendrogram
plot(hc_adhd, hang = -1, 
     main = "Hierarchical Clustering of ADHD Terms (Ward's Method)",
     xlab = "Terms", 
     sub = "",
     cex = 0.6)

# draw rectangles around 5 clusters
rect.hclust(hc_adhd, k = 5, border = "red")

# Inspect clusters

# cut tree into 5 clusters
clusters_5 <- cutree(hc_adhd, k = 5)

# show terms in each cluster
for (i in 1:5) {
  cat("\n--- Cluster", i, "---\n")
  cluster_terms <- names(clusters_5[clusters_5 == i])
  cat(paste(cluster_terms[1:min(20, length(cluster_terms))], collapse = ", "), "\n")
}
## 
## --- Cluster 1 ---
## ’ll, ’re, activ, adult, advic, affect, ago, alarm, almost, alreadi, anger, angri, anoth, answer, anxious, anymor, anyon, anyway, articl, ask 
## 
## --- Cluster 2 ---
## ’ve, abl, actual, add, also, alway, anxieti, anyth, appoint, around, assign, back, bad, better, book, brain, cant, class, colleg, come 
## 
## --- Cluster 3 ---
## abil, absolut, abus, academ, accept, accomplish, account, achiev, across, act, action, addict, addit, adhder, adjust, admit, advanc, afford, afraid, afternoon 
## 
## --- Cluster 4 ---
## adderal, concerta, day, dose, effect, hour, med, medic, night, ritalin, side, sleep, take, vyvans, work 
## 
## --- Cluster 5 ---
## ampxb

3.4.2 K-Means Clustering

# K-Means Clustering

# set number of clusters
k <- 4

kmeans_adhd <- kmeans(dtm_tfidf_adhd_m, centers = k, nstart = 25)

# summary
cat("K-means clustering with k =", k, "\n")
## K-means clustering with k = 4
cat("Cluster sizes:\n")
## Cluster sizes:
print(table(kmeans_adhd$cluster))
## 
##     1     2     3     4 
## 19887  5426  2683   939
# visualize K-means clusters

fviz_cluster(list(data = dtm_tfidf_adhd_m, cluster = kmeans_adhd$cluster),
             geom = "point",
             ellipse.type = "convex",
             main = paste("K-means clustering of ADHD posts (k =", k, ")"),
             subtitle = "Based on TF-IDF weighted document-term matrix",
             ggtheme = theme_minimal(base_size = 14))

# top words per cluster
cluster_info <- lapply(1:k, function(i) {
  cluster_docs_idx <- which(kmeans_adhd$cluster == i)
  cluster_docs <- dtm_tfidf_adhd_m[cluster_docs_idx, , drop = FALSE]
  word_freq <- sort(colSums(cluster_docs), decreasing = TRUE)
  top_words <- paste(names(word_freq)[1:10], collapse = ", ")
  data.frame(
    Cluster = i,
    Number_of_documents = length(cluster_docs_idx),
    Top_10_words = top_words,
    stringsAsFactors = FALSE
  )
})

cluster_info_df <- do.call(rbind, cluster_info)

datatable(cluster_info_df,
          caption = "K-means clusters",
          rownames = FALSE,
          options = list(pageLength = 10, dom = 't'))
# word clouds for each cluster

par(mfrow = c(2, 2), mar = c(2, 2, 3, 2))

for (i in 1:k) {
  cluster_docs_idx <- which(kmeans_adhd$cluster == i)
  
  cluster_docs <- dtm_tfidf_adhd_m[cluster_docs_idx, , drop = FALSE]
  word_freq <- colSums(cluster_docs)
  
  wordcloud(names(word_freq), freq = word_freq, 
            max.words = 30, 
            colors = brewer.pal(8, "Dark2"),
            scale = c(3, 0.5))
  title(paste("Cluster", i))
}

par(mfrow = c(1, 1))

3.4.3 Cosine similarity between periods

# centroid vectors for each period - average discussion in each time period

is_2018 <- adhd_combined$period == "control_2018"
is_2019 <- adhd_combined$period == "control_2019"
is_2020 <- adhd_combined$period == "mid_pandemic"

# calculate column means (centroids) for each period
vec_2018_adhd <- colMeans(dtm_tfidf_adhd_m[is_2018, ])
vec_2019_adhd <- colMeans(dtm_tfidf_adhd_m[is_2019, ])
vec_2020_adhd <- colMeans(dtm_tfidf_adhd_m[is_2020, ])

# calculate pairwise cosine similarities
sim_2018_2019 <- cosine(vec_2018_adhd, vec_2019_adhd)
sim_2018_2020 <- cosine(vec_2018_adhd, vec_2020_adhd)
sim_2019_2020 <- cosine(vec_2019_adhd, vec_2020_adhd)

cat("Cosine similarity between periods:\n")
## Cosine similarity between periods:
cat("2018 vs 2019:", round(sim_2018_2019, 4), "\n")
## 2018 vs 2019: 0.9891
cat("2018 vs mid-pandemic 2020:", round(sim_2018_2020, 4), "\n")
## 2018 vs mid-pandemic 2020: 0.9918
cat("2019 vs mid-pandemic 2020:", round(sim_2019_2020, 4), "\n")
## 2019 vs mid-pandemic 2020: 0.9923
cat("Values close to 1 indicate similar discussion patterns across periods\n")
## Values close to 1 indicate similar discussion patterns across periods
cat("1.0 = identical discussions\n")
## 1.0 = identical discussions
cat("0.0 = completely different discussions\n")
## 0.0 = completely different discussions

3.5 Topic Modeling (LDA)

3.5.1 Prepare data for LDA

# 4. VECTORIZATION (INTEGER WORD COUNTS)

# LDA requires raw counts, not TF-IDF
# Remove documents with zero terms (required for LDA)
row_totals <- apply(dtm_sparse_adhd, 1, sum)
dtm_clean_adhd <- dtm_sparse_adhd[row_totals > 0, ]

cat("DTM for LDA:\n")
## DTM for LDA:
print(dtm_clean_adhd)
## <<DocumentTermMatrix (documents: 28931, terms: 1019)>>
## Non-/sparse entries: 1378978/28101711
## Sparsity           : 95%
## Maximal term length: 12
## Weighting          : term frequency (tf)

3.5.2 LDA with k=4 topics

# APPLY LDA TOPIC MODELING

# set number of topics k = 4
k <- 4

set.seed(1234)
lda_model_adhd <- LDA(dtm_clean_adhd, k = k, method = "Gibbs", 
                      control = list(seed = 1234, iter = 1000, burnin = 500))

cat("LDA model fitted with k =", k, "topics\n")
## LDA model fitted with k = 4 topics
# Show top keywords per topic k = 4

# extract topic-word probabilities
lda_topics_adhd <- tidy(lda_model_adhd, matrix = "beta")

# prepare top terms per topic
top_terms_adhd <- lda_topics_adhd %>%
  group_by(topic) %>%
  slice_max(beta, n = 10) %>%
  ungroup() %>%
  arrange(topic, desc(beta))

# visualize top words per topic
ggplot(top_terms_adhd, aes(reorder_within(term, beta, topic), beta, fill = as.factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", labeller = label_both) +
  coord_flip() +
  scale_x_reordered() +
  scale_fill_brewer(palette = "Dark2") +
  theme_minimal(base_size = 12) +
  labs(title = "Top 10 words per topic (LDA with k = 4)",
       subtitle = "ADHD posts",
       x = NULL, 
       y = "Word importance (beta)") +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold", size = 11)
  )

3.5.3 LDA with k=6 topics

# set the number of topics k = 6 

k <- 6

set.seed(1234)
lda_model_adhd_6 <- LDA(dtm_clean_adhd, k = k, method = "Gibbs",
                        control = list(seed = 1234, iter = 1000, burnin = 500))

cat("LDA model fitted with k =", k, "topics\n")
## LDA model fitted with k = 6 topics
# show top keywords per topic k = 6

# extract topic-word probabilities
lda_topics_adhd_6 <- tidy(lda_model_adhd_6, matrix = "beta")

# prepare top terms per topic
top_terms_adhd_6 <- lda_topics_adhd_6 %>%
  group_by(topic) %>%
  slice_max(beta, n = 10) %>%
  ungroup() %>%
  arrange(topic, desc(beta))

# visualize
ggplot(top_terms_adhd_6, aes(reorder_within(term, beta, topic), beta, fill = as.factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", labeller = label_both, ncol = 3) +
  coord_flip() +
  scale_x_reordered() +
  scale_fill_brewer(palette = "Dark2") +
  theme_minimal(base_size = 11) +
  labs(title = "Top 10 words per topic (LDA with k = 6)",
       subtitle = "ADHD posts across all periods",
       x = NULL, 
       y = "Word importance (beta)") +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold", size = 10)
  )

3.5.4 LDA with k=8 topics

# set the number of topics k = 8 

k <- 8

set.seed(1234)
lda_model_adhd_8 <- LDA(dtm_clean_adhd, k = k, method = "Gibbs",
                        control = list(seed = 1234, iter = 1000, burnin = 500))

cat("LDA model fitted with k =", k, "topics\n")
## LDA model fitted with k = 8 topics
# show top keywords per topic k = 8

# extract topic-word probabilities
lda_topics_adhd_8 <- tidy(lda_model_adhd_8, matrix = "beta")

# prepare top terms per topic
top_terms_adhd_8 <- lda_topics_adhd_8 %>%
  group_by(topic) %>%
  slice_max(beta, n = 10) %>%
  ungroup() %>%
  arrange(topic, desc(beta))

# visualize
ggplot(top_terms_adhd_8, aes(reorder_within(term, beta, topic), beta, fill = as.factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", labeller = label_both, ncol = 3) +
  coord_flip() +
  scale_x_reordered() +
  scale_fill_brewer(palette = "Dark2") +
  theme_minimal(base_size = 11) +
  labs(title = "Top 10 words per topic (LDA with k = 8)",
       subtitle = "ADHD posts across all periods",
       x = NULL, 
       y = "Word importance (beta)") +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold", size = 10)
  )

analyze_lda_k <- function(lda_model, k_val) {
  
  adhd_topics <- tidy(lda_model, matrix = "beta")
  
  top_terms <- adhd_topics %>%
    group_by(topic) %>%
    slice_max(beta, n = 10) %>%
    ungroup() %>%
    arrange(topic, -beta)
  
  doc_topics <- posterior(lda_model)$topics
  dominant_topic <- apply(doc_topics, 1, which.max)
  total_docs <- nrow(doc_topics)
  
  # summary table
  lda_summary_df <- top_terms %>%
    group_by(topic) %>%
    summarise(
      Num_docs = sum(dominant_topic == topic),
      Percent_of_docs = paste0(round(100 * sum(dominant_topic == topic) / total_docs, 1), "%"),
      Top_10_Words = paste(term[1:10], collapse = ", "),
      .groups = "drop"
    )
  
  cat("\nSummary for k =", k_val, "\n")
  print(
    datatable(
      lda_summary_df,
      caption = paste("LDA Summary for k =", k_val),
      options = list(dom = 't', paging = FALSE, autoWidth = TRUE),
      rownames = FALSE,
      class = "compact stripe hover"
    )
  )
  
  # get dominant topic counts for plotting
  lda_gamma <- tidy(lda_model, matrix = "gamma") %>%
    group_by(document) %>%
    slice_max(gamma, n = 1) %>%
    ungroup()
  
  p <- ggplot(lda_gamma, aes(x = factor(topic))) +
    geom_bar(fill = "#2a9d8f", width = 0.6, alpha = 0.8) +
    geom_text(stat = 'count', aes(label = after_stat(count)), vjust = -0.5, size = 4) +
    labs(
      title = paste("Number of posts per dominant topic (k =", k_val, ")"),
      x = "Topic",
      y = "Count of posts"
    ) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))
  
  print(p)
}

analyze_lda_k(lda_model_adhd, 4)
## 
## Summary for k = 4

analyze_lda_k(lda_model_adhd_6, 6)
## 
## Summary for k = 6

analyze_lda_k(lda_model_adhd_8, 8)
## 
## Summary for k = 8

3.5.5 Model Selection - Perplexity

# 6. EXPLORE TOPICS

# TEST MULTIPLE k VALUES
# set range of k
k_values <- seq(2, 10, 2)
cat("Running LDA for k =", k, "...\n")
## Running LDA for k = 8 ...
# initialize vector of results
perplexities_adhd <- numeric(length(k_values))

# calculate perplexity for each model
for (i in seq_along(k_values)) {
  k <- k_values[i]
  cat("Running LDA for k =", k, "...\n")
  
  # fit the model
  model <- LDA(
    dtm_clean_adhd,
    k = k,
    method = "Gibbs",
    control = list(seed = 1234, iter = 1000, burnin = 500, thin = 100)
  )
  
  # calculate perplexity
  ll <- as.numeric(logLik(model))
  n_tokens <- sum(dtm_clean_adhd)
  perplexities_adhd[i] <- exp(-ll / n_tokens)
}
## Running LDA for k = 2 ...
## Running LDA for k = 4 ...
## Running LDA for k = 6 ...
## Running LDA for k = 8 ...
## Running LDA for k = 10 ...
# create data frame for visualization
perplexity_df_adhd <- data.frame(
  k = k_values,
  perplexity = perplexities_adhd
)
# visualization
ggplot(perplexity_df_adhd, aes(x = k, y = perplexity)) +
  geom_line(color = "#2a9d8f", linewidth = 1) +
  geom_point(color = "#2a9d8f", size = 3) +
  geom_text(aes(label = k), vjust = -1, size = 3.5) +
  theme_minimal(base_size = 12) +
  labs(
    title = "LDA model fit by number of topics",
    subtitle = "Perplexity diagnostic (lower is better)",
    x = "Number of topics (k)",
    y = "Perplexity"
  ) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank()
  )

Number of topics between 6 and 8 appear to provide a good balance between model complexity and interpretability based on perplexity scores.

Suggested topics for k=8
Topic Suggested label
Topic 1 Social relationships and communication
Topic 2 Medication management
Topic 3 Academic/school life
Topic 4 Employment and career
Topic 5 Diagnosis and medical access
Topic 6 Cognitive symptoms (focus/memory)
Topic 7 Time management and routine
Topic 8 Emotional regulation 
# extract topic probabilities
topics_gamma <- tidy(lda_model_adhd_8, matrix = "gamma")

topics_gamma <- topics_gamma %>%
  mutate(document_id = as.integer(document)) %>%
  mutate(period = adhd_combined$period[document_id])

topic_labels <- tibble(
  topic = 1:8,
  topic_name = c("1. Social/relationships", 
                 "2. Medication management", 
                 "3. School/academic", 
                 "4. Employment/jobs", 
                 "5. Diagnosis process", 
                 "6. Cognitive symptoms", 
                 "7. Routine/sleep", 
                 "8. Emotional regulation")
)

plot_data <- topics_gamma %>%
  left_join(topic_labels, by = "topic") %>%
  mutate(period = factor(period, levels = c("control_2018", "control_2019", "mid_pandemic"))) %>%
  mutate(topic_name = factor(topic_name, levels = topic_labels$topic_name))

plot_data %>%
  group_by(period, topic_name) %>%
  summarize(mean_gamma = mean(gamma), .groups = "drop") %>%
  ggplot(aes(x = topic_name, y = mean_gamma, fill = period)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  
  scale_x_discrete(limits = rev(levels(plot_data$topic_name))) +
  
  labs(title = "Shift of ADHD topics",
       subtitle = "Comparison of topic prevalence (LDA, k=8)",
       x = NULL,
       y = "Topic prevalence (mean gamma)") +
  
  coord_flip() + 
  theme_minimal() +
  
  scale_fill_manual(values = c("control_2018" = "#66c2a5", 
                               "control_2019" = "#fc8d62", 
                               "mid_pandemic" = "#8da0cb")) +
  theme(legend.position = "bottom")

3.6 Sentiment Analysis

3.6.1 Tokenization for Sentiment Analysis

# For dictionary-based sentiment analysis:
# - no corpus object needed
# - no DTM needed  
# - no stemming (must keep full words for lexicon matching)
# - minimal cleaning - remove stopwords

# custom stopwords
custom_stop_words_adhd <- stop_words

# domain-specific words to ignore
custom_stop_words_adhd <- bind_rows(
  tibble(word = c("adhd", "im", "ive", "dont", "thats"), 
         lexicon = "custom"),
  custom_stop_words_adhd
)
# TOKENIZATION

# Tokenize using tidytext
# unnest_tokens():
# - lowercases text
# - removes punctuation
# - standardizes whitespace
# - extracts one-token-per-row

tidy_posts_adhd <- adhd_combined %>%
  unnest_tokens(output = word, input = post)

cat("Number of tokens before stopword removal:", nrow(tidy_posts_adhd), "\n")
## Number of tokens before stopword removal: 5071266
# remove stopwords
tidy_posts_adhd <- tidy_posts_adhd %>%
  anti_join(custom_stop_words_adhd, by = "word")

cat("Number of tokens after stopword removal:", nrow(tidy_posts_adhd), "\n")
## Number of tokens after stopword removal: 1610263
cat("Tokens removed:", nrow(adhd_combined %>% unnest_tokens(word, post)) - nrow(tidy_posts_adhd), "\n")
## Tokens removed: 3461003

3.6.2 Sentiment lexicons

cat("1. Bing: binary positive/negative classification\n")
## 1. Bing: binary positive/negative classification
bing_summary <- get_sentiments("bing") %>% count(sentiment)
print(bing_summary)
## # A tibble: 2 × 2
##   sentiment     n
##   <chr>     <int>
## 1 negative   4781
## 2 positive   2005
cat("\n2. AFINN: numeric scale from -5 (negative) to +5 (positive)\n")
## 
## 2. AFINN: numeric scale from -5 (negative) to +5 (positive)
afinn_summary <- get_sentiments("afinn") %>%
  summarize(min_value = min(value), max_value = max(value), n_words = n())
print(afinn_summary)
## # A tibble: 1 × 3
##   min_value max_value n_words
##       <dbl>     <dbl>   <int>
## 1        -5         5    2477
cat("\n3. NRC: emotions + positive/negative\n")
## 
## 3. NRC: emotions + positive/negative
nrc_summary <- get_sentiments("nrc") %>% count(sentiment) %>% arrange(desc(n))
print(nrc_summary)
## # A tibble: 10 × 2
##    sentiment        n
##    <chr>        <int>
##  1 negative      3316
##  2 positive      2308
##  3 fear          1474
##  4 anger         1245
##  5 trust         1230
##  6 sadness       1187
##  7 disgust       1056
##  8 anticipation   837
##  9 joy            687
## 10 surprise       532
# Visualize lexicon distributions

# Bing sentiment counts
p1 <- get_sentiments("bing") %>%
  count(sentiment) %>%
  mutate(sentiment = fct_reorder(sentiment, n)) %>%
  ggplot(aes(x = sentiment, y = n, fill = sentiment)) +
  geom_col(width = 0.5, show.legend = FALSE) +
  geom_text(aes(label = n), vjust = -0.5, size = 3.5) +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = "Bing lexicon", x = NULL, y = "Word count") +
  theme_minimal(base_size = 11)

# NRC sentiment counts
p2 <- get_sentiments("nrc") %>%
  count(sentiment) %>%
  mutate(sentiment = fct_reorder(sentiment, n)) %>%
  ggplot(aes(x = sentiment, y = n, fill = sentiment)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  coord_flip() +
  labs(title = "NRC lexicon", x = NULL, y = "Word count") +
  theme_minimal(base_size = 11)

# AFINN value distribution
p3 <- get_sentiments("afinn") %>%
  count(value) %>%
  mutate(value = fct_reorder(as.factor(value), value, .desc = FALSE)) %>%
  ggplot(aes(x = value, y = n, fill = value)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  coord_flip() +
  labs(title = "AFINN lexicon", x = "Sentiment value", y = "Word count") +
  theme_minimal(base_size = 11)

grid.arrange(p1, p2, p3, ncol = 3, 
             top = "Sentiment lexicon distributions")

3.6.3 Bing

# Bing dictionary analysis
bing_dict <- get_sentiments("bing") %>% 
  distinct(word, .keep_all = TRUE)

bing_results_adhd <- tidy_posts_adhd %>%
  inner_join(bing_dict, by = "word")

cat("Number of sentiment-bearing words found:", nrow(bing_results_adhd), "\n")
## Number of sentiment-bearing words found: 251040
bing_word_counts_adhd <- bing_results_adhd %>%
  count(word, sentiment, sort = TRUE)

cat("\nTop 10 positive words:\n")
## 
## Top 10 positive words:
print(bing_word_counts_adhd %>% filter(sentiment == "positive") %>% head(10))
##          word sentiment    n
## 1      pretty  positive 3254
## 2        love  positive 2593
## 3       super  positive 1734
## 4      helped  positive 1715
## 5       happy  positive 1519
## 6        fine  positive 1176
## 7  productive  positive 1176
## 8        fast  positive 1124
## 9        easy  positive 1031
## 10      clean  positive  996
cat("\nTop 10 negative words:\n")
## 
## Top 10 negative words:
print(bing_word_counts_adhd %>% filter(sentiment == "negative") %>% head(10))
##          word sentiment    n
## 1     anxiety  negative 5539
## 2        hard  negative 5262
## 3         bad  negative 4039
## 4    symptoms  negative 3778
## 5  depression  negative 3189
## 6      issues  negative 2908
## 7        shit  negative 2212
## 8       worse  negative 2039
## 9    struggle  negative 2028
## 10      tired  negative 2023
# visualize top sentiment words

bing_top_words_adhd <- bing_word_counts_adhd %>%
  group_by(sentiment) %>%
  slice_max(n, n = 10) %>%
  ungroup() %>%
  mutate(word = reorder_within(word, n, sentiment))

ggplot(bing_top_words_adhd, aes(x = n, y = word, fill = sentiment)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ sentiment, scales = "free_y") +
  scale_y_reordered() +
  scale_fill_manual(values = c("negative" = "#d7191c", "positive" = "#2b83ba")) +
  labs(
    title = "Top 10 sentiment words in ADHD posts",
    subtitle = "Bing dictionary",
    x = "Frequency",
    y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold", size = 12)
  )

3.6.3.1 Sentiment by period

sentiment_by_period_adhd <- bing_results_adhd %>%
  count(period, sentiment) %>%
  pivot_wider(
    names_from = sentiment,
    values_from = n,
    values_fill = 0
  ) %>%
  mutate(
    total = positive + negative,
    net_sentiment = positive - negative,
    sentiment_ratio = positive / (positive + negative)
  )

print(sentiment_by_period_adhd)
## # A tibble: 3 × 6
##   period       negative positive  total net_sentiment sentiment_ratio
##   <chr>           <int>    <int>  <int>         <int>           <dbl>
## 1 control_2018    34212    14891  49103        -19321           0.303
## 2 control_2019    62516    27028  89544        -35488           0.302
## 3 mid_pandemic    79038    33355 112393        -45683           0.297
sentiment_plot_data_adhd <- sentiment_by_period_adhd %>%
  pivot_longer(
    cols = c(positive, negative),
    names_to = "sentiment",
    values_to = "count"
  ) %>%
  mutate(
    period = factor(period, levels = c("control_2018", "control_2019", "mid_pandemic")),
    sentiment = factor(sentiment, levels = c("positive", "negative"))
  )

# stacked bar chart
ggplot(sentiment_plot_data_adhd, aes(x = period, y = count, fill = sentiment)) +
  geom_col(position = "stack", width = 0.6) +
  scale_fill_manual(
    values = c("positive" = "#2b83ba", "negative" = "#d7191c"),
    labels = c("Positive", "Negative")
  ) +
  labs(
    title = "Sentiment distribution",
    subtitle = "ADHD posts using Bing",
    x = "Period",
    y = "Count of words",
    fill = "Sentiment"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "right"
  )

# net sentiment by period

ggplot(sentiment_by_period_adhd, aes(x = period, y = net_sentiment)) +
  geom_col(aes(fill = net_sentiment > 0), width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = net_sentiment), vjust = -0.5, size = 4) +
  scale_fill_manual(values = c("TRUE" = "#2b83ba", "FALSE" = "#d7191c")) +
  labs(
    title = "Net sentiment score",
    subtitle = "Positive minus negative (Bing)",
    x = "Period",
    y = "Net sentiment score"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14)
  )

3.6.4 AFINN

# AFINN dictionary analysis

afinn_results_adhd <- tidy_posts_adhd %>%
  inner_join(get_sentiments("afinn"), by = "word", relationship = "many-to-one")

afinn_by_period_adhd <- afinn_results_adhd %>%
  group_by(period) %>%
  summarise(
    mean_sentiment = mean(value),
    median_sentiment = median(value),
    total_sentiment = sum(value),
    n_words = n(),
    .groups = "drop"
  ) %>%
  mutate(period = factor(period, levels = c("control_2018", "control_2019", "mid_pandemic")))

print(afinn_by_period_adhd)
## # A tibble: 3 × 5
##   period       mean_sentiment median_sentiment total_sentiment n_words
##   <fct>                 <dbl>            <dbl>           <dbl>   <int>
## 1 control_2018         -0.632               -1          -28439   44965
## 2 control_2019         -0.592               -1          -48466   81832
## 3 mid_pandemic         -0.617               -1          -63911  103635
# visualize AFINN sentiment

ggplot(afinn_by_period_adhd, aes(x = period, y = mean_sentiment)) +
  geom_col(aes(fill = mean_sentiment > 0), width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = round(mean_sentiment, 2)), vjust = -0.5, size = 4) +
  scale_fill_manual(values = c("TRUE" = "#2b83ba", "FALSE" = "#d7191c")) +
  labs(
    title = "Average sentiment score period",
    subtitle = "AFINN lexicon (scale: -5 to +5)",
    x = "Period",
    y = "Mean sentiment score"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14)
  )

3.6.5 NRC emotion analysis

nrc_results_adhd <- tidy_posts_adhd %>%
  inner_join(get_sentiments("nrc"), by = "word", relationship = "many-to-many")

# count emotions by period
nrc_by_period_adhd <- nrc_results_adhd %>%
  count(period, sentiment) %>%
  mutate(period = factor(period, levels = c("control_2018", "control_2019", "mid_pandemic")))

cat("Emotion counts by period:\n")
## Emotion counts by period:
print(nrc_by_period_adhd)
##          period    sentiment     n
## 1  control_2018        anger 12696
## 2  control_2018 anticipation 20099
## 3  control_2018      disgust  8874
## 4  control_2018         fear 15349
## 5  control_2018          joy 11923
## 6  control_2018     negative 29709
## 7  control_2018     positive 34615
## 8  control_2018      sadness 15569
## 9  control_2018     surprise  7002
## 10 control_2018        trust 21354
## 11 control_2019        anger 22624
## 12 control_2019 anticipation 37240
## 13 control_2019      disgust 15846
## 14 control_2019         fear 27488
## 15 control_2019          joy 22196
## 16 control_2019     negative 54079
## 17 control_2019     positive 63641
## 18 control_2019      sadness 28204
## 19 control_2019     surprise 13062
## 20 control_2019        trust 39785
## 21 mid_pandemic        anger 28611
## 22 mid_pandemic anticipation 45663
## 23 mid_pandemic      disgust 20303
## 24 mid_pandemic         fear 35405
## 25 mid_pandemic          joy 27575
## 26 mid_pandemic     negative 68657
## 27 mid_pandemic     positive 78296
## 28 mid_pandemic      sadness 35655
## 29 mid_pandemic     surprise 16298
## 30 mid_pandemic        trust 48651
# visualize NRC emotions

ggplot(nrc_by_period_adhd, aes(x = sentiment, y = n, fill = period)) +
  geom_col(position = "dodge", width = 0.7) +
  coord_flip() +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    title = "Emotion distribution",
    subtitle = "NRC emotion lexicon - ADHD posts",
    x = "Emotion",
    y = "Frequency",
    fill = "Period"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "right"
  )

3.6.6 Sentiment word clouds

# WORDCLOUD OF SENTIMENT WORDS

par(mfrow = c(1, 2))

# Positive words
positive_words_adhd <- bing_word_counts_adhd %>%
  filter(sentiment == "positive")

wordcloud(words = positive_words_adhd$word, 
          freq = positive_words_adhd$n,
          min.freq = 5,
          max.words = 100,
          random.order = FALSE,
          colors = brewer.pal(8, "Blues")[4:8],
          scale = c(4, 0.5))
title("Positive sentiment words")

# Negative words
negative_words_adhd <- bing_word_counts_adhd %>%
  filter(sentiment == "negative")

wordcloud(words = negative_words_adhd$word, 
          freq = negative_words_adhd$n,
          min.freq = 5,
          max.words = 100,
          random.order = FALSE,
          colors = brewer.pal(8, "Reds")[4:8],
          scale = c(4, 0.5))
title("Negative sentiment words")

par(mfrow = c(1, 1))

3.6.7 Comparison period

clean_comparison_data <- adhd_combined %>%
  unnest_tokens(word, post) %>%
  mutate(word = str_remove_all(word, "[[:punct:]]")) %>%
  anti_join(stop_words, by = "word") %>%
  filter(!word %in% c("adhd", "im", "dont", "can", "just", "like", "get", "amp", "ve", "ive", "cant", "wont", "didnt", "x200b", "ve", "ill", "cant", "adhd", "x200b", "amp", "im", "dont", "can", "ive", "its", "ive")) %>%
  filter(!str_detect(word, "^[0-9]+$")) %>%
  count(period, word) %>%
  acast(word ~ period, value.var = "n", fill = 0)

comparison.cloud(clean_comparison_data,
                 max.words = 100,
                 scale = c(3.5, 0.5),
                 colors = c("#1b9e77", "#d95f02", "#7570b3"),
                 title.size = 1,
                 random.order = FALSE,
                 match.colors = TRUE)

title("Comparison word cloud by period", line = 1)

3.7 Text Networks and Bigrams

3.7.1 Build Bigrams

# create bigrams (word pairs)
bigrams_adhd <- adhd_combined %>%
  unnest_tokens(output = bigram, 
                input = post, 
                token = "ngrams", 
                n = 2) %>%
  separate(bigram, into = c("word1", "word2"), sep = " ")

cat("Total bigrams created:", nrow(bigrams_adhd), "\n")
## Total bigrams created: 5042331
# filter bigrams

# custom stopwords (keeping negators)
custom_stop_bigrams <- stop_words

# add domain-specific stopwords
custom_stop_bigrams <- bind_rows(
  tibble(word = c("adhd", "i'm", "x200b", "tl", "dr", "hey", "guys", "ago"), 
         lexicon = "custom"),
  custom_stop_bigrams
)

# filter bigrams: remove stopwords from both positions
bigrams_filtered_adhd <- bigrams_adhd %>%
  filter(!word1 %in% custom_stop_bigrams$word,
         !word2 %in% custom_stop_bigrams$word,
         !str_detect(word1, "^[0-9]+$"),
         !str_detect(word2, "^[0-9]+$"),
         !is.na(word1),
         !is.na(word2),
         )

cat("Bigrams after filtering:", nrow(bigrams_filtered_adhd), "\n")
## Bigrams after filtering: 383240
# count bigrams

# count most frequent bigrams
bigram_counts_adhd <- bigrams_filtered_adhd %>%
  count(word1, word2, sort = TRUE)

cat("Top 20 bigrams:\n")
## Top 20 bigrams:
print(head(bigram_counts_adhd, 20))
##       word1      word2   n
## 1    mental     health 816
## 2  adderall         xr 808
## 3      hard       time 759
## 4   started     taking 652
## 5  recently  diagnosed 557
## 6       pay  attention 553
## 7      fall     asleep 418
## 8     video      games 403
## 9    entire       life 401
## 10   taking   adderall 398
## 11  stopped     taking 387
## 12    heart       rate 334
## 13   paying  attention 331
## 14   social    anxiety 331
## 15    makes      sense 290
## 16 adderall         ir 285
## 17     fast    forward 284
## 18   taking medication 282
## 19   taking       meds 274
## 20    brain        fog 272
# visualize most frequent bigrams

bigrams_top20_adhd <- bigram_counts_adhd %>%
  slice_max(n, n = 20) %>%
  mutate(
    bigram = paste(word1, word2),
    bigram = fct_reorder(bigram, n)
  )

ggplot(bigrams_top20_adhd, aes(x = n, y = bigram)) +
  geom_col(fill = "#2a9d8f", width = 0.7) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  labs(
    title = "Most frequent bigrams in ADHD pstss",
    x = "Frequency",
    y = "Bigram"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14)
  )

3.7.2 Bigram Network visualization

# create bigram network

# filter for bigrams
bigram_graph_adhd <- bigram_counts_adhd %>%
  filter(n >= 100) %>%  # appearing at least 100 times
  graph_from_data_frame()

# create network visualization
ggraph(bigram_graph_adhd, layout = "fr") +
  
  # EDGES
  geom_edge_link(aes(edge_alpha = n, edge_width = n), 
                 edge_colour = "cyan4",
                 show.legend = FALSE) +
  
  # NODES
  geom_node_point(color = "darkblue", size = 4) +
  
  # LABELS
  geom_node_text(aes(label = name), 
                 vjust = 1.5, hjust = 0.5, 
                 size = 3.5,
                 repel = TRUE) +
  
  scale_edge_width(range = c(0.5, 3)) +
  scale_edge_alpha(range = c(0.3, 0.8)) +
  
  theme_void() +
  labs(
    title = "Bigram Network of ADHD posts",
    subtitle = "Minimum 20 occurrences"
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5, margin = ggplot2::margin(b = 20))
  )

3.7.3 Network by period

# start from the original text not to lose the year labels
bigrams_with_period <- adhd_combined %>%
  unnest_tokens(bigram, post, token = "ngrams", n = 2) %>%
  separate(bigram, c("word1", "word2"), sep = " ") %>%
  
  # remove stop words and numbers
  filter(!word1 %in% stop_words$word) %>%
  filter(!word2 %in% stop_words$word) %>%
  filter(!word1 %in% c("adhd", "amp", "x200b", "people", "time", "don't", "i'm")) %>%
  filter(!word2 %in% c("adhd", "amp", "x200b", "people", "time", "don't", "i'm")) %>%
  filter(!str_detect(word1, "^[0-9]+$")) %>%
  filter(!str_detect(word2, "^[0-9]+$")) %>%

  count(period, word1, word2, sort = TRUE)


plot_network <- function(data, period_name, min_n = 15) {
  
  plot_data <- data %>%
    filter(period == period_name) %>%
    filter(n >= min_n) %>%
    select(word1, word2, n)
  
  if(nrow(plot_data) == 0) return(NULL)
  
  # plot
  graph_from_data_frame(plot_data) %>%
    ggraph(layout = "fr") +
    geom_edge_link(aes(edge_alpha = n, edge_width = n), 
                   edge_colour = "cyan4", show.legend = FALSE) +
    geom_node_point(color = "darkblue", size = 3) +
    geom_node_text(aes(label = name), vjust = 1.5, size = 3, repel = TRUE) +
    labs(title = paste("Network:", period_name), 
         subtitle = paste("Connections appearing >=", min_n, "times")) +
    theme_void() + 
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p1 <- plot_network(bigrams_with_period, "control_2018", min_n = 15)
p2 <- plot_network(bigrams_with_period, "control_2019", min_n = 15)
p3 <- plot_network(bigrams_with_period, "mid_pandemic", min_n = 15)

print(p1)

print(p2)

print(p3)

3.7.4 Word Associations Using TDM

# ASSOCIATIONS – FINDING CO-OCCURRING WORDS

# create TDM for association analysis
tdm_adhd <- TermDocumentMatrix(corpus_adhd_stemmed)

# find words associated with key terms
cat("Associations with 'take':\n")
## Associations with 'take':
print(findAssocs(tdm_adhd, "take", 0.3))
## $take
##  day 
## 0.31
cat("\nAssociations with 'effect':\n")
## 
## Associations with 'effect':
print(findAssocs(tdm_adhd, "effect", 0.3))
## $effect
## side 
## 0.62
cat("\nAssociations with 'sleep':\n")
## 
## Associations with 'sleep':
print(findAssocs(tdm_adhd, "sleep", 0.3))
## $sleep
##          “sleep      adenosine—    boardcertifi          drive—     importance— 
##            0.40            0.40            0.40            0.40            0.40 
##    neurologist’ proinflammatori        society’          depriv            wake 
##            0.40            0.40            0.40            0.35            0.33 
##          advent        sleepwak 
##            0.30            0.30
# visualization of associations

target_word <- "take"
cor_limit <- 0.2

# calculate associations
associations_adhd <- findAssocs(tdm_adhd, target_word, corlimit = cor_limit)
assoc_vector_adhd <- associations_adhd[[target_word]]

# if no associations found
if (length(assoc_vector_adhd) > 0) {
  assoc_sorted_adhd <- sort(assoc_vector_adhd, decreasing = TRUE)
  
  assoc_df_adhd <- data.frame(
    word = names(assoc_sorted_adhd),
    score = assoc_sorted_adhd,
    row.names = NULL
  )
  
  # lollipop chart
  ggplot(assoc_df_adhd, aes(x = score, y = reorder(word, score), color = score)) +
    geom_segment(aes(x = 0, xend = score, y = word, yend = word), 
                 linewidth = 1.2) +
    geom_point(size = 4) +
    geom_text(aes(label = round(score, 2)), 
              hjust = -0.3, size = 3.5, color = "black") +
    scale_color_gradient(low = "#a6bddb", high = "#08306b") +
    scale_x_continuous(
      limits = c(0, max(assoc_df_adhd$score) + 0.1),
      expand = expansion(mult = c(0, 0.15))
    ) +
    theme_minimal(base_size = 12) +
    labs(
      title = paste0("Word associations with '", target_word, "'"),
      subtitle = paste0("Pearson correlation threshold r ≥ ", cor_limit),
      x = "Correlation coefficient",
      y = "Associated word",
      color = "Association\nStrength"
    ) +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      legend.position = "right"
    )
} else {
  cat("No associations found for", target_word, "at threshold", cor_limit, "\n")
}

# visualization of associations

target_word <- "effect"
cor_limit <- 0.2

# calculate associations
associations_adhd <- findAssocs(tdm_adhd, target_word, corlimit = cor_limit)
assoc_vector_adhd <- associations_adhd[[target_word]]

# if no associations found
if (length(assoc_vector_adhd) > 0) {
  assoc_sorted_adhd <- sort(assoc_vector_adhd, decreasing = TRUE)
  
  assoc_df_adhd <- data.frame(
    word = names(assoc_sorted_adhd),
    score = assoc_sorted_adhd,
    row.names = NULL
  )
  
  # lollipop chart
  ggplot(assoc_df_adhd, aes(x = score, y = reorder(word, score), color = score)) +
    geom_segment(aes(x = 0, xend = score, y = word, yend = word), 
                 linewidth = 1.2) +
    geom_point(size = 4) +
    geom_text(aes(label = round(score, 2)), 
              hjust = -0.3, size = 3.5, color = "black") +
    scale_color_gradient(low = "#a6bddb", high = "#08306b") +
    scale_x_continuous(
      limits = c(0, max(assoc_df_adhd$score) + 0.1),
      expand = expansion(mult = c(0, 0.15))
    ) +
    theme_minimal(base_size = 12) +
    labs(
      title = paste0("Word associations with '", target_word, "'"),
      subtitle = paste0("Pearson correlation threshold r ≥ ", cor_limit),
      x = "Correlation coefficient",
      y = "Associated word",
      color = "Association\nStrength"
    ) +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      legend.position = "right"
    )
} else {
  cat("No associations found for", target_word, "at threshold", cor_limit, "\n")
}

# visualization of associations

target_word <- "sleep"
cor_limit <- 0.2

# calculate associations
associations_adhd <- findAssocs(tdm_adhd, target_word, corlimit = cor_limit)
assoc_vector_adhd <- associations_adhd[[target_word]]

# if no associations found
if (length(assoc_vector_adhd) > 0) {
  assoc_sorted_adhd <- sort(assoc_vector_adhd, decreasing = TRUE)
  
  assoc_df_adhd <- data.frame(
    word = names(assoc_sorted_adhd),
    score = assoc_sorted_adhd,
    row.names = NULL
  )
  
  # lollipop chart
  ggplot(assoc_df_adhd, aes(x = score, y = reorder(word, score), color = score)) +
    geom_segment(aes(x = 0, xend = score, y = word, yend = word), 
                 linewidth = 1.2) +
    geom_point(size = 4) +
    geom_text(aes(label = round(score, 2)), 
              hjust = -0.3, size = 3.5, color = "black") +
    scale_color_gradient(low = "#a6bddb", high = "#08306b") +
    scale_x_continuous(
      limits = c(0, max(assoc_df_adhd$score) + 0.1),
      expand = expansion(mult = c(0, 0.15))
    ) +
    theme_minimal(base_size = 12) +
    labs(
      title = paste0("Word associations with '", target_word, "'"),
      subtitle = paste0("Pearson correlation threshold r ≥ ", cor_limit),
      x = "Correlation coefficient",
      y = "Associated word",
      color = "Association\nStrength"
    ) +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      legend.position = "right"
    )
} else {
  cat("No associations found for", target_word, "at threshold", cor_limit, "\n")
}

4 Summary

This analysis successfully applied multiple text mining techniques to examine mental health discussions on Reddit across different time periods. Key achievements:

  1. Exploratory Data Analysis - identified the most frequent terms and themes in ADHD subreddit
  2. Clustering - grouped similar terms and posts to reveal underlying patterns in mental health discourse
  3. Topic Modeling - discovered latent topics using LDA, revealing distinct discussion themes
  4. Sentiment Analysis - tracked emotional tone across periods using multiple lexicons (Bing, AFINN, NRC)
  5. Text Networks - visualized word co-occurrence patterns through bigram networks

Examied assumptions:

  1. Control periods appropriately account for seasonal variation - cosine similarity analysis between 2018 and 2019 ADHD discussions showed high similarity (> 0.9), supporting the validity of using these as control periods.
  2. Mid-pandemic period represents a distinct change point - sentiment analysis revealed changes in emotional tone during the mid-pandemic period.
  3. Text preprocessing preserves meaningful content - word frequency analysis and topic modeling produced interpretable results and domain-specific stopword removal (e.g. “adhd”) improved signal-to-noise ratio.
  4. Bigram networks capture contextual relationships - network visualizations revealed meaningful word associations and period comparison of networks showed evolution of discussion patterns.
  5. Dictionary-based sentiment is valid for mental health discourse - multiple lexicons (Bing, AFINN, NRC) showed consistent directional trends, however general-purpose lexicons may miss domain-specific nuances.

Key findings in ADHD discussions:

  1. Clustering - terms grouped into themes related to general, medication, emotional and life impact, routine
  2. Topics - LDA revealed distinct topics including e.g. medication management, daily functioning, and emotional experiences
  3. Sentiment - relatively stable sentiment across periods with slight variations
  4. Networks - co-occurrence of medication-related terms and symptom descriptions

Methodological strengths:

  1. Multi-method approach - combining unsupervised (clustering, topic modeling) and lexicon-based (sentiment) methods
  2. Temporal comparison - using multiple control periods strengthens conclusions about pandemic-related changes
  3. Reproducibility - seed setting and documented preprocessing steps enable replication

Limitations:

  1. Dictionary limitations - general-purpose sentiment lexicons may not capture mental health-specific language nuances
  2. Context loss - bag-of-words approaches lose sentence-level context
  3. Selection bias - analysis limited to specific subreddits and timeframes
  4. Stemming trade-offs - stemming reduces vocabulary but may merge distinct concepts

5 References

  1. Low, D. M., Rumker, L., Torous, J., Cecchi, G., Ghosh, S. S., & Talkar, T. (2020). Natural Language Processing Reveals Vulnerable Mental Health Support Groups and Heightened Health Anxiety on Reddit During COVID-19: Observational Study. Journal of Medical Internet Research, 22(10), e22635. https://doi.org/10.2196/2263
  2. “Text Mining and Social Media Mining” course materials