rm(list=ls(all=TRUE))
setwd("C:/Users/luigi/Dropbox/TOPIC MODEL")
getwd()
library(quanteda)
library(readtext)
library(ggplot2)
library(stm)
library(igraph)
#########################################
#########################################
################################################################
# Let's estimate a STM
################################################################
#########################################
#########################################
# We will now use a dataset that contains the lead paragraph of around 5,000 articles about the economy published
# in the New York Times between 1980 and 2014.
nyt <- readtext("nyt.csv", text_field = "lead_paragraph")
str(nyt)
nyt$datetime <- as.Date(nyt$datetime, "%m/%d/%Y") # let's convert the string of the date in a date
str(nyt)
# let's create a new variable according to the party of the president in charge during the temporal period considered
year <- as.numeric(substr(nyt$datetime, 1, 4)) # let's extract the year [i.e., the four first numbers that appear]
table(year)
repub <- ifelse(year %in% c(1981:1992, 2000:2008), 1, 0) # 1 for repub; 0 for dems
table(repub)
nyt$repub <- repub
nyt$year <- year
str(nyt)
library(dplyr)
# let's extract a random sample of 1,000 texts (to make things faster in the Lab!)
set.seed(123)
nyt2 <- sample_n(nyt, 1000)
str(nyt2)
myCorpus <- corpus(nyt2)
head(summary(myCorpus))
tok2 <- tokens(myCorpus , remove_punct = TRUE, remove_numbers=TRUE, remove_symbols = TRUE, split_hyphens = TRUE, remove_separators = TRUE)
tok2 <- tokens_remove(tok2, stopwords("en"))
tok2 <- tokens_wordstem (tok2)
myDfm <- dfm(tok2)
# Some form of trim is usually suggested when you have many documents and you want to run a STM (exactly as it happens with any topic-models)
# In our case, given the relative low number of documents and features, we are going to proceed with our original dfm (i.e., myDfm)
# REMEMBER: if you trim the dfm, be sure that each of the documents (i.e., rows) in the Dfm does not report just 0s!
myDfm [ntoken(myDfm ) == 0,]
# everything is fine with our example, otherwise drop delete such documents (i.e., rows) from the dfm before passing them to "stm".
# myDfm <- news_dfm[ntoken(myDfm) > 0,]
# Convert the dfm from Quanteda to STM: this is a crucial step.
# When coverting the dfm, you have also to list the document variables that are present in the corpus
# in our case: datetime, repub, year
head(docvars(myCorpus))
DfmStm <- convert(myDfm, to = "stm", docvars = docvars(myCorpus))
str(DfmStm$meta)
################################################################
# STM run
################################################################
# 1) we want to extract 15 topics (i.e., K=15)
# 2) we want to control for the following covariates in affecting topic prevalence: repub and year.
# We therefore need to list such variables in the comman below via "prevalence = ~ repub + s(year)"
# Note that by writing "s(year)" we use a spline functions for non-linear transformations of the time-variable.
# if you want to add just a linear relationship between "year" and topics, just write "year"
# We could have ran a model with covariates affecting topical content as well (see below for an example).
# While including more covariates in topic prevalence will rarely affect the speed of the model,
# including additional levels of the content covariates (topical content) 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) it is generally advisable 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 of 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 results you get out of the estimation procedure may depend
# on starting values of the parameters (e.g., the distribution over words for a particular topic) - as it is always the case with a LDA (remember?).
# So remember always to define a set seed to replicate your analysis in this latter case
# You can also relax the number of maximum required iterations if you want, but that would require you (much) more
# time especially with large datasets. However 75 is a very small number! We apply that here just to save time!
# The default value is 500.
system.time(stmFitted <- stm(DfmStm $documents, DfmStm $vocab, K = 15, max.em.its = 75,
prevalence = ~ repub + s(year), data = DfmStm $meta, init.type = "Spectral")) # around 36 seconds on my laptop
##############################################################################################################
# NB: running a model with or without covariates for topic prevalance produces of course different results;
# the same is true if you run a model with covariates just for the topic prevalance part, or with both the
# topic prevalance part as well as the topical content part
##############################################################################################################
################################################################
# Interpreting the STM by plotting and inspecting results
################################################################
# Summary visualization: understanding the extracted topics through words and example documents
# First: which words/topic association?
labelTopics(stmFitted, n=7) # 7 features for each topic (including FREX words!)
# The frequency/exclusivity (FREX) scoring summarizes words according to their probability of appearance under a topic and the exclusivity
# to that topic. More in details: the FREX for a word is computed as the harmonic mean of the word's rank in terms of exclusivity and frequency.
# [Harmonic mean? The harmonic mean can be expressed as the reciprocal of the arithmetic mean of the reciprocals of the given set of observations. Say you have 2, 3 and 4,
# the harmonic mean is 3/(1/2+1/3+1/4)]
# The harmonic mean is attractive here because it does not allow a high rank along one of the dimensions to compensate for the lower rank in another (as it happens with an
# arithmetic mean. Thus words with a high score must be high along both dimensions (exclusivity AND frequency).
# As a result, FREX words provide quite often more semantically intuitive representations of a given topic.
labelTopics(stmFitted, n=7, topics=1) # the same just for topic 1
plot(stmFitted, type = "labels", labeltype = c("frex"), n=5) # plot just frex words
plot(stmFitted, type = "summary", labeltype = c("frex"), n=5) # topic 6 is the most frequent one
# topic meaning: according to frex words, topic 6 seems to be related to the trend in the economy; topic 10 to politics;
# topic 11 to inflation; etc.
plot(stmFitted, type = "hist", labeltype = c("frex")) # Here topic 6 appears as more "evenly" distributed across documents than
# for example topic 10 for example
# Let's read the documents most associated with each topic
# In particular, let's identify the most representative documents for a particular topic.
# We can also use this in order to get a better sense of the content of actual documents with a high topical content.
docs <- nyt2$text
str(docs)
# Let's focus on topic 6 for example and let's identify the 3 texts with the highest theta for that topic
# and let's read the first 200 words from each of them
docs2 <- substr(docs, start = 1, stop = 200)
thought6 <- findThoughts(stmFitted, texts=docs2, topics=6, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought6, width = 30, main = "Topic 6 - Economy")
# Let's focus on topic 10 and let's do the same
thought10 <- findThoughts(stmFitted, texts=docs2, topics=10, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought10, width = 30, main = "Topic 10 - Politics")
# Let's plot Topic 6 and 10 together
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought6, width = 30, main = "Topic 6 - Economy")
plotQuote(thought10, width = 30, main = "Topic 10 - Politics")
# which are the documents with the highest theta for each k?
# Note: when you write "2", the apply function works over columns; "1" over rows. We want here to identify the highest
# theta over columns (i.e., the highest theta for each single topic extracted)
apply(stmFitted$theta,2,which.max)
# let's read the entire documents with the highest theta for topic=6
strwrap(as.character(myCorpus)[93])
# same one as below
thought6 <- findThoughts(stmFitted, texts=docs2, topics=6, n=1)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought6, width = 30, main = "Topic 6")
# which is the document with the 2nd largest value for each topic?
maxn <- function(n) function(x) order(x, decreasing = TRUE)[n]
apply(stmFitted$theta, 2, maxn(2))
# let's read the entire documents with the second highest theta for topic=6
strwrap(as.character(myCorpus)[67])
# which is the document with the 3rd largest value for each topic?
apply(stmFitted$theta, 2, maxn(3))
# which are the most likely topics across our documents?
# NOTE: here we employ the apply function over rows (i.e., over documents)
apply(stmFitted$theta,1,which.max)
table(apply(stmFitted$theta,1,which.max) )
# you can save them as a document-level variable.
docvars(myDfm , 'STMtopic') <- apply(stmFitted$theta,1,which.max)
str(myDfm)
head(docvars(myDfm))
# or you can also save them back in the original dataframe
nyt2$topic <- apply(stmFitted$theta,1,which.max)
str(nyt2)
# Topic 6 - 5 random documents associated to it
set.seed(123)
sample(nyt2$text[nyt2$topic==6], 5)
# But we already know that STM actually computes a distribution over topics.
# In other words, each document is considered to be about a mixture of topics as we have already discussed!
# This information is included in the matrix thera in the STM object
# For example, news 1 is 44% about topic 9 and 34% about topic 7 for example
round(stmFitted$theta[1,], 2)
################################################################
# Checking for topic correlation
################################################################
# 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!
mod.out.corr <- topicCorr(stmFitted)
plot(mod.out.corr) # in our case topics 8, 10 and 15 are positively correlated
plot(stmFitted, type = "summary", labeltype = c("frex"), n=5)
# i.e. when an article discusses about Reagan (topic 8) there is a positive
# probability of discussing about public budget (topic 15) and politics (topic 10)
#########################################
# Estimating topic prevalence using as IV: repub and year
#########################################
# The estimateEffect function performs a regression where topic-proportions are the outcome variable.
# The function can handle factors and numeric variables. Dates should be converted to numeric variables before analysis.
# When you estimate the effects, you need always to seet a seed
set.seed(123)
prep <- estimateEffect(1:15 ~ repub + s(year), stmFitted, meta = DfmStm $meta)
set.seed(123)
summary(prep)
# for example there is a tendency to talk more about topic 15 under a Republican Presidency
# Method used for plotting topic prevalence:
# a) "continuous" estimates how topic proportions vary over the support of a continuous covariate.
# b) "pointestimate" estimates mean topic proportions for each value of the categorical covariate
# c) "difference" estimates the mean difference in topic proportions for two different values of the covariate (cov.value1 and cov.value2 must be specified)
# Let's start with year (i.e., a continuum variable):
# let's plot it! as you can see topic 10 prevalence (politics) is stable but at the end it increases (after the 2008
# economic crisis)
plot(prep, "year", method = "continuous", topics = 10,
model = stmFitted, printlegend = FALSE, xaxt = "n", xlab = "Year")
seq <- seq(from = as.numeric("1980"), to = as.numeric("2014"))
axis(1, at = seq)
title("Topic 10 - Politics")
abline(h=0, col="blue")
# the opposite happens with topics 11 (inflation) - although the decrease does not appear to be significant
plot(prep, "year", method = "continuous", topics = 11,
model = stmFitted, printlegend = FALSE, xaxt = "n", xlab = "Year")
seq <- seq(from = as.numeric("1980"), to = as.numeric("2014"))
axis(1, at = seq)
title("Topic 11 - Inflation")
abline(h=0, col="blue")
# Let's now move to repub (i.e., a categorical variable) for topic 11 (inflation) and topic 15 (budget issues)
plot(stmFitted, type = "summary", labeltype = c("frex"), n=5)
# the salience to topic 11 seems to be higher under a Democratic presidency. While there is a tendency to discuss more
# about topic 15 under a Republican presidency
plot(prep, "repub", model=stmFitted, method="pointestimate", topic=c(11, 15))
# let's see if the difference between the two coefficients is statistically different from 0.
# Always remember that you plot "coef.value1-coef.value2". In this case (see below the script)
# it implies "coeff. of republican-coeff. of democratic" given that republic is a 0/1 variable and cov.value1 = "1" (i.e., Republican)
# and cov.value2 = "0" (i.e., Democratic)
plot(prep, covariate = "repub", topics = c(11, 15),
model = stmFitted, method = "difference",
cov.value1 = "1", cov.value2 = "0", xlim = c(-1, 1),
xlab = "More Democratic... More Repubblican",
main = "Effect of Democratic vs. Repubblican Presidency",
labeltype = "custom",
custom.labels = c('Topic 11 - Inflation', 'Topic 15 - Public Budget'))
# and what if we plot all the topics together?
a=c()
for (i in 1:15) {
d=paste("Topic ", i)
a=c(a,d)
}
a
plot(prep, covariate = "repub", topics = c(1:15),
model = stmFitted, method = "difference",
cov.value1 = "1", cov.value2 = "0", xlim = c(-1, 1),
xlab = "More Democratic... More Repubblican",
main = "Effect of Democratic vs. Repubblican Presidency",
labeltype = "custom",
custom.labels = a)
#########################################
# Estimating topical content using as IV: repub
#########################################
# Let's run a topical content analysis by including as a covariate repub.
# As currently implemented the variable used for topical content must be a single variable which defines a discrete partition of the dataset
# (each document is in one and only one group - so only categorical variables can be employed)
# You will need around 4 mins to run it. So let's call the RData where the model is already fitted
# system.time(stmFitted_topical <- stm(DfmStm $documents, DfmStm $vocab, K = 15, max.em.its = 75,
# prevalence = ~ repub + s(year), content = ~ repub, data = DfmStm $meta, init.type = "Spectral")) # around 191 seconds [3 mins]
load("C:\\Users\\luigi\\Dropbox\\TOPIC MODEL\\Stm.RData")
ls()
# You see that the results you get with and w/o topical content are somehow different
table(apply(stmFitted$theta,1,which.max) )
table(apply(stmFitted_topical$theta,1,which.max) )
apply(stmFitted$theta,2,mean) # mean distribution of topics
apply(stmFitted_topical$theta,2,mean) # mean distribution of topics
# high - but not perfect - correlation
cor(apply(stmFitted$theta,2,mean), apply(stmFitted_topical$theta,2,mean))
# This figure shows how the nyt articles talk about topic 11 (inflation) and topic 10 (politics) under a democratic or repubblican presidency
plot(stmFitted_topical, type = "perspectives", topics = 11)
plot(stmFitted_topical, type = "perspectives", topics = 10)
# 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 selected two topics, normalized by the maximum
# difference in probability of any word between the two topics (here topics 11 and 10)
plot(stmFitted_topical, type = "perspectives", topics = c(11, 10))
#########################################
# 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.
# For example, one could estimate a STM model for 3, 4 and 5 topics (just as an example!)
# and compare the results along each of the criteria
# But REMEMBER: nothing replaces your interpretation!!!!
# Add always a set.seed here (despite you employ a Spectral model, given that in this case there is a probabilistic part of the analysis. Why?
# Cause the function also calculate the heldout-likelihood - and therefore documents have proportion of them heldout at random
set.seed(02138)
K <-c(3,4,5)
system.time(storage <- searchK(DfmStm $documents, DfmStm $vocab, K, max.em.its = 75,
prevalence = ~ repub + s(year), data = DfmStm $meta, init.type = "Spectral")) # around 1 minute
#####################################################
# The semantic coherence-exclusivity ‘frontier’
#####################################################
# You should focus, as we did in the topic model class, on those models that lie on the semantic coherence-exclusivity ‘frontier’,
# that is, where no model strictly dominates another in terms of semantic coherence and
# exclusivity (i.e., models with average scores towards the upper right side of the plot).
# The figure below identifies model 5 as the one that better satisfies the
# ‘frontier’ criterion, i.e., a model with desirable properties in both dimensions
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)
# 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.
# Do not run in the Lab!!!
# library(geometry)
# library(Rtsne)
# library(rsvd)
# set.seed(123)
# system.time(total <- stm(DfmStm $documents, DfmStm $vocab, K=0, max.em.its = 500,
# prevalence = ~ Party +Year, data = DfmStm $meta, init.type = "Spectral"))
# labelTopics(total)
#########################################
##### graphic extensions about STM. There are several of them!
##### Take a look at library(stmCorrViz) and library(stminsights) among the others
#########################################
# Let's take a look at stmBrowser
library(stmBrowser)
str(nyt2)
DfmStm $meta$docs <- nyt2$text
str(DfmStm $meta)
stmBrowser(stmFitted, data= DfmStm $meta, c("repub", "year"), text="docs", labeltype = "frex")
# you have an error! Why? cause in the original nyt2$text you have some characters that are not valid UTF-8
# This creates problem when transforming the object to JSON (using the rjson package).
# Therefore: let's use iconv
nyt2$text2 <- iconv(nyt2$text, "UTF-8", "UTF-8",sub='')
str(nyt2)
DfmStm $meta$docs <- nyt2$text2
str(DfmStm $meta)
stmBrowser(stmFitted, data= DfmStm $meta, c("repub", "datetime"), text="docs", labeltype = "frex")
# as an example: select in the HTML file you produced "X-axis=Topic 15" (public budget), "Y-axis=datetime", "Radius=none", "Color=repub"
# blue=Democratic Presidency; red=Republican Presidency