```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(quanteda)
library(quanteda.textstats)
library(quanteda.textplots)
library(quanteda.textmodels)
library(caret)
library(stm)
library(LSX) 
```

# Bags-of-words

Let's load our input data. For the tutorial today and tomorrow, we are using data from the Manifesto Project that codes quasi-sentences of party programmes into different categories. Its main purpose is to trace the emphasis parties put on different issues during their campaigns. However, as it has directional parameters as well (e.g. *Protectionism: Positive* and *Protectionism: Negative*), it is also often used to infer party positions. 

In an additional script, I downloaded the quasi-sentences of all English-speaking countries, and transformed it into a data frame. Let's load it into our environment.
```{r}
df <- readRDS("./data_prep.RDS")

# let us remove texts without an issue
df <- subset(df, issue!="NA")

# create a new id variable that has a unique identifier for each row
df$id <- 1:nrow(df)
summary(df)
table(df$country)
table(df$parfam_n)
```

This data frame consists of different forms of sentences and has some meta information (country names, party names, party family, issue of a sentence, year).

Let's look into the first five texts
```{r}
df$sentence[1:5]
```


# Data preparation
Let's first transform our data frame into a corpus. 
```{r}
corp <- corpus(df$sentence)
head(corp)
summary(head(corp))
```

In this example, each element represents one text. 

The text are still not separated at all. For bags-of-words, we need a tokenizer that throws all tokens of a text into one bag. The tokenizer we use here separates each word, for other applications, we may need a sentence tokenizer that splits texts into sentences. 
```{r}
toks <- tokens(corp, what="word")

# let's look into the first 20 tokens of our first document
head(toks[[1]], 20)
```

Now, that's good but it still contains a bit of garbage. 

Generally, in all text-as-data tasks, we want to ensure that we work with the smallest object possible for our purpose (without losing any information). For transformer models, that does not really work as we need all of the context available. But for bags-of-words, we want to get rid of all the unnecessary tokens (like white space, punctuation, links). We can do this in two different ways: either by using regular expressions (https://cran.r-project.org/web/packages/stringr/vignettes/regular-expressions.html) or by using quanteda's built-in structures.

For simplicity, we rely on quanteda's functions first but generally, a profound knowledge of regular expressions is very helpful!

```{r}
toks <- tokens(corp, what="word",
               remove_punct = T,
               remove_numbers = T,
               remove_symbols = T,
               remove_url = T)
 
# again, 20 tokens in first document
head(toks[[1]], 20)
```
That looks better already. But there are still some terms which have very little meaning in substance. We call these terms stopwords (like "the", "and", "in", etc.). Luckily, quanteda has a built-in function to remove these as well!

*Be careful! Sometimes, you want to keep stopwords, and as this function depends on a list of words, some of them may be missing*

```{r}
# use quanteda's list of stopwords
head(quanteda::stopwords(language="en"))

# if you need to add stopwords
stopwords2 <- c(quanteda::stopwords(language="en"), "can")

# remove stopwords
toks <- toks %>% tokens_remove(stopwords2)
 
# again, 20 tokens in first document
head(toks[[1]], 20)
```

Now, we are ready to transform our character vectors into numeric representations. For most bags-of-words models, we want to transform it into a *d*ocument-*f*requency-*m*atrix (*dfm*). A dfm is a matrix whose rows correspond to the document and its columns to tokens. The cells denote how often a token appears in a text. 

Let's do that, and transform characters to lower case at the same time.
```{r}
m_dfm <- dfm(toks, tolower = T)

m_dfm
```

This is how a *dfm* looks like. The meta information is helpful here: a dfm is always quite sparse, that means that most of the cells are empty (because many texts do no contain the token). The more we can reduce the sparseness, the better. One common pre-processing procedure is therefore to remove tokens that are extremely rare. At the same time, those tokens which are extremely common may also be not that informative (part of many documents). So let's remove those as well.

Let's drop all of the terms that occur less than five times, and those tokens that are used in more than 50 percent of the documents.
```{r}
m_dfm <- dfm(toks, tolower = T) %>% dfm_trim(min_termfreq = 5,
                                             termfreq_type = "count",
                                             max_docfreq = 0.5,
                                             docfreq_type = "prop")
m_dfm
```

The dfm is still very sparse but we reduced the number of features from 78k to 15k, that's a huge reduction in size. 

What's still missing are meta information. We did not define any docvars (like issue, party, etc.). Let's do that now.

```{r}
docvars(m_dfm, "text_id") <- df$id
docvars(m_dfm, "country") <- df$country
docvars(m_dfm, "issue") <- df$issue
docvars(m_dfm, "parfam") <- df$parfam_n
docvars(m_dfm, "year") <- df$year

m_dfm


# let's do the same for the corpus
docvars(corp, "text_id") <- df$id
docvars(corp, "country") <- df$country
docvars(corp, "issue") <- df$issue
docvars(corp, "parfam") <- df$parfam_n
docvars(corp, "year") <- df$year
```

As we removed so many features, it is possible that some documents are empty now. This is not a problem but when doing analysis later, this may affect our evaluation. So, let's only keep those documents containing more than 0 tokens. Drop empty texts also from data frame. This facilitates later merging of the data with our data frame.

```{r}
# create vector that contains ids of texts that are not empty
drop_ids <- m_dfm@docvars$docname_[ntoken(m_dfm)==0] %>% str_extract(., "\\d{1,}") %>% as.numeric()

# drop from data frame
df <- df[!row.names(df) %in% drop_ids, ]

# drop from dfm
m_dfm <- dfm_subset(m_dfm, ntoken(m_dfm)>0)
m_dfm

# drop from corp
corp <- corpus_subset(corp, !text_id %in% drop_ids)
```

We see that almost 1,000 documents have actually been removed because of that. 


Now that our dfm is prepared, let's inspect it. 
What are the most frequent terms?
```{r}
topfeatures(m_dfm)
```

...and which terms are used in most documents?
```{r}
topfeatures(m_dfm, scheme="docfreq")
```


# Descriptive visualizations

It's often quite fruitful to visualize some of the relations between words without doing any complex analysis. 

In the past, many people did wordclouds. Although they are not that informative, it may look good on a cover page or in a report.
```{r}
textplot_wordcloud(m_dfm, max_words = 300, color="firebrick")
```


## Readability
We can also visualize the complexity of the texts which is often denoted by the "Flesch" score. The higher the score, the easier a text is to read. 
```{r}
complexity <- textstat_readability(corp)

# add some meta information
complexity <- cbind(complexity, df$partyname, df$year)

head(complexity, 20)
```

Generally, there is quite a variation in the readability of sentences (a score between 30-50 is usually considered as college-level difficulty).

## Keywords in context
Now, let's get a glimpse of the data by looking at how different words are used in context (something we will come back to tomorrow when talking about word embeddings). 

The kwic method is unfortunately only available for tokens objects, so we have to use the tokens object that includes a few more garbage.

```{r}
kwic(toks, "country",window=4)
```
This gives us an idea of how different actors speak about their country. 


## Dictionary 
We can pre-define a dictionary with dictionary(). Each list object denotes a category which is comprised of a character vector of words associated. We can use quantifiers for word stems (refugee* includes, among others, "refuge", "refugee", "refugees"). 

After using dfm_lookup, we can group the dfm by a docvar (here party family) to find out who is talking how often about a specific category.

```{r}
# define dictionary
dict <- dictionary(list(economy = c("deficit", "budget", "growth", "employment"),
                        migration = c("immigrants", "refuge*", "asylum")))

# apply dictionary to dfm
(dict_results <- dfm_lookup(m_dfm, dict))

# retrieve category prevalence by group
dfmat_parfam_dict <- dfm(dict_results) %>% 
  dfm_group(groups = parfam)
dfmat_parfam_dict
```


## Keyness
Keyness statistics help us to know more about how different terms are used differently across groups. 
For this purpose, the meta-information we stored as docvars is useful. 
In this example, we want to know what distinguishes mainstream left parties from all other parties. 

```{r}
keyness <- dfm_group(m_dfm, parfam) %>%
 textstat_keyness("Mainstream Left")

head(keyness)

keyness %>%
 textplot_keyness()
```
We see some of the typical terms associated with centre-left parties, like workers and multinationals. But there are still many terms that are just self-references (labour) or other party references (tories).


# Unsupervised approaches

The first step is done, so let's do a proper analysis of our data. 

We start with a simple scaling task. Some of the first applications in polsci of CTA have done exactly that (e.g. Laver et al. 2003).

Wordfish is usually done on larger documents, not on single sentences. So, for this purpose, we should change the structure of our data. That's a good time to practice what you've learned before. 

Furthermore, for our showcase, let us focus on the UK only.

```{r}
df_text <- df %>%
  filter(country=="United Kingdom" & partyname %in% c("Labour Party", "Conservative Party")) %>%
  group_by(partyname, parfam_n, year) %>%
  summarise(text = toString(sentence)) 
```


***Your Turn***
Pre-process the new data frame.
```{r}

```






...and now do it together to have the same data frame
```{r}
corp_text <- corpus(df_text$text)
toks_text <- tokens(corp_text, what="word",
               remove_punct = T,
               remove_numbers = T,
               remove_symbols = T,
               remove_url = T)


m_dfm2 <- dfm(toks_text, tolower = T) # for the moment, let us not be too strict, and remove any more features. After all, our data frame is small already.


docvars(m_dfm2, "parfam") <- df_text$parfam_n
docvars(m_dfm2, "year") <- df_text$year
docvars(m_dfm2, "partyname") <- df_text$partyname
docnames(m_dfm2) <- df_text$year

m_dfm2 <- dfm_subset(m_dfm2, ntoken(m_dfm2)>0)
```



First, let's start with wordfish.

```{r}
m_wordfish <- textmodel_wordfish(m_dfm2) 
summary(m_wordfish)
textplot_scale1d(m_wordfish, groups=df_text$partyname)
```
Well, we've got the results. What do they depict? Generally, we see the theta score for each party for each year it has published a manifesto. The theta increases with time. But within the years, there are some differences between parties as well. 

For instance, comparing the theta score for the Conservative and Labour Party in 1964. 
```{r}
df_wordfish <- cbind(theta=m_wordfish$theta, m_wordfish$x@docvars) 

print(paste0(df_wordfish$partyname[df_wordfish$year=="1964"], ": ", df_wordfish$theta[df_wordfish$year=="1964"]))
```
The Labour Party has a smaller score. 

```{r}
df_wordfish %>%
  group_by(partyname) %>%
  summarise(position = mean(theta, na.rm=T))
```
On average, the scores are very similar, though. We see that they are symmetric.

Let's look into the scores each term has been assigned to.
```{r}
# create a data frame with the scores for each word
words_wordfish <- data.frame(term = m_wordfish$features, score = m_wordfish$beta)

words_wordfish %>%
  arrange(score) %>%
  head(20)

words_wordfish %>%
  arrange(score) %>%
  tail(20)
```
We see that words we might associate with left or right are not clearly scaled as such. Maybe, we simply have too many different categories in our data? After all, parties do talk about many issues. Let's leverage on the fact that we have information about the issues parties talk about and create a new data frame on the basis of economic issues.


```{r}
df_text <- df %>%
  filter(countryname=="United Kingdom" & partyname %in% c("Labour Party", "Conservative Party") & issue %in% as.character(c(401:417))) %>%
  group_by(parfam_n, partyname, year) %>%
  summarise(text = toString(sentence)) 


corp_text <- corpus(df_text$text)
toks_text <- tokens(corp_text, what="word",
               remove_punct = T,
               remove_numbers = T,
               remove_symbols = T,
               remove_url = T)


m_dfm2 <- dfm(toks_text, tolower = T) # for the moment, let us not be too strict, and remove any more features. After all, our data frame is small already.


docvars(m_dfm2, "parfam") <- df_text$parfam_n
docvars(m_dfm2, "partyname") <- df_text$partyname
docvars(m_dfm2, "year") <- df_text$year
docnames(m_dfm2) <- df_text$year

m_dfm2 <- dfm_subset(m_dfm2, ntoken(m_dfm2)>0)
```
A new attempt:
```{r}
m_wordfish <- textmodel_wordfish(m_dfm2) 
summary(m_wordfish)
textplot_scale1d(m_wordfish, groups=df_text$partyname)
```
This looks a bit more specific, there are notable difference between Labour and Conservatives. There's still a general time trend towards higher theta scores but within the years. Labour has always been assigned a bigger score.

If we now look into the terms associated with each side:

```{r}
# create a data frame with the scores for each word
words_wordfish <- data.frame(term = m_wordfish$features, score = m_wordfish$beta)

words_wordfish %>%
  arrange(score) %>%
  head(20)

words_wordfish %>%
  arrange(score) %>%
  tail(20)
```
These terms are still a bit messy, but some of them tell us a bit more about ideology (like tax-free, recovery, supplies, recession on the right; and greed and social groups like bame on the left).


# Topic models
Unsupervised scaling works mediocre, after all it's a very simplistic approach. We will later learn better approaches to predict the positions of a text. Bear with me, let's move to something else. A very common approach of CTA is to get a better idea which topics a text covers. 

Let's do a topic model on a small subset of text first (5% of the text). This still takes 3 minutes.
```{r}
start <- Sys.time()
set.seed(3145)
m_dfm_sub <- dfm_sample(m_dfm, size=0.05*ndoc(corp))

# drop empty features 
m_dfm_sub <- dfm_trim(m_dfm_sub, min_termfreq = 1)

# estimate stm
stm_sub <- stm(m_dfm_sub, K = 20, seed=421, emtol=0.001)
end <- Sys.time()
end-start
```


## Interpretation
I've estimated a topic model with the full amount of data, if you want load it and compare to your own findings
```{r}
stm1 <- readRDS("./stm_marpor.RDS")
labelTopics(stm_sub)
labelTopics(stm1)
```

Examples of these topics in the text would help, right?
```{r}
findThoughts(stm1,
             df$sentence,topics=11,n=15)
```

We see some texts several times (that's an issue with Australia in which parties contested as electoral alliances on same programmes). But apart from that, it gives us a good idea on what the topic might be about.

We can also use some in-built functions of stm to plot words typical terms coming with certain topics. Let us plot topic 8 (which seems to be about the economy) and topic 11 (about health).

```{r}
plot(stm1, 
     type="perspectives", 
     topics=c(8,11), 
     plabels = c("Economy","Health"))
```


Many of these topics are quite intuitive, so let's give them a name and add the probability each document has for topic k to our data frame. 
```{r}
df_topic <- bind_cols(df, stm1$theta)

# now assign plausible names to columns
colnames(df_topic)[15:34] <- tolower(c("climate", "agriculture", "banking", "energy", "development", "infrastructure", "minority", "innovation", "employment", "ext_relations", "health", "economy", "party", "housing", "sea", "education","crime" ,"rights" ,  "trash", "trash2"))
```


Every text has 20 topics assigned, having a probability for each of these. Usually, we want to classify texts into *one* topic. Let us do that by assigning the topic with the highest probability to each text. This is computationally intense. We do not really have time to go into different ways to speed up your code (like parallelization). But generally, it's always good to look for a *data.table* solution. Data table is a similar structure than data frame but usually much faster (if you got time, you can check the difference to the *dplyr* version below).

```{r}
library(data.table)
# convert to data table
df_topic <- data.table(df_topic)

# return the column name which has the highest number in column 15-34 (our columns with topic probabilities)
df_topic[, topic := names(.SD)[max.col(.SD)], .SDcols = 15:34]

# convert back to data frame
df_topic <- data.frame(df_topic)

# dplyr code: takes ages!
#df_topic <- df_topic %>%
#  rowwise() %>%
#  mutate(topic = names(.)[which.max(c_across(climate:trash2))])
```


This leaves us with an object we can deal with. A proper data frame. 
Let's look how many topics we have, and visualize their frequency. 
```{r}
df_topic %>%
  ggplot(aes(topic))+ 
  geom_bar()  + theme_light()+ theme(axis.text.x=element_text(angle=90))
```

We've got many posts on health and the economy, few on trash - that's good! 

You can now run all sorts of analysis.
```{r}
summary(m1 <- lm(economy ~ parfam_n + as.factor(year), df_topic))
```

For instance, we can see that the mainstream left and right are talking more about the economy than Green parties (and surprisingly the radical left). This is of course just a superficial view into the data. Usually, it makes sense to further aggregate the topics.


## Advancements to the standard topic model

As we often want to know how different actors talk about the same topic, we can use the *content* option by *stm*. Let us re-estimate the topic model but add the party family as a covariate. 


```{r}
stm2 <- stm(m_dfm, K = 20, seed=421, emtol=0.001, content=~parfam)
saveRDS(stm2, "./stm_marpor_parfam.RDS")
labelTopics(stm2)
```

In this case, results are rather poor. The covariate words are a bit informative with the Radical Left using the term "neoliberal" more often while the Greens talk more about tax deduction. But the interactions between topics and terms lead to very specific keywords. 


## Validation

For the technical validation, we could try out different numbers of k to evaluate which lead to the lowest error.
To do this, we can use the function searchK which requires documents and vocabulary as inputs. This is only stored in a stm object. However, we can transform our dfm object into stm. 


(takes more than 6 minutes)
```{r}
start <- Sys.time()
set.seed(145)
# convert dfm to stm
m_stm_sub <- convert(m_dfm_sub,"stm")

# create vector with different numbers for k (here, 5, 10 and 15)
K<- c(5,10,15) 

# run validation
best_k <- searchK(m_stm_sub$documents, m_stm_sub$vocab, K)
end <- Sys.time()
end-start
saveRDS(best_k, "./searchK.RDS")
```

Let's plot to see which number of topics is best.
```{r}
plot(best_k)
```
Based on the metrics, we should select K with high likehlihood (the closer to zero, the better), high semantic coherence, high lower bound and low residuals. It's (always) a tough choice, ultimately rather depending on theoretical considerations. Semantic coherence is usually high when number of topics is low, residuals rather decrease with more topics. 


# Semi-supervised scaling
In the chunks above, we have already tried out wordfish to scale the positions of political parties. These results were very poor. Wordfish is agnostic about any distribution of terms. As social scientists, we might have certain expectations of how parties from the Left and Right communicate. We can use these expectations as priors and feed them into a model that's called Latent Semantic Scaling, developed by Watanabe. 

We start by building a dictionary on left- and right-wing terms (let's focus on the economy here). 
```{r}
lr_dict <- dictionary(list(left = c("unemployment", "justice", "wage", "employee", "bargaining"),
                           right = c("budget", "merit", "deficit", "business", "growth")))
seed <- as.seedwords(lr_dict)
```

Let us construct the dfm of interest again, here filtering to UK and topics on the economy. We do not need to aggregate to a party level here, we want a sentence level.  
```{r}
df_text <- df %>%
  filter(issue %in% as.character(c(401:417))) 


corp_text <- corpus(df_text$sentence_context)
toks_text <- tokens(corp_text, what="word",
               remove_punct = T,
               remove_numbers = T,
               remove_symbols = T,
               remove_url = T)


m_dfm2 <- dfm(toks_text, tolower = T) # for the moment, let us not be too strict, and remove any more features. After all, our data frame is small already.


docvars(m_dfm2, "parfam") <- df_text$parfam_n
docvars(m_dfm2, "partyname") <- df_text$partyname
docvars(m_dfm2, "year") <- df_text$year
docvars(m_dfm2, "text_id") <- df_text$id
docvars(m_dfm2, "country") <- df_text$country

m_dfm2 <- dfm_subset(m_dfm2, ntoken(m_dfm2)>0)
```

Now, we are already ready to feed them into our model. 
```{r}  
lss_model <- textmodel_lss(m_dfm2, seeds = seed,
                         k = 300, auto_weight = T)
```

Let's visualize the words surrounding our polarising words 
```{r}
textplot_terms(lss_model, lr_dict, max_words = 1000)
```

The LSS package has generated a weight for every term (whether it's rather left or right wing). Generally, the closer terms are to the polarizing words, the larger their weight. We can use these weights to predict a general score along the two-dimensional continuum for each of our sentences.
```{r}
dfmat_doc <- dfm_group(m_dfm2)
dat <- docvars(m_dfm2)
dat$fit <- predict(lss_model, newdata = dfmat_doc, rescale = T)
```

...now we can summarize these scores: for example, by party family.
```{r}
dat %>%
  group_by(parfam) %>%
  summarise(position = mean(fit, na.rm=T))
```
...or plot them
```{r}
dat %>%
  filter(partyname %in% c("Labour Party", "Conservative Party") & country=="United Kingdom") %>%
  group_by(partyname, year) %>%
  summarise(position = mean(fit, na.rm=T)) %>%
  ggplot(aes(year, position, color=partyname)) +
  geom_line() + geom_point() + theme_light()
```

These results look a bit more plausible. The jumps are still a bit dubious. Generally, we could improve our measure by a lot more data. 


# Supervised prediction task 

The unsupervised task has given us more information on the content of our corpus. We know however already quite a lot about this since it has been coded by humans. 

Preparation  of our data frame. Let's create a binary classification task in which we predict whether a text is about the welfare state (per504 or per505) or migration (per601 or per602). After this, we need to subset our data frame to align it with our dfm. 
```{r}
df_sub <- subset(df, issue %in% c("504", "505", "601", "602"))

df_sub$label <- ifelse(df_sub$issue %in% c("504", "505"), 1, 0)
saveRDS(df_sub, "./pred_df.RDS") # save for later

# now that we have a smaller data frame, we need to adjust our dfm accordingly
keep <- df_sub$id
m_dfm_sub <- dfm_subset(m_dfm, text_id %in% keep)

# assign label as docvar to dfm
docvars(m_dfm_sub, "label") <- df_sub$label
```


Split the data into training and test data. The ratio is up to you, often people choose 80 (training)/20 (test) splits. 
```{r}
# make a test training split (this involves taking a random sample, for reproducibility, use a seed!)
set.seed(250)

# optional: further reducing the size of our dfm to increase speed
# m_dfm_small <- dfm_sample(m_dfm_sub,0.3*ndoc(corp))


train <- dfm_sample(m_dfm_sub,0.8*nrow(df_sub))
test <- dfm_subset(m_dfm_sub,
  !(docnames(m_dfm_sub) %in% docnames(train)))
```


Now, we've got two dfms. As we sampled only some texts from these, some features (words) can be dropped from the dfm. Morevoer, it's possible that they have different document frequency matrices. Let's match them.
```{r}
train <- train %>% dfm_trim(1)
test <- dfm_match(test, featnames(train))
```


We are ready to use the classification algorithms to fit our data. 
We follow a Naive Bayes approach first. It follows Bayes Theorem - given a word w in a text, how likely is it that it belongs to category A or B. It's very straightforward and extremely fast. 
```{r}
start <- Sys.time()

# Model
nb_model<-textmodel_nb(train,docvars(train, "label"))
nb_model

end <- Sys.time()
end-start
```

Predicting test
```{r}
test_predictions<-predict(nb_model, newdata=test)

head(test_predictions,5)
(eval_mat <- table(docvars(test,"label"),test_predictions))
```

How did we perform?
```{r}
print(paste0("There are ", eval_mat[[2]] + eval_mat[[3]], " wrong predictions"))
```

With caret, we can also show more sophisticated evaluation measures
```{r}
confusionMatrix(test_predictions, as.factor(docvars(test,"label")))
```

How predictive are certain features of a class?
```{r}
coefs_nb <- coef(nb_model) %>% 
  data.frame() %>%
  mutate(feature=rownames(.))

# words that are more predictive of migration
coefs_nb[coefs_nb$feature=="migration",]
coefs_nb[coefs_nb$feature=="asylum",]

# words that are more predictive of welfare
coefs_nb[coefs_nb$feature=="pension",]
coefs_nb[coefs_nb$feature=="health",]
```

Ok, that's not too bad. To have a better comparison for other supervised tasks we will do next week, let's keep the train test ids.
```{r}
train_ids <- docvars(train)["text_id"] %>% unlist()
test_ids <- docvars(test)["text_id"] %>% unlist

saveRDS(train_ids, "./train_ids.RDS")
saveRDS(test_ids, "./test_ids.RDS")
```


There are other models as well we could try. One of these models is support vector machines. 
```{r}
start <- Sys.time()

# Model
nb_model<-textmodel_svm(train,docvars(train,
  "label"))
nb_model

end <- Sys.time()
end-start



test_predictions<-predict(nb_model,
  newdata=test)

head(test_predictions,5)
(eval_mat2 <- table(docvars(test,"label"),test_predictions))
```
```{r}
print(paste0("There are ", eval_mat[[2]] + eval_mat[[3]], " wrong predictions"))
```

That's already an improvement, and that concludes our session. 