Food and animals

Suppose you have the following set of sentences:

1. I ate a banana and spinach smoothie for breakfast.
2. I like to eat broccoli and bananas.
3. Chinchillas and kittens are cute.
4. My sister adopted a kitten yesterday.
5. Look at this cute hamster munching on a piece of broccoli.

Latent Dirichlet allocation is a way of automatically discovering topics that these sentences contain. For example, given these sentences and asked for 2 topics, LDA might produce something like

• Sentences 1 and 2: 100% Topic A
• Sentences 3 and 4: 100% Topic B

• Sentence 5: 60% Topic A, 40% Topic B

• Topic A: 30% broccoli, 15% bananas, 10% breakfast, 10% munching, …

• Topic B: 20% chinchillas, 20% kittens, 20% cute, 15% hamster, …

You could infer that topic A is a topic about food, and topic B is a topic about cute animals. But LDA does not explicitly identify topics in this manner. All it can do is tell you the probability that specific words are associated with the topic.

An LDA document structure

LDA represents documents as mixtures of topics that spit out words with certain probabilities. It assumes that documents are produced in the following fashion: when writing each document, you

• Decide on the number of words \(N\) the document will have
• Choose a topic mixture for the document (according to a Dirichlet probability distribution over a fixed set of \(K\) topics). For example, assuming that we have the two food and cute animal topics above, you might choose the document to consist of 1/3 food and 2/3 cute animals.
• Generate each word in the document by:
  • First picking a topic (according to the distribution that you sampled above; for example, you might pick the food topic with 1/3 probability and the cute animals topic with 2/3 probability).
  • Then using the topic to generate the word itself (according to the topic’s multinomial distribution). For instance, the food topic might output the word “broccoli” with 30% probability, “bananas” with 15% probability, and so on.

Assuming this generative model for a collection of documents, LDA then tries to backtrack from the documents to find a set of topics that are likely to have generated the collection.

Latent Dirichlet allocation with the topicmodels package

Right now this data frame is in a tidy form, with one-term-per-document-per-row. However, the topicmodels package requires a DocumentTermMatrix (from the tm package). We can cast a one-token-per-row table into a DocumentTermMatrix with cast_dtm():

chapters_dtm <- word_counts %>%
  cast_dtm(title_chapter, word, n)

chapters_dtm

## <<DocumentTermMatrix (documents: 193, terms: 18215)>>
## Non-/sparse entries: 104722/3410773
## Sparsity           : 97%
## Maximal term length: 19
## Weighting          : term frequency (tf)

Now we are ready to use the topicmodels package to create a four topic LDA model.

library(topicmodels)
chapters_lda <- LDA(chapters_dtm, k = 4, control = list(seed = 1234))
chapters_lda

## A LDA_VEM topic model with 4 topics.

• In this case we know there are four topics because there are four books; this is the value of knowing the latent topic structure.
• seed = 1234 sets the starting point for the random iteration process. If we don’t set a consistent seed, each time we run the script we may estimate slightly different models.

Now tidytext gives us the option of returning to a tidy analysis, using the tidy() and augment() verbs borrowed from the broom package. In particular, we start with the tidy() verb.

library(tidytext)

chapters_lda_td <- tidy(chapters_lda)
chapters_lda_td

## # A tibble: 72,860 x 3
##    topic    term         beta
##    <int>   <chr>        <dbl>
##  1     1     joe 1.436612e-17
##  2     2     joe 5.962111e-61
##  3     3     joe 9.881855e-25
##  4     4     joe 1.447329e-02
##  5     1   biddy 5.139275e-28
##  6     2   biddy 5.022015e-73
##  7     3   biddy 4.307280e-48
##  8     4   biddy 4.775557e-03
##  9     1 estella 2.431464e-06
## 10     2 estella 4.323253e-68
## # ... with 72,850 more rows

Notice that this has turned the model into a one-topic-per-term-per-row format. For each combination the model has beta (\(\beta\)), the probability of that term being generated from that topic.

We could use top_n() from dplyr to find the top 5 terms within each topic:

top_terms <- chapters_lda_td %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
top_terms

## # A tibble: 20 x 3
##    topic      term        beta
##    <int>     <chr>       <dbl>
##  1     1 elizabeth 0.014101270
##  2     1     darcy 0.008810341
##  3     1      miss 0.008708777
##  4     1    bennet 0.006944344
##  5     1      jane 0.006494613
##  6     2   captain 0.015510635
##  7     2  nautilus 0.013051927
##  8     2       sea 0.008843483
##  9     2      nemo 0.008709651
## 10     2       ned 0.008031955
## 11     3    people 0.006785987
## 12     3  martians 0.006456394
## 13     3      time 0.005343667
## 14     3     black 0.005277449
## 15     3     night 0.004491174
## 16     4       joe 0.014473289
## 17     4      time 0.006852889
## 18     4       pip 0.006828209
## 19     4    looked 0.006366418
## 20     4      miss 0.006232761

This model lends itself to a visualization:

top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_bar(alpha = 0.8, stat = "identity", show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free") +
  coord_flip()

• These topics are pretty clearly associated with the four books:
  • “nemo”, “sea”, and “nautilus” belongs to Twenty Thousand Leagues Under the Sea
  • “jane”, “darcy”, and “elizabeth” belongs to Pride and Prejudice
  • “pip” and “joe” from Great Expectations
  • “martians”, “black”, and “night” from The War of the Worlds
• Also note that LDA() does not assign any label to each topic. They are simply topics 1, 2, 3, and 4. We can infer these are associated with each book, but it is merely our inference.

Per-document classification

Each chapter was a “document” in this analysis. Thus, we may want to know which topics are associated with each document. Can we put the chapters back together in the correct books?

chapters_lda_gamma <- tidy(chapters_lda, matrix = "gamma")
chapters_lda_gamma

## # A tibble: 772 x 3
##                    document topic        gamma
##                       <chr> <int>        <dbl>
##  1    Great Expectations_57     1 1.338547e-05
##  2     Great Expectations_7     1 1.456215e-05
##  3    Great Expectations_17     1 2.096237e-05
##  4    Great Expectations_27     1 1.900804e-05
##  5    Great Expectations_38     1 3.552749e-01
##  6     Great Expectations_2     1 1.706715e-05
##  7    Great Expectations_23     1 5.470853e-01
##  8    Great Expectations_15     1 1.243917e-02
##  9    Great Expectations_18     1 1.259492e-05
## 10 The War of the Worlds_16     1 1.073638e-05
## # ... with 762 more rows

Setting matrix = "gamma" returns a tidied version with one-document-per-topic-per-row. Now that we have these document classifiations, we can see how well our unsupervised learning did at distinguishing the four books. First we re-separate the document name into title and chapter:

chapters_lda_gamma <- chapters_lda_gamma %>%
  separate(document, c("title", "chapter"), sep = "_", convert = TRUE)
chapters_lda_gamma

## # A tibble: 772 x 4
##                    title chapter topic        gamma
##  *                 <chr>   <int> <int>        <dbl>
##  1    Great Expectations      57     1 1.338547e-05
##  2    Great Expectations       7     1 1.456215e-05
##  3    Great Expectations      17     1 2.096237e-05
##  4    Great Expectations      27     1 1.900804e-05
##  5    Great Expectations      38     1 3.552749e-01
##  6    Great Expectations       2     1 1.706715e-05
##  7    Great Expectations      23     1 5.470853e-01
##  8    Great Expectations      15     1 1.243917e-02
##  9    Great Expectations      18     1 1.259492e-05
## 10 The War of the Worlds      16     1 1.073638e-05
## # ... with 762 more rows

Then we examine what fraction of chapters we got right for each:

# reorder titles in order of topic 1, topic 2, etc before plotting
chapters_lda_gamma %>%
  mutate(title = reorder(title, gamma * topic)) %>%
  ggplot(aes(factor(topic), gamma)) +
  geom_boxplot() +
  facet_wrap(~ title)

We notice that almost all of the chapters from Pride and Prejudice, War of the Worlds, and Twenty Thousand Leagues Under the Sea were uniquely identified as a single topic each.

chapter_classifications <- chapters_lda_gamma %>%
  group_by(title, chapter) %>%
  top_n(1, gamma) %>%
  ungroup() %>%
  arrange(gamma)

chapter_classifications

## # A tibble: 193 x 4
##                 title chapter topic     gamma
##                 <chr>   <int> <int>     <dbl>
##  1 Great Expectations      54     3 0.4812041
##  2 Great Expectations      22     4 0.5362582
##  3 Great Expectations      23     1 0.5470853
##  4 Great Expectations      31     4 0.5473876
##  5 Great Expectations      33     4 0.5689223
##  6 Great Expectations      47     4 0.5799048
##  7 Great Expectations      56     4 0.6060511
##  8 Great Expectations      38     4 0.6446998
##  9 Great Expectations       3     4 0.6602340
## 10 Great Expectations      11     4 0.6675174
## # ... with 183 more rows

We can determine this by finding the consensus book for each, which we note is correct based on our earlier visualization:

book_topics <- chapter_classifications %>%
  count(title, topic) %>%
  group_by(title) %>%
  top_n(1, n) %>%
  ungroup() %>%
  transmute(consensus = title, topic)

book_topics

## # A tibble: 4 x 2
##                               consensus topic
##                                   <chr> <int>
## 1                    Great Expectations     4
## 2                   Pride and Prejudice     1
## 3                 The War of the Worlds     3
## 4 Twenty Thousand Leagues under the Sea     2

Then we see which chapters were misidentified:

chapter_classifications %>%
  inner_join(book_topics, by = "topic") %>%
  count(title, consensus) %>%
  knitr::kable()

We see that only a few chapters from Great Expectations were misclassified.

By word assignments: augment

One important step in the topic modeling expectation-maximization algorithm is assigning each word in each document to a topic. The more words in a document are assigned to that topic, generally, the more weight (gamma) will go on that document-topic classification.

We may want to take the original document-word pairs and find which words in each document were assigned to which topic. This is the job of the augment() verb.

assignments <- augment(chapters_lda, data = chapters_dtm)

We can combine this with the consensus book titles to find which words were incorrectly classified.

assignments <- assignments %>%
  separate(document, c("title", "chapter"), sep = "_", convert = TRUE) %>%
  inner_join(book_topics, by = c(".topic" = "topic"))

assignments

## # A tibble: 29,412 x 6
##                  title chapter    term count .topic           consensus
##                  <chr>   <int>   <chr> <dbl>  <dbl>               <chr>
##  1 Pride and Prejudice      47  pocket     1      1 Pride and Prejudice
##  2 Pride and Prejudice      50  pocket     1      1 Pride and Prejudice
##  3 Pride and Prejudice      49  pocket     1      1 Pride and Prejudice
##  4  Great Expectations      38 brother     2      1 Pride and Prejudice
##  5 Pride and Prejudice      43 brother     4      1 Pride and Prejudice
##  6 Pride and Prejudice      18 brother     1      1 Pride and Prejudice
##  7  Great Expectations      22 brother     4      1 Pride and Prejudice
##  8 Pride and Prejudice      45 brother     2      1 Pride and Prejudice
##  9 Pride and Prejudice      16 brother     1      1 Pride and Prejudice
## 10 Pride and Prejudice      10 brother     2      1 Pride and Prejudice
## # ... with 29,402 more rows

We can, for example, create a “confusion matrix” using dplyr::count() and tidyr::spread:

assignments %>%
  count(title, consensus, wt = count) %>%
  spread(consensus, n, fill = 0) %>%
  knitr::kable()

We notice that almost all the words for Pride and Prejudice, Twenty Thousand Leagues Under the Sea, and War of the Worlds were correctly assigned, while Great Expectations had a fair amount of misassignment.

What were the most commonly mistaken words?

wrong_words <- assignments %>%
  filter(title != consensus)

wrong_words

## # A tibble: 3,140 x 6
##                                    title chapter     term count .topic
##                                    <chr>   <int>    <chr> <dbl>  <dbl>
##  1                    Great Expectations      38  brother     2      1
##  2                    Great Expectations      22  brother     4      1
##  3                    Great Expectations      23     miss     2      1
##  4                    Great Expectations      22     miss    23      1
##  5 Twenty Thousand Leagues under the Sea       8     miss     1      1
##  6                    Great Expectations      31     miss     1      1
##  7                    Great Expectations       5 sergeant    37      1
##  8                    Great Expectations      22   people     3      1
##  9                    Great Expectations      31   people     2      1
## 10                    Great Expectations      23     dear    10      1
## # ... with 3,130 more rows, and 1 more variables: consensus <chr>

wrong_words %>%
  count(title, consensus, term, wt = count) %>%
  ungroup() %>%
  arrange(desc(n))

## # A tibble: 2,291 x 4
##                 title           consensus      term     n
##                 <chr>               <chr>     <chr> <dbl>
##  1 Great Expectations Pride and Prejudice      love    44
##  2 Great Expectations Pride and Prejudice  sergeant    37
##  3 Great Expectations Pride and Prejudice      lady    32
##  4 Great Expectations Pride and Prejudice      miss    26
##  5 Great Expectations Pride and Prejudice    father    19
##  6 Great Expectations Pride and Prejudice      baby    18
##  7 Great Expectations Pride and Prejudice   flopson    18
##  8 Great Expectations Pride and Prejudice    family    16
##  9 Great Expectations Pride and Prejudice   replied    16
## 10 Great Expectations Pride and Prejudice attention    14
## # ... with 2,281 more rows

Notice the word “flopson” here; these wrong words do not necessarily appear in the novels they were misassigned to. Indeed, we can confirm “flopson” appears only in Great Expectations:

word_counts %>%
  filter(word == "flopson")

## # A tibble: 3 x 3
##           title_chapter    word     n
##                   <chr>   <chr> <int>
## 1 Great Expectations_22 flopson    10
## 2 Great Expectations_23 flopson     7
## 3 Great Expectations_33 flopson     1

The algorithm is stochastic and iterative, and it can accidentally land on a topic that spans multiple books.

Associated Press articles

The topicmodels package includes a document-term matrix of a sample of articles published by the Associated Press in 1992. Let’s load them into R and convert them to a tidy format.

data("AssociatedPress", package = "topicmodels")

ap_td <- tidy(AssociatedPress)
ap_td

## # A tibble: 302,031 x 3
##    document       term count
##       <int>      <chr> <dbl>
##  1        1     adding     1
##  2        1      adult     2
##  3        1        ago     1
##  4        1    alcohol     1
##  5        1  allegedly     1
##  6        1      allen     1
##  7        1 apparently     2
##  8        1   appeared     1
##  9        1   arrested     1
## 10        1    assault     1
## # ... with 302,021 more rows

AssociatedPress is originally in a document-term matrix, exactly what we need for topic modeling. Why tidy it first? Because the original document-term matrix contains stop words - we want to remove them before modeling the data. Let’s remove the stop words, then cast the data back into a document-term matrix.

ap_dtm <- ap_td %>%
  anti_join(stop_words, by = c(term = "word")) %>%
  cast_dtm(document, term, count)
ap_dtm

## <<DocumentTermMatrix (documents: 2246, terms: 10134)>>
## Non-/sparse entries: 259208/22501756
## Sparsity           : 99%
## Maximal term length: 18
## Weighting          : term frequency (tf)

Selecting \(k\)

Remember that for LDA, you need to specify in advance the number of topics in the underlying topic structure.

ap_lda <- LDA(ap_dtm, k = 4, control = list(seed = 1234))
ap_lda

## A LDA_VEM topic model with 4 topics.

ap_lda_td <- tidy(ap_lda)

top_terms <- ap_lda_td %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
top_terms

## # A tibble: 20 x 3
##    topic       term        beta
##    <int>      <chr>       <dbl>
##  1     1     police 0.010192765
##  2     1     people 0.008563645
##  3     1  officials 0.004779668
##  4     1       city 0.003985224
##  5     1     killed 0.003730190
##  6     2     soviet 0.011316539
##  7     2 government 0.009824318
##  8     2  president 0.008803274
##  9     2     united 0.007751739
## 10     2      party 0.005982017
## 11     3    percent 0.019093039
## 12     3    million 0.012657406
## 13     3    billion 0.008941798
## 14     3     market 0.006493595
## 15     3    company 0.006010625
## 16     4      court 0.006149978
## 17     4       bush 0.004824762
## 18     4     people 0.004306685
## 19     4    dukakis 0.004189072
## 20     4  president 0.004169274

top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_bar(alpha = 0.8, stat = "identity", show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", ncol = 2) +
  coord_flip()

ap_lda <- LDA(ap_dtm, k = 12, control = list(seed = 1234))
ap_lda

## A LDA_VEM topic model with 12 topics.

ap_lda_td <- tidy(ap_lda)

top_terms <- ap_lda_td %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
top_terms

## # A tibble: 60 x 3
##    topic       term        beta
##    <int>      <chr>       <dbl>
##  1     1     people 0.004966579
##  2     1       dont 0.003962532
##  3     1        air 0.003910778
##  4     1       time 0.003768088
##  5     1       york 0.003672212
##  6     2     soviet 0.017475311
##  7     2        aid 0.010831426
##  8     2    million 0.006403843
##  9     2 government 0.006218165
## 10     2       corn 0.005717200
## # ... with 50 more rows

top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_bar(alpha = 0.8, stat = "identity", show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", ncol = 3) +
  coord_flip()

Perplexity

Well, sort of. Some aspects of LDA are driven by gut-thinking (or perhaps truthiness). However we can have some help. Perplexity is a statistical measure of how well a probability model predicts a sample. As applied to LDA, for a given value of \(k\), you estimate the LDA model. Then given the theoretical word distributions represented by the topics, compare that to the actual topic mixtures, or distribution of words in your documents.

topicmodels includes the function perplexity() which calculates this value for a given model.

perplexity(ap_lda)

## [1] 2277.876

However, the statistic is somewhat meaningless on its own. The benefit of this statistic comes in comparing perplexity across different models with varying \(k\)s. The model with the lowest perplexity is generally considered the “best”.

Let’s estimate a series of LDA models on the Associated Press dataset. Here I make use of purrr and the map() functions to iteratively generate a series of LDA models for the AP corpus, using a different number of topics in each model.¹

n_topics <- c(2, 4, 10, 20, 50, 100)
ap_lda_compare <- n_topics %>%
  map(LDA, x = ap_dtm, control = list(seed = 1234))

n_topics <- c(2, 4, 10, 20, 50, 100)

# takes forever to estimate this model - store results and use if available
if(file.exists("extras/ap_lda_compare.Rdata")){
  load(file = "extras/ap_lda_compare.Rdata")
} else{
  ap_lda_compare <- n_topics %>%
    map(LDA, x = ap_dtm, control = list(seed = 1234))
  save(ap_lda_compare, file = "extras/ap_lda_compare.Rdata")
}

data_frame(k = n_topics,
           perplex = map_dbl(ap_lda_compare, perplexity)) %>%
  ggplot(aes(k, perplex)) +
  geom_point() +
  geom_line() +
  labs(title = "Evaluating LDA topic models",
       subtitle = "Optimal number of topics (smaller is better)",
       x = "Number of topics",
       y = "Perplexity")

It looks like the 100-topic model has the lowest perplexity score. What kind of topics does this generate? Let’s look just at the first 12 topics produced by the model (ggplot2 has difficulty rendering a graph for 100 separate facets):

ap_lda_td <- tidy(ap_lda_compare[[6]])

top_terms <- ap_lda_td %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
top_terms

## # A tibble: 500 x 3
##    topic       term        beta
##    <int>      <chr>       <dbl>
##  1     1  president 0.008017375
##  2     1        oil 0.005617035
##  3     1     people 0.005532095
##  4     1    embassy 0.005264867
##  5     1 television 0.005183352
##  6     2 convention 0.016274726
##  7     2       york 0.010238948
##  8     2    dukakis 0.008491059
##  9     2   national 0.006934871
## 10     2    jackson 0.006472480
## # ... with 490 more rows

top_terms %>%
  filter(topic <= 12) %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_bar(alpha = 0.8, stat = "identity", show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", ncol = 3) +
  coord_flip()

We are getting even more specific topics now. The question becomes how would we present these results and use them in an informative way? Not to mention perplexity was still dropping at \(k=100\) - would \(k=200\) generate an even lower perplexity score?²

Again, this is where your intuition and domain knowledge as a researcher is important. You can use perplexity as one data point in your decision process, but a lot of the time it helps to simply look at the topics themselves and the highest probability words associated with each one to determine if the structure makes sense. If you have a known topic structure you can compare it to (such as the books example above), this can also be useful.