rm(list=ls(all=TRUE))
setwd("C:/Users/luigi/Dropbox/TOPIC MODEL/FOR NEXT YEAR 2020/Newsmap")
getwd()
library(quanteda)
library(readtext)
library(ggplot2)
library(ldatuning)
library(topicmodels)
library(lubridate)
library(topicdoc)
library(cowplot)
library(stm)
library(igraph)
#########################################
#########################################
################################################################
########## Let's estimate a Topic-model (LDA)
################################################################
#########################################
#########################################
# This data frame contains 1,959 Guardian news articles published in 2016
myText <- read.csv("guardian.csv", stringsAsFactors=FALSE)
str(myText)
# creating a new variable "text2". Later it will be useful
myText$text2 <- myText$text
news_corp <- corpus(myText)
head(summary(news_corp))
# After removal of function words and punctuations in dfm(), we keep only the top 5% of the most frequent features
# (min_termfreq = 0.95) that appear less than 10% in all the documents (max_docfreq = 0.1) using dfm_trim()
# to focus on common but distinguishing features.
news_dfm <- dfm(news_corp, remove =(stopwords("english")), tolower = TRUE, stem = TRUE,
remove_punct = TRUE, remove_numbers=TRUE)
news_dfm <- dfm_remove(news_dfm, c('*-time', '*-timeUpdated', 'GMT', 'BST'))
news_dfm <- dfm_trim(news_dfm, min_termfreq = 0.95, termfreq_type = "quantile",
max_docfreq = 0.1, docfreq_type = "prop")
news_dfm[1:2, 1:5]
# 20 top words
topfeatures(news_dfm, 20)
# keeping only documents with number of tokens >0 [given that after the above reduction, some articles could actually be empty - 8 here]
news_dfm[ntoken(news_dfm) == 0,]
news_dfm <- news_dfm[ntoken(news_dfm) > 0,]
# as docvars there is the text as well (text2) - it will be helpful later on!
str(news_dfm)
# quanteda does not implement own topic models, but you can easily access to LDA() from the topicmodel package
# through convert(). k = 5 specifies the number of topics to be discovered.
dtm <- convert(news_dfm, to = "topicmodels")
# let's identify k=5
# being a probabilist model, always a good idea to declare the seed!
set.seed(123)
system.time(lda <- LDA(dtm, method= "Gibbs", k = 5))
# alternative way to define the seed:
# system.time(lda <- LDA(dtm, method= "Gibbs", k = 5, control = list(seed = 123)))
# We can use get_terms to the top n terms from the topic model (i.e., terms with the highest beta associated with each topic)
# and get_topics to predict the top k topic for each document. This will help us interpret the results of the model.
# Which the most likely terms for each topic?
terms <- get_terms(lda, 10)
# 10 terms for topic 1
terms[,1]
# 10 terms for topic 2, and so on
terms[,2]
terms[,3]
terms[,4]
terms[,5]
# or we can extract directly the most important terms from the model using terms().
terms(lda, 15)
# which meaning for each topic according to you??
# which is the most likely topic for each document?
topics <- get_topics(lda, 1)
head(topics, 10)
# or you can obtain the most likely topics using topics() and save them as a document-level variable.
head(topics(lda))
docvars(news_dfm, 'pred_topic') <- topics(lda)
str(news_dfm)
# Let’s take a closer look at some of these topics. To help us interpret the output, we can look
# at the words associated with each topic and take a random sample of documents highly associated
# with each topic.
# Topic 1
str(news_dfm)
paste(terms[,1], collapse=", ")
sample(news_dfm@docvars$text2[topics==1], 2)
# Topic 3
paste(terms[,3], collapse=", ")
sample(news_dfm@docvars$text2[topics==3], 2)
# Topic 3
paste(terms[,5], collapse=", ")
sample(news_dfm@docvars$text2[topics==5], 3)
# Looking at the evolution of certain topics over time can be interesting
# Let’s look for example at Topic 5, which appears to be related to Isis
str(news_dfm@docvars$date)
news_dfm@docvars$month <- substr(news_dfm@docvars$date, 1, 7) # extract year and month (first 7 characters)
table(news_dfm@docvars$month )
# frequency table with articles about Isis as their MAIN topic, per month in 2016
tab <- table(news_dfm@docvars$month [news_dfm@docvars$pred_topic==5])
plot(tab)
# But we can actually do better than this. LDA is a probabilistic model,
# which means that for each document, it 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 gamma in the LDA object (theta in the notation we used for the slides)
# For example, news 1 is 22% about topic 1, while news 3 is 38% about topic 1
round(lda@gamma[1,], 2)
round(lda@gamma[3,], 2)
round(lda@gamma[7,], 2)
# So we can actually take the information in the matrix and aggregate it to compute the average probability that a news each month
# is about a particular topic. Let’s now choose once again Topic 5
# add probability to dfm
news_dfm@docvars$topic5 <- lda@gamma[,5]
str(news_dfm)
# now aggregate at the month level
agg <- aggregate(news_dfm@docvars$topic5, by=list(month=news_dfm@docvars$month), FUN=mean)
str(agg)
agg$date <- as.Date(paste(agg$month,"-01",sep=""))
str(agg)
# and plot it
plot(agg$date , agg$x, type="l", xlab="Month", ylab="Avg. prob. of article about topic 5",
main="Estimated proportion of articles discussing about ISIS")
##################################
# identifying the optimal number of topics between for example: 4 and 9
##################################
# other alternatives you could consider are focusing on the perplexity metric or on the log-likelihood of your topic model given different values of K
# let's focus on topic coherence and exclusivity of our models with k=5
topic_diagnostics(lda, dtm)
topic_coherence(lda, dtm)
topic_exclusivity(lda)
mean(topic_coherence(lda, dtm))
mean(topic_exclusivity(lda))
# now let's make k changes between 4 and 9
top <- c(4:15)
top
# let's create an empty data frame that we will fill later on
results <- data.frame(first=vector(), second=vector(), third=vector())
results
# with add iter=100 where iter is the number of Gibbs iterations to make LDA faster; the default value is 2000 [so this is not an advisable strategy
# for a serious project!]
system.time(
for (i in top)
{
set.seed(123)
lda <- LDA(dtm, method= "Gibbs", k = (i), control=list(verbose=50L, iter=100))
topic <- (i)
coherence <- mean(topic_coherence(lda, dtm))
exclusivity <- mean(topic_exclusivity(lda))
results <- rbind(results , cbind(topic, coherence, exclusivity ))
}
)
results
str(results)
# k=9 seems to be the best choice
plot(results$coherence, results$exclusivity, main="Scatterplot Example",
xlab="Semantic Coherence", ylab="Exclusivity ", pch=19)
text(results$coherence, results$exclusivity, labels=results$topic, cex= 1, pos=4)
#########################################
#########################################
################################################################
# Let's estimate a STM
################################################################
#########################################
#########################################
setwd("C:/Users/luigi/Dropbox/TOPIC MODEL")
getwd()
# 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)
# 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))
myDfm <- dfm(myCorpus, remove = stopwords("english"), tolower = TRUE, stem = TRUE,
remove_punct = TRUE, remove_numbers=TRUE)
# Some form of trim is usually suggested when you have many documents and you want to run a STM.
# For example, we can decide to keep only words occurring in at least 2 documents
myDfm.trim <-dfm_trim(myDfm, min_docfreq = 2, verbose=TRUE)
length(myDfm@Dimnames$features) # 6553 features
length(myDfm.trim@Dimnames$features) # 3332 features
# In our case, given the relative low number of documents and features, we are going to proceed with our
# original dfm (i.e., myDfm)
# 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
DfmStm <- convert(myDfm, to = "stm", docvars = docvars(myCorpus))
head(docvars(myCorpus))
str(DfmStm)
################################################################
# 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 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. REMEMBER however: the results you get when you add variables for topic prevalence but no
# variables fro topical content will be (slighly or more) different to the ones when you add variables for both
# topic prevalence and topical content together
# 3) it is 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 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 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 remember always to define a set seed to replicate your analysis in this latter case (as we did with Topic Models)!
# 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"
# 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 28 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. These words provide more semantically intuitive representations of each 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 5 is the most frequent one
# topic meaning: according to frex words, topic 5 seems to be related to the trend in the economy; topic 10 to inflation;
# topic 9 about politics; etc.
plot(stmFitted, type = "hist", labeltype = c("frex")) # Here topic 5 appears as more "evenly" distributed across documents than
# for example topic 11 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 5 for example and let's identify the three 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)
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")
# Let's focus on topic 9 and let's do the same
thought9 <- findThoughts(stmFitted, texts=docs2, topics=9, n=3)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought9, width = 30, main = "Topic 9")
# Let's plot Topic 5 and 9 together
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought5, width = 30, main = "Topic 5 - Economy")
plotQuote(thought9, width = 30, main = "Topic 9 - Politics")
# which are the documents with the highest theta for each k?
apply(stmFitted$theta,2,which.max)
# let's read the entire documents with the highest theta for topic=5
strwrap(texts(myCorpus)[485])
# same one as below
docs2 <- substr(docs, start = 1, stop = 200)
thought5 <- findThoughts(stmFitted, texts=docs2, topics=5, n=1)$docs[[1]]
par(mfrow = c(1, 2),mar = c(.5, .5, 1, .5))
plotQuote(thought5, width = 30, main = "Topic 5")
# 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=5
strwrap(texts(myCorpus)[93])
# which is the document with the 3rd largest value for each topic?
apply(stmFitted$theta, 2, maxn(3))
# let's read the entire documents with the third highest theta for topic=5
strwrap(texts(myCorpus)[67])
# which are the most likely topics across our 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 5 - 5 random documents associated to it
set.seed(123)
sample(nyt2$text[nyt2$topic==5], 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 43% about topic 7 and 28% about topic 4 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 and 15 are positively correlated
plot(stmFitted, type = "summary", labeltype = c("frex"), n=5) # i.e. when an article discusses about Reagan there is a positive
# probability of also discussing about public budget
mod.out.corr$cor
#########################################
# Estimating topic prevalence using as IV: repub and year
#########################################
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 less about topic 10 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 (inflation) is stable but at the end it declines becoming
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 - Inflation")
abline(h=0, col="blue")
# the opposite happens with topics 9 (politics)
plot(prep, "year", method = "continuous", topics = 9,
model = stmFitted, printlegend = FALSE, xaxt = "n", xlab = "Year")
seq <- seq(from = as.numeric("1980"), to = as.numeric("2014"))
axis(1, at = seq)
title("Topic 9 - Politics")
abline(h=0, col="blue")
# Let's now move to repub (i.e., a categorical variable) for topic 10 (inflation) and topic 14 (stock-market)
plot(stmFitted, type = "summary", labeltype = c("frex"), n=5)
# the salience to topic 14 seems to be the same under different partisan presidency. While there is a tendency to discuss more
# about topic 10 under a democratic presidency
plot(prep, "repub", model=stmFitted, method="pointestimate", topic=c(10, 14))
# 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(10, 14),
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 10 - Inflation', 'Topic 14 - Stock-Market'))
#########################################
# 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 222 seconds [4 mins]
# load("C:\\Users\\mw\\Dropbox (VOICES)\\TOPIC MODEL\\Lecture 5\\Material 2020\\STM.RData")
load("C:\\Users\\luigi\\Dropbox\\TOPIC MODEL\\Lecture 5\\Material 2020\\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 10 (inflation) under a democratic or repubblican presidency
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 two topics, normalized by the maximum
# difference in probability of any word between the two topics.
plot(stmFitted_topical, type = "perspectives", topics = c(10, 14))
#########################################
# 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
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
# important: whenever the documents you want to analyze are pretty long, whenever you employ the searchK
# command, it could be advisable also adding this line: N = floor(0.2 * length(Dfmstm2$documents)) -
# where Dfmstm2$documents are the documents that you want to analyze. See below for an example
# storage2 <- searchK(DfmStm $documents, DfmStm $vocab, K, data = DfmStm $meta, max.em.its = 500,
# N = floor(0.2 * length(DfmStm $documents)), prevalence = ~ repub + s(year), init.type = "Spectral")
# str(storage2)
#####################################################
# The semantic coherence-exclusivity ‘frontier’
#####################################################
# You should focus 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 4 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.
# 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.
# 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 see two of them:
library(LDAvis)
library(servr)
toLDAvis(mod=stmFitted, docs=DfmStm $documents, reorder.topics = FALSE)
# Distance between the topics is an approximation of semantic relationship between the topics
# The topic which shares common words will be overlapping (closer in distance) in comparison
# to the non-overlapping topic
# On selecting a topic (clicking on a topic bubble), top 30 words (with the red-shaded area) are shown
# Finally, decreasing the lambda parameter increases the weight of the ratio of the frequency of word given a specific topic.
# Moving the slider allows to adjust the rank of terms based on much discriminatory (or "relevant") are for the
# specific topic: try to move lambda from 1 to 0 and check the results!
# Second: stmBrowser
library(stmBrowser)
str(DfmStm)
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 probelm 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 10", "Y-axis=datetime", "Radius=none", "Color=repub"