06-topic-models.Rmd

# 主题模型 {#topicmodeling}

```{r echo = FALSE}
library(knitr)
opts_chunk$set(message = FALSE, warning = FALSE, cache = TRUE)
options(width = 100, dplyr.width = 150)
library(tidytext)
library(tidyr)
library(dplyr)
library(jiebaR)
library(ggplot2)
library(methods)
library(scales)
library(showtext)
showtext_auto(enable = TRUE)
pdf()
theme_zh <- theme_light() +
  theme(text = element_text(family = "wqy-microhei"))
theme_set(theme_zh)
``` 

在文本挖掘中，我们经常有多个包含多个文档的集合，如博客或新闻文章，我们想将其自然分组，以便分开阅读。主题模型是对这样的文档进行无监督分类的一种方法，类似于数字型数据的聚类，即使不知道需要寻找什么也可以找到自然的组别。

隐含狄利克雷分布（Latent Dirichlet allocation，LDA）是一种特别流行的拟合主题模型的方法。它将每个文档视为多个主题的混合，而每个主题又是多个词的混合。LDA 允许文档在内容层面互有重合，而不是被分成离散的多个组，这在某种意义上体现了自然语言的典型用法。

```{r tidyflowchartch6, echo = FALSE, out.width = '100%', fig.cap = "应用主题模型进行文本分析的流程图。"}
knitr::include_graphics("images/tidyflow-ch-6.png")
```

如图 \@ref(fig:tidyflowchartch6) 所示，我们可以使用 tidy 文本原则得到主题模型，通过在全书中使用的同一套 tidy 工具。在本章中，我们将学习 [topicmodels 包](https://cran.r-project.org/package=topicmodels) [@R-topicmodels] 中的 `LDA` 对象，特别是将这些模型 tidy 化，以便在 ggplot2 和 dplyr 中加以操作。我们还将探索一个多本书籍章节聚类的例子，从中可以看到主题模型能基于文本内容“学习”到不同的书有何不同。

## 隐含狄利克雷分布

隐含狄利克雷分布是主题模型最通行的算法之一。无需触及其模型背后的数学，我们可以从两个原则出发理解 LDA。

* **每个文档都是主题的混合** 设想每个文档可能含有的词来自特定比例的几个主题。比如，针对一个双主题的模型，我们可以说“文档 1 是 90% 主题 A 和 10% 主题 B，而文档 2 是 30% 主题 A 和 70% 主题 B。”
* **每个主题都是词的混合** 比如，可以设想一个美国新闻的双主题模型，一个主题是“政治”，一个是“娱乐”。政治主题中最常见的词可能是“总统”“议院”“政府”，而组成娱乐主题的词如“电影”“电视”“演员”。重要的是，词可以被主题共用，比如“预算”就可能同等地出现在两个主题中。

LDA 是同时估算这两件事的数学方法：找到与每个主题相关联的词的混合，同时确定描述每个文档的主题的混合。这个算法的实现已经有很多种，我们将深度探索其中之一。

在章 \@ref(dtm) 中我们简要地介绍了 topicmodels 包提供的 `AssociatedPress` 数据集作为 DocumentTermMatrix 的一个例子。这是一个2246篇新闻文章的集合，来自美国的一个通讯社，主要发表于1988年前后。

```{r}
library(topicmodels)

data("AssociatedPress")
AssociatedPress
```

我们可以使用 topicmodels 包中的 `LDA()` 函数，设定 `k = 2` 以创建一个双主题的 LDA 模型。

```{block, type = "rmdnote"}
Almost any topic model in practice will use a larger `k`, but we will soon see that this analysis approach extends to a larger number of topics.
```

这个函数返回一个对象，包含拟合模型的全部细节，词如何与主题相关联，而主题如何与文档相关联。

```{r ap_lda}
# set a seed so that the output of the model is predictable
ap_lda <- LDA(AssociatedPress, k = 2, control = list(seed = 1234))
ap_lda
```

拟合模型比较简单，余下的分析将包括使用 tidytext 包探索和解释模型。

### 词-主题概率

在章 \@ref(dtm) 中我们介绍了 `tidy()` 方法，最初来自 broom 包 [@R-broom]，可以将模型对象 tidy 化。tidytext 包提供这个方法从模型中提取每主题每词的概率，称为 $\beta$（beta）。

```{r ap_topics}
library(tidytext)

ap_topics <- tidy(ap_lda, matrix = "beta")
ap_topics
```

注意这把模型变成了一个主题每术语每行的格式。对于每个组合，模型计算了该术语来自该主题的概率。比如，术语 `r ap_topics$term[1]` 有 $`r ap_topics$beta[1]`$ 概率来自主题 `r ap_topics$topic[1]`，而有 $`r ap_topics$beta[2]`$ 概率来自主题 `r ap_topics$topic[2]`。

我们可以使用 dplyr 的 `top_n()` 列出每个主题最常见的10个术语。作为一个 tidy 数据框，这很容易由 ggplot2 可视化（图 \@ref(fig:aptoptermsplot)）。

```{r aptoptermsplot, dependson = "ap_topics", fig.height=5, fig.width=7, fig.cap = "每个主题中这些术语最常见"}
ap_top_terms <- ap_topics %>%
  group_by(topic) %>%
  top_n(10, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
ap_top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free") +
  coord_flip()
```

可视化可以帮助我们理解从文章中提取出的两个主题。主题1中最常见的词包括“percent”“million”“billion”“company”，提示可能代表了商业或金融新闻。主题2最常见的词有“president”“government”“soviet”，提示这个主题代表了政治新闻。一个重要的观察结果是两个主题中有相同的词，如“new”和“people”，在两个主题中都常见。这是主题模型相比“硬性聚类”方法的一个优点：使用自然语言的主题在用词上可能会有交叉。

另外，我们可以考虑 主题1和主题2间 $\beta$ 有 *最大距离* 的术语。这可以通过对数比例估算：$\log_2(\frac{\beta_2}{\beta_1})$ （对数比例可以使距离均匀化：两倍的 $\beta_2$ 即对数比例1，而两倍 $\beta_1$ 则是-1）。要限制为比较相关的词的集合，我们可以过滤相对常见的词，如在至少一个主题中 $\beta$ 超过1/1000。

```{r beta_spread}
beta_spread <- ap_topics %>%
  mutate(topic = paste0("topic", topic)) %>%
  spread(topic, beta) %>%
  filter(topic1 > .001 | topic2 > .001) %>%
  mutate(log_ratio = log2(topic2 / topic1))
beta_spread
```

在两个主题中距离最大的词见图 \@ref(fig:topiccompare)。

(ref:topiccap) 主题2与主题1 $\beta$ 差值最大的词

```{r topiccompare, dependson = "beta_spread", fig.cap = "(ref:topiccap)", echo = FALSE}
beta_spread %>%
  group_by(direction = log_ratio > 0) %>%
  top_n(10, abs(log_ratio)) %>%
  ungroup() %>%
  mutate(term = reorder(term, log_ratio)) %>%
  ggplot(aes(term, log_ratio)) +
  geom_col() +
  labs(y = "主题2和主题1 beta 值以2为底的对数比例") +
  coord_flip()
```

我们可以看到，主题2的常见词包括政党如“democratic”和“republican”，以及政治家的名字如“dukakis”和“gorbachev”。主题1的特征词更多的是货币如“yen”和“dollar”，还有金融术语如“index”“prices”和“rates”。这帮助我们进一步确认算法识别出的两个主题是政治和金融新闻。

### 文档-主题概率

除了按词的混合估计每个主题，LDA 也建立了作为主题混合的文档模型。我们检查一下每文档每主题的概率，称作 $\gamma$（gamma），作为 `matrix = "gamma"` 参数传给 `tidy()`。

```{r ap_documents}
ap_documents <- tidy(ap_lda, matrix = "gamma")
ap_documents
```

这里的每一个值都是一个估算的来自该文档的词有多大比例来自该主题。比如，模型估算文档 `r ap_documents$document[1]` 中大约 `r ap_documents$gamma[1]` 的词来自主题 `r ap_documents$topic[1]`。

我们可以看到这些文档中很多都是两个主题的混合，但文档6几乎全部来自主题2，来自主题1的 $\gamma$ 接近0。要检验这个答案，我们可以 `tidy()` 其文档-术语矩阵（见章 \@ref(tidy-dtm)）并检查该文档中最常见的词。

```{r ap_document_6}
tidy(AssociatedPress) %>%
  filter(document == 6) %>%
  arrange(desc(count))
```

基于最常见的词，看起来这篇文章是关于美国政府与巴拿马当时的统治者 Noriega 的，这意味着算法正确地把它分在了主题2中（政治新闻）。

## 例：图书馆大捣乱 {#library-heist}

要检查一个统计方法，很有用的办法是在一个很简单的例子上尝试使用，而且你知道“正确答案”。比如，我们可以收集一组文档，恰好属于6个相互分离的主题，然后运行主题模型，看看算法能否正确区分出6个组。这可以让我们确认方法是有效的，并且感知到其如何运作，何时又会失效。我们用来自古典文学的一些数据进行尝试。

假设有个强盗闯入你的书斋，把几大名著及其（非原作者的）续作撕破了：

* 《水滸傳》《水滸後傳》
* 《西遊記》《後西游記》
* 《紅樓夢》《補紅樓夢》

这个强盗把书沿章回的边缘撕开堆成了一堆。我们怎么才能把这些杂乱的章节按原书整理好呢？这个问题具有挑战性，因为每个单章都 **无标注**：我们不知道什么词可以用来分组。因此我们使用主题模型来发现章如何聚类为可分辨的主题，每个主题（想必）就代表一本书。

我们用章 \@ref(tfidf) 中介绍的 gutenbergr 包获取文本。

```{r titles}
titles <- c("水滸傳", "西遊記", "紅樓夢",
            "水滸後傳", "後西游記", "補紅樓夢")
```

```{r eval = FALSE}
library(gutenbergr)

books <- gutenberg_works(title %in% titles, languages = 'zh') %>%
  gutenberg_download(meta_fields = "title")
```

```{r topic_books, echo = FALSE}
load("data/books.rda")
```

作为预处理，我们把每部小说都按回分开，再用 tidytext 的 `unnest_tokens()` 分成每词一行。不幸的是，汉语分词的结果对最终模型影响很大，在不同系统或版本上可能会得到不甚一致的结论，有时差异还非常大。我们把每回当作独立的“文档”，命名为类似 `紅樓夢_1` 或 `水滸傳_11` 这样。在其它应用中，每个文档可以是一篇报纸上的文章，或者博客上的一篇博文。

```{r word_counts, dependson = "topic_books"}
library(stringr)

# 移除停止词，也可尝试不移除或调整停止词列表以适应文档特性
cutter <- worker(bylines = TRUE, stop_word = "data/stop_word_zh.utf8")
# 分成文档，每个文档为一回（章）
by_chapter <- books %>%
  group_by(title) %>%
  mutate(text = sapply(segment(text, cutter), function(x){paste(x, collapse = " ")})) %>%
  mutate(chapter = cumsum(str_detect(text, "^第[零一二三四五六七八九十百 ]+回"))) %>%
  ungroup() %>%
  filter(chapter > 0) %>%
  unite(document, title, chapter)
# 按词切分
by_chapter_word <- by_chapter %>%
  unnest_tokens(word, text)
# 获得文档-词计数
word_counts <- by_chapter_word %>%
  count(document, word, sort = TRUE) %>%
  ungroup()
word_counts
```

### 按章得到 LDA

现在我们的 `word_counts` 数据框为 tidy 形式，每行每文档一个术语，但 topicmodels 包需要一个 `DocumentTermMatrix`。如在章 \@ref(cast-dtm) 中所叙述的，我们可以用 tidytext 的 `cast_dtm()` 把一个每行一个符号的表格映射为 `DocumentTermMatrix`。

```{r chapters_dtm}
chapters_dtm <- word_counts %>%
  cast_dtm(document, word, n)
chapters_dtm
```

然后即可使用 `LDA()` 函数创建一个六主题模型。在这个例子里我们知道要寻找6个主题，因为有6本书；在其它问题里可以尝试几个不同的 `k` 值以找到比较合适的。有意思的是，对于本数据集三主题也很准确，但那三本并不精彩的续书就要分别和其原著混起来了！作为练习，读者可以自己尝试。

```{r chapters_lda}
chapters_lda <- LDA(chapters_dtm, k = 6, control = list(seed = 1234))
chapters_lda
```

就像对 Associated Press 数据所做的一样，我们可以查看每主题每词概率。

```{r chapter_topics}
chapter_topics <- tidy(chapters_lda, matrix = "beta")
chapter_topics
```

注意，这将模型变成了每行术语每术语一个主题的格式。对于每个组合，模型计算了该术语来自该主题的概率。比如，术语“`r chapter_topics$term[2]`”在主题`r chapter_topics$topic[2]`中的概率有`r chapter_topics$beta[2]`。

我们可以用 dplyr 的 `top_n()` 得到每个主题概率最高的术语。

```{r top_terms}
top_terms <- chapter_topics %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)
top_terms
```

这个 tidy 输出可直接进行 ggplot2 可视化（图 \@ref(fig:toptermsplot)）。

```{r toptermsplot, fig.height=6, fig.width=7, fig.cap = "每个主题中最常见的术语"}
top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free") +
  coord_flip()
```

有些主题很清晰地与不同书联系了起来！“賈”属于《紅樓夢》，“行者”应该属于《西遊記》等。我们还可注意到，既然 LDA 属于“模糊聚类”方法，有些词为多个主题共有，如“道”“有”“人”等。

### 按文档分类 {#per-document}

此分析中的每个文档都表示了一个单章。因此，我们可能想要知道每个文档关联哪个主题。我们能把章放回到正确的书里吗？我们可以查看每主题每文档概率，$\gamma$（gamma）。

```{r chapters_gamma_raw}
chapters_gamma <- tidy(chapters_lda, matrix = "gamma")
chapters_gamma
```

每个值都是来自该文档的词属于该主题的估测比例。

现在有了主题概率，可以看看我们的无监督学习区分6本书的表现如何。我们预期一本书中的每章基本（或全部）出现在生成的对应主题中。

首先我们再把文档名拆分成书名与章回，然后绘制每主题每文档概率图 \@ref(fig:chaptersldagamma)。

```{r chapters_gamma}
chapters_gamma <- chapters_gamma %>%
  separate(document, c("title", "chapter"), sep = "_", convert = TRUE)
chapters_gamma
```

(ref:chaptersldagammacap) 每本书中每章的 $\gamma$ 概率

```{r chaptersldagamma, fig.width=8, fig.height=8, fig.cap = "(ref:chaptersldagammacap)"}
# reorder titles in order of topic 1, topic 2, etc before plotting
chapters_gamma %>%
  mutate(title = reorder(title, gamma * topic)) %>%
  ggplot(aes(factor(topic), gamma)) +
  geom_boxplot() +
  facet_wrap(~ title)
```

注意，看起来每本书的多数章都被识别为唯一的单个主题；同时，续作与原著有一定相似性。有没有哪个章节被识别到另一本书？首先我们用 `top_n()` 找到每章被关联最多的主题，即该章的“分类”。

```{r chapter_classifications, dependson = "chapters_gamma"}
chapter_classifications <- chapters_gamma %>%
  group_by(title, chapter) %>%
  top_n(1, gamma) %>%
  ungroup()
chapter_classifications
```

之后可以与每本书的“共识”主题进行比较（该书各章最多见的主题），看看哪些最常被误识别。

```{r book_topics, dependson = "chapter_classifications"}
book_topics <- chapter_classifications %>%
  count(title, topic) %>%
  group_by(title) %>%
  top_n(1, n) %>%
  ungroup() %>%
  transmute(consensus = title, topic)
chapter_classifications %>%
  inner_join(book_topics, by = "topic") %>%
  filter(title != consensus)
```

可以看到400章中只有几章识别错误，这个结果不算差。考虑到续作与原著的相似性，我们只挑出完全识别错的书。

```{r book_topics_filtered, dependson = "book_topics"}
library(stringi)

chapter_classifications %>%
  inner_join(book_topics, by = "topic") %>%
  filter(str_sub(title, -1, -1) != str_sub(consensus, -1, -1))
```

只剩两个章节，如《紅樓夢》第一回情节散漫，并无特定主题。

### 词的分配：`augment`

LDA 算法的一步是将每个文档里的每个词分配给一个主题。一个文档里越多词被分配给该主题，一般来说在文档-主题分类中的权重（`gamma`）就越高。

我们可能想要取出原始的文档-词对，找出每个文档里的哪些词分配到了哪个主题。这是 `augment()` 函数的功能，也来自 broom 包，也是 tidy 化模型输出的一部分。`tidy()` 获取模型的统计组件，`augment()` 使用模型在原始数据的每个观察上添加信息。

```{r assignments, dependson = "chapters_lda"}
assignments <- augment(chapters_lda, data = chapters_dtm)
assignments
```

这返回书-术语计数的 tidy 数据框，但额外加了一列：`.topic`，每个文档中每个术语所分配的主题。（`augment` 添加的额外列总以 `.` 开头，以避免覆盖已有的列）。我们可以把共识书名与此 `assignments` 表格结合，找出未正确分类的词。

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

这个真实书名（`title`）和被分配的书名（`consensus`）的组合对进一步探索很有用。例如我们可以绘制 **混淆矩阵**，显示一本书里的词有多经常被分配给另一本，使用 dplyr 的 `count()` 和 ggplot2 的 `geom_tile`（图 \@ref(fig:confusionmatrix)。

```{r confusionmatrix, dependson = "assignments2", fig.width = 10, fig.height = 8, fig.cap = "展示 LDA 对每本书中的词进行分配的混淆矩阵。表格的每行代表词来自原书，每列代表被分配的书。"}
assignments %>%
  count(title, consensus, wt = count) %>%
  group_by(title) %>%
  mutate(percent = n / sum(n)) %>%
  ggplot(aes(consensus, title, fill = percent)) +
  geom_tile() +
  scale_fill_gradient2(high = "red", label = percent_format()) +
  theme_minimal() +
  labs(x = "词被分配到的书",
       y = "词来自的书",
       fill = "正确分配的百分比") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        panel.grid = element_blank(),
        text = element_text(family = "wqy-microhei"))
```

注意到几乎所有词都正确分配了，但不少书里的词都倾向于了《水滸後傳》。最常出错的词有哪些？

```{r wrong_words, dependson = "assignments2"}
library(stringi)

wrong_words <- assignments %>%
  filter(str_sub(title, -1, -1) != str_sub(consensus, -1, -1))
wrong_words %>%
  count(title, consensus, term, wt = count) %>%
  ungroup() %>%
  arrange(desc(n))
```

这里我们不看原著与续作互相混淆的词。長老、和尚、行者等称呼在多部书中都出现不少，错误分配可以理解。

另外，有可能有的词在一本书中从未出现过但仍然被分配给该书。如“礙於”只见于《後西游記》一次，却被分配给了《水滸後傳》。

```{r dependson = "word_counts"}
word_counts %>%
  filter(word == "礙於")
```

LDA 算法有一定随意性，有可能把一个主题覆盖到多本书。

## 另一种 LDA 实现

topicmodels 包中的 `LDA()` 函数只是隐含狄利克雷分布分配算法的一种实现。例如 [mallet](https://cran.r-project.org/package=mallet) 包 [@R-mallet] 就实现了对 [MALLET](http://mallet.cs.umass.edu/) Java 文本分类工具包的包装，tidytext 包也提供了对这种模型的输出的 tidy 化工具。如果无法正常调用 `library(mallet)` （主要原因是其 rJava 依赖），可以尝试在终端中运行 `R CMD javareconf`，成功后重新打开 RStudio。

mallet 包采用了不太一样的输入格式。比如输入未经符号化的文档并自行符号化，且需要一个单独的停止词文件。这意味着我们在运行 LDA 之前需要把每个文档的文本都合并成一个字符串。

```{r mallet_model, results = "hide", eval = FALSE}
library(mallet)
# create a vector with one string per chapter
collapsed <- by_chapter_word %>%
  anti_join(stop_words, by = "word") %>%
  mutate(word = str_replace(word, "'", "")) %>%
  group_by(document) %>%
  summarize(text = paste(word, collapse = " "))
# create an empty file of "stopwords"
file.create(empty_file <- tempfile())
docs <- mallet.import(collapsed$document, collapsed$text, empty_file)
mallet_model <- MalletLDA(num.topics = 4)
mallet_model$loadDocuments(docs)
mallet_model$train(100)
```

然而，模型创建后我们就可以像本章其它部分一样以几乎一致的方式使用 `tidy()` 和 `augment()` 函数。包括提取每个主题中词或每个文档中主题的概率。

```{r eval = FALSE}
# word-topic pairs
tidy(mallet_model)
# document-topic pairs
tidy(mallet_model, matrix = "gamma")
# column needs to be named "term" for "augment"
term_counts <- rename(word_counts, term = word)
augment(mallet_model, term_counts)
```

可以以同样的方式用 ggplot2 探索和绘制 LDA 输出的模型。

## 小结

本章介绍了主题模型，用于寻找对一组文档具有代表性的词的聚类，展示了 `tidy()` 功能如何帮我们使用 dplyr 和 ggplot2 探索和理解这些模型。这是用 tidy 方法探索模型的一个优势：输出格式不同的挑战可由 tidy 化的函数处理，我们可以用一组标准工具探索模型的结果。特别地，我们看到主题模型能够把来自4本分立的书混合的章节拆分并识别，也探索了模型的局限性，找到了错误分配的词等。