stm01.Rmd

---
title: "Structural Topic Modeling of the Anger Narratives"
author: "Pooya Razavi"
date: "2022-12-06"
output: 
  html_document:
    toc: true
    theme: united
editor_options: 
  chunk_output_type: console
---



```{r setup, include=FALSE, warning=FALSE}
#load libraries
package_list <- c("dplyr", "tidyr", "ggplot2", "tidytext", "topicmodels", "stm")
lapply(package_list, require, character.only = TRUE)


df <- readxl::read_xlsx("C:/Users/pooya/Dropbox (University of Oregon)/Anger Dissertation/Prototype study analysis/ProcessedData_F21_W22_S22_F22.xlsx")

knitr::opts_chunk$set(echo = TRUE)
```

# Functions

```{r topic-prop-plot-function}
stm_plot <- function(model_name, perc_bar_lim = .3, ylabel = "yes") {
  
  #get the expected proportions for each topic
    exp_proportions <- make.dt(model_name) %>% 
                              summarise_all(mean) %>% 
                                t() %>% as.data.frame() %>% 
                                tibble::rownames_to_column("Topic") %>% 
                                filter(Topic != "docnum")
    #number of topics
    t_number <- nrow(exp_proportions)
    
    #get the top word list
    exp_proportions <- exp_proportions %>% 
                                cbind((labelTopics(model_name, c(1:t_number), n = 7))[["frex"]])
    
    #attach it to the proportions df
    top_words <- paste(exp_proportions$`1`, exp_proportions$`2`,
                        exp_proportions$`3`, exp_proportions$`4`,
                        exp_proportions$`5`, exp_proportions$`6`,
                        exp_proportions$`7`, sep = ", ")
    
    exp_proportions <- exp_proportions %>% 
                                mutate(top_words = top_words)
    
    #ggplot
    exp_proportions %>% 
      ggplot(., aes(x = reorder(Topic, V1), y = V1, label = top_words, fill = Topic)) +
        geom_col(show.legend = FALSE) +
        geom_text(hjust = 0, nudge_y = 0.002, size = 4,
                  family = "serif") +
        coord_flip() +
        scale_y_continuous(expand = c(0,0),
                           limits = c(0, perc_bar_lim),
                           labels = scales::percent_format())  +
      ggthemes::theme_tufte(base_family = "serif", ticks = FALSE) +
      theme(plot.title = element_text(size = 16,
                                      family="serif"),
            plot.subtitle = element_text(size = 13),
        axis.line.x = element_line(color='darkgrey')) +
      labs(x = NULL, y = if_else(ylabel == "yes", "Expected Topic Proportions", ""),
           title = paste0(t_number, " Topics"))
           
  
  
}


```



```{r, report-covariate-function}
#function to report the comparisons in a single table

report_stm_cov_effect <- function(prepped_model, save_as_csv = "no") {
      
      model_results <- data.frame(Topic = NA,
                                  Var = NA,
                                  Estimate = NA,
                                  SE = NA,
                                  t = NA,
                                  p = NA)
      t_number <- max(prepped_model[["topics"]])
      
      for(topic in 1:t_number){
        #print(topic)
        
        topic_summary <- (summary(prepped_model, topics = topic))[["tables"]] %>% 
                              as.data.frame() %>% 
                              tibble::rownames_to_column("Var")
        
        colnames(topic_summary) <- c("Var", "Estimate", "SE", "t", "p")
        Topic <- c(topic, topic)
        topic_summary <- cbind(Topic, topic_summary)
        
        model_results <- rbind(model_results, topic_summary)
      }
      
      if(save_as_csv == "no") {
      model_results[-1,] %>% format(scientific = FALSE) %>% 
          knitr::kable(digits = 3) %>% kableExtra::kable_styling()
      } else{
        model_results[-1,] %>% format(scientific = FALSE) %>%
          write.csv(., save_as_csv)
      }
}
```


```{r, topic-quality-function}
#generates exclusivity-coherence plots using ggplot

topic_quality_plot <- function(model_name){
  n_t <- length(stm::semanticCoherence(model_name, documents = out$documents))
  topic_q_df <- data.frame(topic = paste0("T", 1:n_t),
                      coher = stm::semanticCoherence(model_name, documents = out$documents),
                      exclus = stm::exclusivity(model_name))
  topic_q_df <- (topic_q_df %>% 
                  rbind(c("Average", mean(.$coher), mean(.$exclus)))) %>% 
                  mutate_at(c('coher', 'exclus'), as.numeric) %>% 
                  mutate(color = if_else(topic == "Average", "red", "darkblue"))
  
  ggplot(topic_q_df, aes(x = coher, y = exclus)) +
    geom_point(color = topic_q_df$color) +
    geom_text(color = topic_q_df$color, label = topic_q_df$topic, nudge_y = 0.1, size = 3, family="serif") +
    labs(title = paste(n_t, "Topics"), x = "Coherence", y = "Exclusivity") +
    xlim(-210, -40) +
    ylim(8.5, 10.25) +
    theme_bw() +
    theme(text=element_text(size = 12, family="serif"))
}
```


# Pre-processing

```{r ingest}
 #apply the preregistered data exclusion
    #assigning values to factor levels
      df$NarrativeWritten <- as.factor(df$NarrativeWritten)
      df$NarrativeRelevant <- as.factor(df$NarrativeRelevant)
      df$Condition <- as.factor(df$Condition)
      
      levels(df$NarrativeWritten) <- c("No", "Yes")
      levels(df$NarrativeRelevant) <- c("No", "Yes", NA, NA) 
      levels(df$Condition) <- c("justified", "nonjustified", NA)
    
    #drop cases following preregistration
      df1 <- df %>% 
        filter(NarrativeWritten != "No") %>% 
        filter(NarrativeRelevant != "No") %>% 
        filter(!is.na(Condition))
      
 
  #keep the relevant variables and drop the rest
      df_stm <- df1 %>% 
        mutate(all_narratives = dplyr::coalesce(right_narrative, nonright_narrative)) %>% 
        select(ResponseId, Condition, all_narratives)
      

```


```{r preprocessing-5threshold}
#stemming, dropping punctuation, numbers, and stop words
  processed <- textProcessor(df_stm$all_narratives, metadata = df_stm) 

#evaluate how many words and documents would be removed from the data set at each word threshold
  plotRemoved(processed$documents, lower.thresh = seq(1, 20, by = 1))
  
#saving the different parts of the processed data
  out <- prepDocuments(processed$documents, processed$vocab, processed$meta, lower.thresh = 5) #including words that appear in at least 5 documents
  docs <- out$documents
  vocab <- out$vocab
  meta <- out$meta

```

# Determining the Optimal K

```{r eval=FALSE}
#Method: using the searchK algorithm
  set.seed(110)
  storage <- searchK(out$documents, out$vocab, K = c(2:40), prevalence = ~Condition, data = meta, seed = 110)
    #this warning was produced: "K=2 is equivalent to a unidimensional scaling model which you may prefer."
  
  plot(storage) #looks like 9, 11, and 14 topics are good candidates (when lowerthreshold = 5)

  #save the results of searchK with frequency threshold = 5 and K = 2:40
    searchk_results <- storage[["results"]] %>% as.data.frame() %>% apply(., 2, as.numeric)
    #write.csv(searchk_results, "searchk_threshold5_topics2to40_results.csv")

```

## Visualize SearchK Results

```{r}
searchk_df <- read.csv("searchk_threshold5_topics2to40_results.csv") %>% 
                      select(-X) #drop the first column

#Exclusivity
  exclus_plot <- ggplot(searchk_df, aes(x = K, y = exclus)) +
                  geom_point() +
                  geom_line() +
                  labs(title = "Topic Exclusivity",
                       x = "Topic",
                       y = "Exclusivity") +
                  theme_bw() +
                  theme(
                  plot.background = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.grid.major = element_blank(),
                  text=element_text(size = 12, family="serif"))

#Semantic Coherence
  coher_plot <- ggplot(searchk_df, aes(x = K, y = semcoh)) +
                  geom_point(color = "darkblue") +
                  geom_line(color = "darkblue") +
                  labs(title = "Semantic Coherence",
                       x = "Topic",
                       y = "Coherence") +
                  theme_bw() +
                  theme(
                  plot.background = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.grid.major = element_blank(),
                  text=element_text(size = 12, family="serif"))

#Residual
  resid_plot <- ggplot(searchk_df, aes(x = K, y = residual)) +
                  geom_point(color = "darkred") +
                  geom_line(color = "darkred") +
                  labs(title = "Residuals",
                       x = "Topic",
                       y = "Residual") +
                  theme_bw() +
                  theme(
                  plot.background = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.grid.major = element_blank(),
                  text=element_text(size = 12, family="serif"))

#Held-out likelihood
  heldout_plot <- ggplot(searchk_df, aes(x = K, y = heldout)) +
                    geom_point(color = "darkgreen") +
                    geom_line(color = "darkgreen") +
                    labs(title = "Held-out Likelihood",
                         x = "Topic",
                         y = "Likelihood") +
                    theme_bw() +
                  theme(
                  plot.background = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.grid.major = element_blank(),
                  text=element_text(size = 12, family="serif"))

  
#four plots in a figure
  four_p <- ggpubr::ggarrange(heldout_plot, resid_plot, exclus_plot, coher_plot,
                            ncol = 2, nrow = 2, align = "hv")
  
  four_p
  ##ggsave(plot = four_p, width = 9, height = 5, dpi = 300, filename = "searchK_threshold5_2to40_ggplot.png")  
```


# STM Estimation 

- lower threshold set as 5
- five models are estimated based on the searchK results


```{r model-est}
#model with 6 topics
anger_fit_6t <- stm(documents = out$documents, vocab = out$vocab, K = 6, prevalence = ~Condition, max.em.its = 100, data = out$meta, init.type = "Spectral", seed = 110)

#model with 12 topics
anger_fit_12t <- stm(documents = out$documents, vocab = out$vocab, K = 12, prevalence = ~Condition, max.em.its = 100, data = out$meta, init.type = "Spectral", seed = 110)

#model with 16 topics
anger_fit_16t <- stm(documents = out$documents, vocab = out$vocab, K = 16, prevalence = ~Condition, max.em.its = 100, data = out$meta, init.type = "Spectral", seed = 110)

#model with 22 topics
anger_fit_22t <- stm(documents = out$documents, vocab = out$vocab, K = 22, prevalence = ~Condition, max.em.its = 100, data = out$meta, init.type = "Spectral", seed = 110)

#model with 25 topics
anger_fit_25t <- stm(documents = out$documents, vocab = out$vocab, K = 25, prevalence = ~Condition, max.em.its = 75, data = out$meta, init.type = "Spectral", seed = 110)

```

# Evaluating Topics

## Topic Quality

```{r}
#using the internal plotting function
  topicQuality(model = anger_fit_6t, documents = out$documents) 
  topicQuality(model = anger_fit_12t, documents = out$documents)
  topicQuality(model = anger_fit_16t, documents = out$documents)
  topicQuality(model = anger_fit_22t, documents = out$documents)
  topicQuality(model = anger_fit_25t, documents = out$documents)

#using my own function

t1_6 <- topic_quality_plot(anger_fit_6t)
t1_12 <- topic_quality_plot(anger_fit_12t)
t1_16 <- topic_quality_plot(anger_fit_16t)
t1_22 <- topic_quality_plot(anger_fit_22t)
t1_25 <- topic_quality_plot(anger_fit_25t)


#five plots in a figure
  five_p <- ggpubr::ggarrange(t1_6, t1_12, t1_16, t1_22, t1_25,
                            ncol = 3, nrow = 2, align = "hv")
  
  five_p
  ##ggsave(plot = five_p, width = 8, height = 7, dpi = 300, filename = "topic_quality_ggplot.png") 

```


## Topic Proportions

```{r}

six_t_plot <- stm_plot(anger_fit_6t, perc_bar_lim = 0.45, ylabel = "no")
twelve_t_plot <- stm_plot(anger_fit_12t, perc_bar_lim = 0.45, ylabel = "no")
sixteen_t_plot <- stm_plot(anger_fit_16t, perc_bar_lim = 0.45, ylabel = "no")
twentytwo_t_plot <- stm_plot(anger_fit_22t, perc_bar_lim = 0.45, ylabel = "yes")
twentyfive_t_plot <- stm_plot(anger_fit_25t, perc_bar_lim = 0.45, ylabel = "yes")

six_t_plot
twelve_t_plot
sixteen_t_plot
twentytwo_t_plot
twentyfive_t_plot

#three plots in a figure
  three_p <- ggpubr::ggarrange(six_t_plot, sixteen_t_plot, twentyfive_t_plot,
                            ncol = 1, nrow = 3, align = "hv",
                            heights = c(0.185, 0.34, 0.525))
  #three_p
  
  ##ggsave(plot = three_p, width = 8.5, height = 10.5, dpi = 300, filename = "three_topwords_w_proportions.png") 
  
#another three plots in a figure
  three_p2 <- ggpubr::ggarrange(six_t_plot, twelve_t_plot, twentytwo_t_plot,
                            ncol = 1, nrow = 3, align = "hv",
                            heights = c(0.20, 0.337, 0.523))
  #three_p2
  
  ##ggsave(plot = three_p2, width = 8.5, height = 11, dpi = 300, filename = "another_three_topwords_w_proportions.png")   
```


# Estimate Covariate Relations (12T Model)

```{r}
out$meta$Condition <- as.factor(out$meta$Condition)

#6-topic model
prep_6t <- estimateEffect(formula = 1:6 ~ Condition, anger_fit_6t, meta = out$meta, uncertainty = "Global")
#summary(prep_6t, topics = c(1:6))

#12-topic model
prep_12t <- estimateEffect(1:12 ~ Condition, anger_fit_12t, meta = out$meta, uncertainty = "Global")
#summary(prep_12t, topics = c(1:12))

#16-topic model
prep_16t <- estimateEffect(1:16 ~ Condition, anger_fit_16t, meta = out$meta, uncertainty = "Global")
#summary(prep_16t, topics = c(1:16))

#22-topic model
prep_22t <- estimateEffect(1:22 ~ Condition, anger_fit_22t, meta = out$meta, uncertainty = "Global")
#summary(prep_22t, topics = c(1:22))

#25-topic model
prep_25t <- estimateEffect(1:25 ~ Condition, anger_fit_25t, meta = out$meta, uncertainty = "Global")
#summary(prep_25t, topics = c(1:25))

report_stm_cov_effect(prep_6t)
report_stm_cov_effect(prep_12t)
#report_stm_cov_effect(prep_12t, save_as_csv = "stm_12t_results.csv") #saves results as csv (for manuscript)
report_stm_cov_effect(prep_16t)
report_stm_cov_effect(prep_22t)
report_stm_cov_effect(prep_25t)

```


## compare thetas (for 12T model) using t-test

```{r}
df_w_stm_12t <- cbind(processed$meta, anger_fit_12t[["theta"]]) 
last_col_number <- ncol(df_w_stm_12t)
t_number <- last_col_number-3
colnames(df_w_stm_12t)[4:last_col_number] <- c(paste0("Topic ", 1:t_number))

glimpse(df_w_stm_12t)

#####Descriptives
#Collapsed across conditions
df_w_stm_12t %>% 
  select(contains("topic")) %>% 
  psych::describe() %>% 
  print(digits = 3) #the outcome matches the Topic Proportions

#For each condition 
df_w_stm_12t %>% 
  select(Condition, contains("topic")) %>% 
  psych::describeBy(group = df_w_stm_12t$Condition) %>% 
  print(digits = 3)

#For each condition (only the means)
df_w_stm_12t %>% 
  group_by(Condition) %>% 
  summarise(across(contains("topic"), mean, na.rm=TRUE))

#####t-tests
for (i in 1:12){
  dv_col_number <- i + 3
  
    iv <- df_w_stm_12t$Condition
    dv <- df_w_stm_12t[,dv_col_number]
  
      ttest <- t.test(dv ~ iv)
      effect_size <- effectsize::cohens_d(dv ~ iv, pooled_sd = FALSE)
      t <- ttest[["statistic"]] %>% round(2)
      df <- ttest[["parameter"]] %>% round(1)
      original_p <- ttest[["p.value"]] %>% round(3)
      p <- if_else(original_p >= .001, paste0("= ", as.character(original_p)), "< .001")
      d <- effect_size[1,1] %>% round(2)    
      topic <- (colnames(df_w_stm_12t))[dv_col_number]
      print(paste0(topic, ": t(", df, ") = ", t, ", p ", p, ", d = ", d))

}
```

## Histogram of thetas

```{r}
df_w_stm_12t_long <- df_w_stm_12t %>% 
                        pivot_longer(cols = contains("topic"),
                                     names_to = "Topic",
                                     values_to = "Theta")

    #all together
    all_histograms <- df_w_stm_12t_long %>% 
                                mutate(across(Topic, factor, levels = c(paste("Topic", 1:12)))) %>% 
                                group_by(Topic) %>% 
                                mutate(med_theta = median(Theta)) %>% 
      ggplot(., aes(x = Theta, fill = Theta)) + 
                      geom_histogram(aes(y = ..density..), binwidth = 0.004, color = "darkblue", fill = "grey", size = 0.3) +
                      #geom_density(fill = "lightblue", color = "grey", alpha = 0.3) +
                      geom_vline(aes(xintercept = med_theta), color = "red", linetype = "dashed", size = 0.3) +
                      facet_wrap ( ~ Topic) +
                      xlim(0, 0.2) +
                      labs(x = expression(paste("Per-Document-Per-Topic Probability (", theta, ")")),
                           y = "Density") +
                      theme_minimal() +
                      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                            text=element_text(size = 12, family="serif"),
                            axis.title.x = element_text(margin=margin(t=10)))
    
    all_histograms
    ##ggsave(plot = all_histograms, width = 6.5, height = 5, dpi = 600, filename = "all_theta_histograms.png") 
    
    
    #different for the two conditions
    ggplot(df_w_stm_12t_long, aes(x = Theta, fill = Condition)) + 
                    geom_histogram(aes(y = ..density..), binwidth = 0.005, alpha = 0.3) +
                    geom_density(alpha = 0.4) +
                    facet_wrap ( ~ Topic) +
                    xlim(0, 0.4) +
                    theme_minimal() +
                    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

```


# Words and Documents highly Associated with Topics (12T)

```{r}
#look at collections of words that are associated with topics
  labels_6t <- labelTopics(anger_fit_6t, c(1:6), n = 10)
  labels_12t <- labelTopics(anger_fit_12t, c(1:12), n = 10)
  labels_16t <- labelTopics(anger_fit_16t, c(1:16), n = 10)
  labels_22t <- labelTopics(anger_fit_22t, c(1:22), n = 10)
  labels_25t <- labelTopics(anger_fit_25t, c(1:25), n = 10)

  labels_6t
  labels_12t
  labels_16t
  labels_22t
  labels_25t
  
  #Example: Plotting and saving highest probable words using the internal function
      #png(file="topic_words_6t.png", width=1100, height=500)
      #par(mfrow = c(1, 2))
      #plot(anger_fit_6t, type = "labels", topics = 1:3, labeltype = "frex", n = 10)
      #plot(anger_fit_6t, type = "labels", topics = 4:6, labeltype = "frex", n = 10)
      #dev.off()
      
  
  
  #Possibly more interpretable: the LIFT method
      #cbind(paste("topic", 1:6), labels_6t[["lift"]]) %>% knitr::kable() %>% kableExtra::kable_styling()

#Create a table to report the top words associated with 12 topics  
    topic_word_list <- data.frame(Topic = NA, Method = NA,
                                  w1 = NA, w2 = NA, w3 = NA, w4 = NA, w5 = NA, 
                                  w6 = NA, w7 = NA, w8 = NA, w9 = NA, w10 = NA)  
        
    for (i in 1:12){  
      topic <- paste("Topic", i)
      prob_words <- labels_12t[["prob"]][i,]
      frex_words <- labels_12t[["frex"]][i,]
      lift_words <- labels_12t[["lift"]][i,]
      
      prob_row <- c(topic, "Prob.", prob_words)
      frex_row <- c(topic, "FREX", frex_words)
      lift_row <- c(topic, "LIFT", lift_words)
      
      topic_word_list <- rbind(topic_word_list,
                               prob_row,
                               frex_row,
                               lift_row)
      
    }  
    
    
    topic_word_list_clean <- topic_word_list[-1,] %>% 
                          tidyr::unite("stem_words", w1:w10, sep = ", ", remove = TRUE)
    
    #save for the manuscript table
    #write.csv(topic_word_list_clean, "top_word_for_12_topics.csv")
    
    
#Examining documents that are highly associated with topics
  
  example_topic1 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 1)$docs[[1]]
  example_topic2 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 2)$docs[[1]]
  example_topic3 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 3)$docs[[1]]
  example_topic4 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 4)$docs[[1]]
  example_topic5 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 5)$docs[[1]] 
  example_topic6 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 6)$docs[[1]]
  example_topic7 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 7)$docs[[1]]
  example_topic8 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 8)$docs[[1]]
  example_topic9 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 9)$docs[[1]]
  example_topic10 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 10)$docs[[1]]
  example_topic11 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 11)$docs[[1]]
  example_topic12 <- findThoughts(anger_fit_12t, texts = meta$all_narratives, n = 3, topics = 12)$docs[[1]]
  
  example_topic1
  example_topic2
  example_topic3
  example_topic4
  example_topic5
  example_topic6
  example_topic7
  example_topic8
  example_topic9
  example_topic10
  example_topic11
  example_topic12

```

# Figures of High Probability Words

## Bar charts

```{r}
top_terms <- tidy(anger_fit_12t, matrix = "beta") %>% 
                group_by(topic) %>% 
                slice_max(beta, n = 10) %>% 
                ungroup() %>% 
                arrange(topic, -beta)

top_terms %>% 
  mutate(term = forcats::fct_reorder(term, beta)) %>% 
      ggplot(aes(beta, term, fill = factor(topic))) + 
      geom_col(show.legend = FALSE) + 
      facet_wrap(~ topic, scales = "free")

```


## Word cloud

```{r}
top_25_terms <- tidy(anger_fit_12t, matrix = "beta") %>% 
                group_by(topic) %>% 
                slice_max(beta, n = 25) %>% 
                ungroup() %>% 
                arrange(topic, -beta)


top_25_terms %>% 
  mutate(term = forcats::fct_reorder(term, beta)) %>% 
    ggplot(aes(label = term, size = beta, color = factor(sample.int(10, nrow(top_25_terms), replace = TRUE)))) +
    ggwordcloud::geom_text_wordcloud(rm_outside = TRUE) +
    scale_size_area(max_size = 12) +
    facet_wrap(~ topic) +
    theme_minimal()
```



# Topic correlations

```{r}
#6t
  topic_correlation_6t <- topicCorr(anger_fit_6t)
  plot(topic_correlation_6t)

#12t
  topic_correlation_12t <- topicCorr(anger_fit_12t)
  plot(topic_correlation_12t)

#16t
  topic_correlation_16t <- topicCorr(anger_fit_16t)
  plot(topic_correlation_16t)

#22t
  topic_correlation_22t <- topicCorr(anger_fit_22t)
  plot(topic_correlation_22t)
    
#25t
  topic_correlation_25t <- topicCorr(anger_fit_25t)
  plot(topic_correlation_25t)

```