rm(list=ls(all=TRUE))
getwd()
# with setwd() you change the current working directory of the R process
setwd("write here your working directory")
getwd()
library(quanteda)
library(readtext)
library(stm)
################################################################
# Example with U.S. Presidential Inaugural Speeches since 1970s
################################################################
summary(data_corpus_inaugural)
# we could add some document-level variables – what quanteda calls docvars – to this corpus.
library(xlsx)
meta.pr <- read.xlsx("C:/Users/mw/Dropbox (VOICES)/TOPIC MODEL/List President USA party.xlsx", 1)
str(meta.pr)
levels(meta.pr$Party)
# do not include factor variable with an unobserved level in the dataset you want to analyze!!!
# for example in our case this happens with "Democratic-Republican" "Federalist" "Unaffiliated" & "Whig"
levels(meta.pr$Party) <- list(Democratic="Dem", Republican="Rep")
levels(meta.pr$Party)
docvars(data_corpus_inaugural, "Party") <- meta.pr$Party
summary(data_corpus_inaugural)
# inspect the document-level variables (that's very important as an option, as we will see..)
head(docvars(data_corpus_inaugural))
myCorpus <- corpus_subset(data_corpus_inaugural, Year > 1970)
summary(myCorpus)
str(myCorpus)
myDfm <- dfm(myCorpus, remove = stopwords("english"), tolower = TRUE, stem = TRUE,
remove_punct = TRUE, remove_numbers=TRUE)
topfeatures(myDfm , 20) # 20 top words
# if you want, you can decide to keep only words occuring >=10 times and in >=2 docs. This is highly suggested
# with many documents.
myDfm.trim <- dfm_trim(myDfm, min_count = 10, min_docfreq = 2)
# Or just drop the words that occurred either too frequently (in more than 99% of documents,
# or too infrequently (in less than 1% of the documents)
myDfm.trim2 <- dfm_trim(myDfm, max_docfreq = 0.99, min_docfreq = 0.01)
topfeatures(myDfm.trim2, 20) # 20 top words
# convert from quanteda to stm: this is a crucial step!!!!
Dfmstm2 <- convert(myDfm, to = "stm", docvars = docvars(myCorpus))
str(Dfmstm2)
# STM run: 1) we want to extract 5 topics; 2) we want to control for the following covariates in affecting topic prevalence: Party and Year;
# we could have ran also at the same time a model with covariates affecting topical content.
# As currently implemented this must be a single variable which defines a discrete partition of the dataset
# (each document is in one and only one group). While including more covariates in topical prevalence will rarely affect
# the speed of the model, including additional levels of the content covariates can make the model much slower to converge
# This is due to the model operating in the much higher dimensional space of words in dictionary
# (which tend to be in the thousands) as opposed to topics. That's why I suggest you to run first a model w/o covariates
# for topical content; 3) you can also run a model w/o covariates. In this case you have an unstructured
# topic model, by fitting a Correlated Topic Model; 4) highly suggested to employ the "spectral" intitialization
# based on the method of moments, which is deterministic and globally consistent under reasonable conditions.
# This means that no matter the seed that is set, the SAME results will be generated. As an alternative you could employ
# a LDA initialization. In this, case, however, the answers the estimation procedure comes up with may depend on starting
# values of the parameters (e.g., the distribution over words for a particular topic). So you should also define
# a set seed to replicate your analysis!
stmFitted <- stm(Dfmstm2$documents, Dfmstm2$vocab, K = 5, max.em.its = 75,
prevalence = ~ Party +Year, data = Dfmstm2$meta, init.type = "Spectral")
# Interpreting the STM by plotting and inspecting results
labelTopics(stmFitted, c(1:5))
# Summary visualization: understanding the extracted topics through words and example documents
# First: which words/topic association?
labelTopics(stmFitted)
plot(stmFitted)
plot(stmFitted, type = "summary", labeltype = c("frex"))
plot(stmFitted, type = "labels", labeltype = c("frex"))
cloud(stmFitted, topic = 3)
# frequency distribution of topics across the 12 documents
plot(stmFitted, type = "hist", labeltype = c("frex"))
# cloud(poliblogPrevFit, topic = 7, scale = c(2,.25))
# Let's read the documents most associated with each topic
docs <- myCorpus$documents$texts
str(docs)
# you cannot read it properly!
thought <- findThoughts(stmFitted, texts=docs, topics=2, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought, width = 30, main = "Topic 2")
# save only the first 200 words from each document and display them!
docs2 <- substr(docs, start = 1, stop = 200)
thought2 <- findThoughts(stmFitted, texts=docs2, topics=2, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought2, width = 30, main = "Topic 2")
thought5 <- findThoughts(stmFitted, texts=docs2, topics=5, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought5, width = 30, main = "Topic 5")
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought2, width = 30, main = "Topic 2")
plotQuote(thought5, width = 30, main = "Topic 5")
# but to whom those speeches belong to?
str(stmFitted$theta)
theta <-as.data.frame(stmFitted$theta)
str(theta )
combine <- paste(myCorpus$documents$President, myCorpus$documents$Year, sep = " ")
combine
rownames(theta ) <- combine
# which is the document with the highest value of Topic 2?
rownames(theta )[which.max(theta $V2)]
# which are the documents with the highest value for each Topic?
rownames(theta )[apply(theta , 2, which.max)]
# and the 2nd amd 3rd largest one?
maxn <- function(n) function(x) order(x, decreasing = TRUE)[n]
rownames(theta )[apply(theta , 2, maxn(1))]
rownames(theta )[apply(theta , 2, maxn(2))]
rownames(theta )[apply(theta , 2, maxn(3))]
# I want to read the Obama speech of 2009 to better understand the content of Topic 4
Obama <- texts(myCorpus)[10]
strwrap(Obama [[1]])
# Plotting the graph to see the distribution across documents of Topic 1
plot(theta$V1,
main= "Proportion of Topic 1 across documents",
xlab= "Documents",
ylab= "% Topic 1 in the document",
col= "blue", pch = 19, cex = 1, lty = "solid", lwd = 2)
text(theta$V1, labels=combine, cex= 1, pos=3)
# Plotting the graph to see the distribution of Topic 1 and 2 across documents
plot(theta$V1, theta$V2,
main= "Proportion of Topic 1 and 2 across documents",
xlab= "% Topic 1 in the document",
ylab= "% Topic 2 in the document",
col= "blue", pch = 19, cex = 1, lty = "solid", lwd = 2)
text(theta$V1, theta$V2, labels=combine, cex= 1, pos=3)
# Additional tools for interpretation and visualization
# Uses a topic correlation graph estimated by topicCorr and the igraph package to plot a network where nodes are topics
# and edges indicate a positive correlation. If no edges = only negative correlations!
library(igraph)
mod.out.corr <- topicCorr(stmFitted)
plot(mod.out.corr)
mod.out.corr$cor
# http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram
library(corrplot)
corrplot(mod.out.corr$cor, method="circle")
corrplot(mod.out.corr$cor, method="circle", type="lower")
corrplot(mod.out.corr$cor, method="circle", type="upper")
corrplot(mod.out.corr$cor, method="pie")
corrplot(mod.out.corr$cor, method="color")
corrplot(mod.out.corr$cor, method="number")
# Estimating topical prevalence
prep <- estimateEffect(1:5 ~ Party+ Year,
stmFitted, meta = Dfmstm2$meta, uncertainty = "Global")
str(prep)
# Method used for plotting topical prevalence:
# a) "pointestimate" estimates mean topic proportions for each value of the covariate.
# b) "difference" estimates the mean difference in topic proportions for two different values of the covariate (cov.value1 and cov.value2 must be specified).
# c) "continuous" estimates how topic proportions vary over the support of a continuous covariate.
plot(prep, "Year", method = "continuous", topics = 5,
model = stmFitted, printlegend = FALSE, xaxt = "n", xlab = "Left-Right")
seq <- seq(from = as.numeric("1970"), to = as.numeric("2017"))
axis(1, at = seq)
title("Topic 5")
abline(h=0, col="blue")
plot(prep, "Party", model=stmFitted,
method="pointestimate")
plot(prep, covariate = "Party", topics = c(1, 2, 3, 4, 5),
model = stmFitted, method = "difference",
cov.value1 = "Democratic", cov.value2 = "Republican", xlim = c(-1, 1),
xlab = "More Repubblican... More Democratic",
main = "Effect of Democratic vs. Repubblican",
labeltype = "custom",
custom.labels = c('Topic 1', 'Topic 2', 'Topic 3', 'Topic 4', 'Topic 5'))
# Plotting covariate interactions
stmFitted_interaction <- stm(Dfmstm2$documents, Dfmstm2$vocab, K = 5, max.em.its = 75,
prevalence = ~ Party * Year, data = Dfmstm2$meta, init.type = "Spectral")
prep <- estimateEffect(c(5) ~ Party * Year,
stmFitted_interaction, meta = Dfmstm2$meta, uncertainty = "Global")
# This plots the interaction between Year and Party (Democratic versus Repubblican). Topic 5 prevalence is plotted as linear
# function of year, holding the rating at either 0 (Democratic) or 1 (Republican). Were other
# variables included in the model, they would be held at their sample medians.
plot(prep, covariate = "Year", model = stmFitted_interaction,
method = "continuous", xlab = "Year", moderator = "Party",
moderator.value = "Democratic", linecol = "blue", ylim = c(-1, 1),
printlegend = F)
plot(prep, covariate = "Year", model = stmFitted_interaction,
method = "continuous", xlab = "Year", moderator = "Party",
moderator.value = "Republican", linecol = "red", add = T,
printlegend = F)
legend(2000, -0.7, c("Democratic", "Republican"),
lwd = 2, col = c("blue", "red"))
abline(h=0, col="black")
# Topical content
stmFitted_content <- stm(Dfmstm2$documents, Dfmstm2$vocab, K = 5, max.em.its = 75,
prevalence = ~ Party +Year, content = ~ Party, data = Dfmstm2$meta, init.type = "Spectral")
# This figure shows how democratic and repubblican Presidents talk about topic 5 differently
plot(stmFitted_content, type = "perspectives", topics = 5)
# This function can also be used to plot the contrast in words across two topics
# This plot calculates the difference in probability of a word for the two topics, normalized by the maximum
# difference in probability of any word between the two topics.
plot(stmFitted, type = "perspectives", topics = c(2, 5))
# plot(stmFitted, type = "perspectives", labeltype = c("frex"), topics = c(2, 5))
# Model selection
# Model search across numbers of topics STM assumes a fixed user-specified number of topics.
# There is not a right answer to the number of topics that are appropriate for a given corpus, but
# the function searchK uses a data-driven approach to selecting the number of topics.
# The function will perform several automated tests to help choose the number of topics including
# calculating frequency, coherence, the held out likelihood and performing a residual analysis. Still nothing replaces your interpretation!!!!
# For example, one could estimate a STM model for 5, 7, 10 and 20 topics and compare the results along each of the criteria
K <-c(5,7,10, 20)
# you can also relax the number of maximum required iterations if you want, but that would require you (much) more
# times especially with large datasets
storage <- searchK(Dfmstm2$documents, Dfmstm2$vocab, K, max.em.its = 75,
prevalence = ~ Party +Year, data = Dfmstm2$meta, init.type = "Spectral")
str(storage)
# Often times users will select a model with desirable properties in both dimensions
# (i.e., models with average scores towards the upper right side of the plot)
plot(storage$results$semcoh, storage$results$exclus,
xlab= "Semantic coherence",
ylab= "Exclusivity",
col= "blue", pch = 19, cex = 1, lty = "solid", lwd = 2)
text(storage$results$semcoh, storage$results$exclus, labels=storage$results$K, cex= 1, pos=4)
plot(storage)
# When init.type="Spectral" and K=0 the number of topics is set using the algorithm in Lee and Mimno (2014).
# This method does not estimate the "true" number of topics and does not necessarily have any particular statistical
# properties for consistently estimating the number of topics. It can however provide a useful starting point.
# The added randomness from the projection means that the algorithm is not deterministic
# as with the standard "Spectral" initialization type. Running it with a different seed can result in not only different
# results but a different number of topics.
total <- stm(Dfmstm2$documents, Dfmstm2$vocab, K=0, max.em.its = 500,
prevalence = ~ Party +Year, data = Dfmstm2$meta, init.type = "Spectral", seed = 8458159)
labelTopics(total)
mod.out.corr2 <- topicCorr(total)
plot(mod.out.corr2)
mod.out.corr2$cor
corrplot(mod.out.corr2$cor, method="pie", type="lower")
# In addition to these functions one can also explore if there words that are extremely highly associated with
# a single topic via the checkBeta function. The ouput gives the user lists of which topics have problems, which words in which topics have problems,
# as well as a count of the total problems in topics and the total number of problem words.
# In our case, we do not have such kind of problem
check <- checkBeta(stmFitted, tolerance = 0.01)
str(check)
##### graphic extension
# Generate STM Correlation Tree
# This function generates an interactive, full-model HTML visualization of topic hierachies (using hierarchical clustering) for a fitted STM model.
# The visualization highlights the correlations among topics, and can be used to view the model at differing levels of complexity
# The visualization takes the form of an indented tree. The leaves of the tree correspond to topics. The leaf nodes are grouped in topic clusters.
# This allows the model to be visualized at differing levels of aggregation
# similar results of plot(mod.out.corr)
library(stmCorrViz)
stmCorrViz(stmFitted, "corrviz.html")
labelTopics(stmFitted)
# to show also the documents example for each topic
stmCorrViz(stmFitted, "corrviz.html", documents_raw=myCorpus$documents$texts, documents_matrix=Dfmstm2$documents)
# with the topics that you get when K=0
stmCorrViz(total, "corrviz.html", documents_raw=myCorpus$documents$texts, documents_matrix=Dfmstm2$documents)
docs <- myCorpus$documents$texts
str(docs)
docs2 <- substr(docs, start = 1, stop = 200)
str(docs2)
thought <- findThoughts(stmFitted_interaction, texts=docs2, topics=1, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought, width = 30, main = "Topic 1")
##### Outputs json file and creates visualization from stm
# This function outputs necessary json files to a working directory, then opens a browser with the visualization
# of the relationship between covariates and topics in the stm.
# better call this library after a STM with only topic prevalence (and not topical content)
library(stmBrowser)
Dfmstm2$meta$docs <- myCorpus$documents$texts
str(Dfmstm2$meta)
stmBrowser(stmFitted, data= Dfmstm2$meta, c("Year", "President", "FirstName", "Party"),
text="docs", labeltype = "frex")
# select in the HTML file you produced "X-axis=Topic 4", "Y-axis=Year", "Radius=non", "Color=party";
# by doing that you can clearly see that topic 4 is more popular on average among Democratic Presidents
# [exactly as found in the STM analysis below!]
plot(prep, covariate = "Party", topics = c(1, 2, 3, 4, 5),
model = stmFitted, method = "difference",
cov.value1 = "Democratic", cov.value2 = "Republican", xlim = c(-1, 1),
xlab = "More Repubblican... More Democratic",
main = "Effect of Democratic vs. Repubblican",
labeltype = "custom",
custom.labels = c('Topic 1', 'Topic 2', 'Topic 3', 'Topic 4', 'Topic 5'))
# if you select "Year" on the Y axis you see however only some numbers, because "Year" in the meta data is a numeric not a date.
# Let's then create a date out of the numeric variable of Year (and let's call it "Year2"), and let's use it as well!
Dfmstm2$meta$Year2 <- Dfmstm2$meta$Year
Dfmstm2$meta$Year2 <- as.character(Dfmstm2$meta$Year2)
str(Dfmstm2$meta)
paste(Dfmstm2$meta$Year2,"-01","-01",sep="")
Dfmstm2$meta$Year2<- as.Date(paste(Dfmstm2$meta$Year2,"-01", "-01",sep=""))
str(Dfmstm2$meta)
# now if you select "Year2" on the Y axis you see the year not numbers! Good work!
stmBrowser(stmFitted, data= Dfmstm2$meta, c("Year2", "Year", "Party"),
text="docs", labeltype = "frex")
##### A third possible graphic visualization
# to learn more about this last visualization, take a look at: https://rdrr.io/cran/stm/man/toLDAvis.html
library(LDAvis)
library(servr)
toLDAvis(mod=stmFitted, docs=Dfmstm2$documents)
# topics displayed according to their frequency
# focus on topic 4 is 5 in toLADavis