rm(list=ls(all=TRUE)) getwd() # with setwd() you change the current working directory of the R process setwd("C:/Users/mw/Dropbox (VOICES)/TOPIC MODEL") getwd() library(quanteda) library(readtext) library(ggplot2) library(stm) library(igraph) library(wordcloud) ################################################################ # 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. meta.pr <- read.csv("C:/Users/mw/Dropbox (VOICES)/TOPIC MODEL/List President USA party.csv", 1) str(meta.pr) levels(meta.pr$Party) table(meta.pr$Party) # Suggestion: if you want to include in a corpus a Factor variable that you are going to include as a covariate # in a STM, be sure that ALL the levels of such categorical variable will be used in the analysis. This is NOT what # happens in our case, given that by focusing just on speeches made after 1970 in our example below, we will have 0 observations # for the levels Democratic-Republican, Federalist, Unaffiliated, Whig! So let's just keep the levels for Democratic and Repubblican! levels(meta.pr$Party) <- list(Democratic="Dem", Republican="Rep") levels(meta.pr$Party) table(meta.pr$Party) docvars(data_corpus_inaugural, "Party") <- meta.pr$Party head(summary(data_corpus_inaugural)) myCorpus <- corpus_subset(data_corpus_inaugural, Year > 1970) summary(myCorpus) # inspect the document-level variables (that's very important as an option, as we will see..) head(docvars(myCorpus )) myDfm <- dfm(myCorpus, remove = stopwords("english"), tolower = TRUE, stem = TRUE, remove_punct = TRUE, remove_numbers=TRUE) # if you want, you can decide to keep only words occuring >=10 times and in >=2 docs # Some form of trim is HIGHLY suggested when you have many documents and you want to run a STM myDfm.trim <- dfm_trim(myDfm, min_termfreq = 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) # 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! # Note that in this case we keep all the words for simplicity as an example (no trimming!) Dfmstm2 <- convert(myDfm, to = "stm", docvars = docvars(myCorpus)) ################################################################ # 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 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) you can also run a model w/o covariates. In this case you have an unstructured # topic model, by fitting a Correlated Topic Model (see below for an example); # 4) 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 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! str(Dfmstm2) 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 ################################################################ # Summary visualization: understanding the extracted topics through words and example documents # First: which words/topic association? labelTopics(stmFitted) plot(stmFitted, type = "summary", labeltype = c("frex")) plot(stmFitted, type = "labels", labeltype = c("frex")) plot(stmFitted, type = "summary") # words with the highest prob. cloud(stmFitted, topic = 3) # it works with the words with the highest prob. # note these grahps also show you the frequency distribution of topics across the 12 documents plot(stmFitted, type = "hist", labeltype = c("frex")) # Let's read the documents most associated with each topic docs <- myCorpus$documents$texts str(docs) # Let's focus on topic 2 for example... 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") # you cannot read it properly! # 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") # Let's focus on topic 5 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 plot Topic 2 and 5 together 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? First let's extract the thetas theta <-as.data.frame(stmFitted$theta) str(theta ) # Then let's combine the thetas with the documents combine <- paste(myCorpus$documents$President, myCorpus$documents$Year, sep = " ") combine rownames(theta ) <- combine rownames(theta ) # which are the the most likely topics across our documents? apply(theta ,1,which.max) # you can save them as a document-level variable. docvars(myDfm , 'STMtopic') <- apply(theta ,1,which.max) docvars(myDfm) # 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 and 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]]) ################################################################ # 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! mod.out.corr <- topicCorr(stmFitted) plot(mod.out.corr) mod.out.corr$cor library(corrplot) 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 topic prevalence using as IV: Party and Year ######################################### prep <- estimateEffect(1:5 ~ Party+ Year, stmFitted, meta = Dfmstm2$meta, uncertainty = "Global") summary(prep) summary(prep, topics=5) # 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) plot(prep, "Year", method = "continuous", topics = 5, model = stmFitted, printlegend = FALSE, xaxt = "n", xlab = "Year") seq <- seq(from = as.numeric("1970"), to = as.numeric("2017")) axis(1, at = seq) title("Topic 5") abline(h=0, col="blue") # Topic 4 appears as a Democratic one, while Topic 5 as a Repubblican one plot(prep, "Party", model=stmFitted, method="pointestimate") # However only for the Topic 4 case there is a statistical significant difference between the 2 parties 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") prep2 <- 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. # As you can see now, the effect of Year on Topic 5 prevalance is significant but ONLY for Repubblican Presidents, and only up to early 90s plot(prep2, covariate = "Year", model = stmFitted_interaction, method = "continuous", xlab = "Year", ylab="Expected Topic 5 Proportion", 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") ######################################### # Estimating topical content using as IV: Party ######################################### # Let's run a topical content analysis by including as a covariate Party. # 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) 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. # For example, one could estimate a STM model for 3, 4 and 5 topics and compare the results along each of the criteria # But REMEMBER: nothing replaces your interpretation!!!! K <-c(3,4,5) # 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 = 500, prevalence = ~ Party +Year, data = Dfmstm2$meta, init.type = "Spectral") str(storage) # important: whenever the documents you want to analyze are pretty long, whenever you employ the searchK # command, it is always 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! # Do not run in the Lab: # storage2 <- searchK(Dfmstm2$documents, Dfmstm2$vocab, K, data = Dfmstm2$meta, max.em.its = 500, # N = floor(0.2 * length(Dfmstm2$documents)), prevalence = ~ Party +Year, init.type = "Spectral") # str(storage2) # 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: # 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) # In addition to these functions one can also explore if there are 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 extensions about STM (if you like them, take a look also at library(stmCorrViz), library(LDAvis) ##### and library(stminsights) among the others ######################################### # 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", "Party"), text="docs", labeltype = "frex") # select in the HTML file you produced "X-axis=Topic 4", "Y-axis=Year", "Radius=none", "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!] # 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! stmBrowser(stmFitted, data= Dfmstm2$meta, c("Year2", "Year", "Party"), text="docs", labeltype = "frex") ######################################### ########## Topic-model (LDA without topic correlations) ######################################### library(topicmodels) # 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(myDfm , to = "topicmodels") lda <- LDA(dtm, k = 5) # You can extract most important terms from the model using terms(). terms(lda, 5) # Let's compare them with the results we got via our STM model labelTopics(stmFitted, topic=1, n=5) labelTopics(stmFitted, topic=2, n=5) labelTopics(stmFitted, topic=3, n=5) labelTopics(stmFitted, topic=4, n=5) labelTopics(stmFitted, topic=5, n=5) # you can then obtain the most likely topics using topics() and save them as a document-level variable. topics(lda) docvars(myDfm , 'LDAtopic') <- topics(lda) # let's compare the results with the ones we got via STM docvars(myDfm) # You can also estimate a correlated topics model from topicmodels by typing the below command. # But do NOT run it, unless you are ready to wait for a couple of hours... # ctm <- CTM(dtm, k = 5) # The above command would be substantially the same as running stm w/o any covariates # stmFitted <- stm(Dfmstm2$documents, Dfmstm2$vocab, K = 5, max.em.its = 75, init.type = "Spectral")