### Regularized regression

Our running example will be a random sample of nearly 5,000 tweets mentioning the names of the candidates to the 2014 EP elections in the UK. We will be using a variable named `communication`, which indicates whether each tweet was hand-coded as being __engaging__ (a tweet that tries to engage with the audience of the account) or __broadcasting__ (just sending a message, without trying to elicit a response).

The source of the dataset is an article co-authored with Yannis Theocharis, Zoltan Fazekas, and Sebastian Popa, published in the Journal of Communication. The link is [here](http://onlinelibrary.wiley.com/doi/10.1111/jcom.12259/abstract). Our goal was to understand to what extent candidates are not engaging voters on Twitter because they're exposed to mostly impolite messages.

Let's start by reading the dataset and creating a dummy variable indicating whether each tweet is engaging

```{r}
library(quanteda)
tweets <- read.csv("../data/UK-tweets.csv", stringsAsFactors=F)
tweets$engaging <- ifelse(tweets$communication=="engaging", 1, 0)
tweets <- tweets[!is.na(tweets$engaging),]
```

We'll do some cleaning as well -- substituting handles with @. Why? We want to provent overfitting.
```{r}
tweets$text <- gsub('@[0-9_A-Za-z]+', '@', tweets$text)
```


As we discussed last week, before we can do any type of automated text analysis, we will need to go through several "preprocessing" steps before it can be passed to a statistical model. We'll use the quanteda package quanteda here.

The basic unit of work for the quanteda package is called a corpus, which represents a collection of text documents with some associated metadata. Documents are the subunits of a corpus. You can use summary to get some information about your corpus.

```{r}
twcorpus <- corpus(tweets$text)
summary(twcorpus)
```

We can then convert a corpus into a document-feature matrix using the dfm function. We will then trim it in order to keep only tokens that appear in 2 or more tweets. Note that we keep punctuation -- it turns out it can be quite informative.


```{r}
twdfm <- dfm(twcorpus, remove=stopwords("english"), remove_url=TRUE, 
             ngrams=1:2, verbose=TRUE)
twdfm <- dfm_trim(twdfm, min_docfreq = 2, verbose=TRUE)
```

Note that other preprocessing options are:

- remove_numbers
- remove_punct
- remove_twitter
- remove_symbols
- remove_separators

You can read more in the `dfm` and `tokens` help pages

Once we have the DFM, we split it into training and test set. We'll go with 80% training and 20% set. Note the use of a random seed to make sure our results are replicable.
```{r}
set.seed(123)
training <- sample(1:nrow(tweets), floor(.80 * nrow(tweets)))
test <- (1:nrow(tweets))[1:nrow(tweets) %in% training == FALSE]
```

Our first step is to train the classifier using cross-validation. There are many packages in R to run machine learning models. For regularized regression, glmnet is in my opinion the best. It's much faster than caret or mlr (in my experience at least), and it has cross-validation already built-in, so we don't need to code it from scratch. We'll start with a ridge regression:

```{r}
library(glmnet)
require(doMC)
registerDoMC(cores=3)
ridge <- cv.glmnet(twdfm[training,], tweets$engaging[training], 
	family="binomial", alpha=0, nfolds=5, parallel=TRUE, intercept=TRUE,
	type.measure="class")
plot(ridge)
```

We can now compute the performance metrics on the test set.
```{r}
## function to compute accuracy
accuracy <- function(ypred, y){
	tab <- table(ypred, y)
	return(sum(diag(tab))/sum(tab))
}
# function to compute precision
precision <- function(ypred, y){
	tab <- table(ypred, y)
	return((tab[2,2])/(tab[2,1]+tab[2,2]))
}
# function to compute recall
recall <- function(ypred, y){
	tab <- table(ypred, y)
	return(tab[2,2]/(tab[1,2]+tab[2,2]))
}
# computing predicted values
preds <- predict(ridge, twdfm[test,], type="class")
# confusion matrix
table(preds, tweets$engaging[test])
# performance metrics
accuracy(preds, tweets$engaging[test])
precision(preds, tweets$engaging[test])
recall(preds, tweets$engaging[test])
```

Something that is often very useful is to look at the actual estimated coefficients and see which of these have the highest or lowest values:

```{r}
# from the different values of lambda, let's pick the highest one that is
# within one standard error of the best one (why? see "one-standard-error"
# rule -- maximizes parsimony)
best.lambda <- which(ridge$lambda==ridge$lambda.1se)
beta <- ridge$glmnet.fit$beta[,best.lambda]
head(beta)

## identifying predictive features
df <- data.frame(coef = as.numeric(beta),
				word = names(beta), stringsAsFactors=F)

df <- df[order(df$coef),]
head(df[,c("coef", "word")], n=30)
paste(df$word[1:30], collapse=", ")
df <- df[order(df$coef, decreasing=TRUE),]
head(df[,c("coef", "word")], n=30)
paste(df$word[1:30], collapse=", ")
```

We can easily modify our code to experiment with Lasso or ElasticNet models:

```{r}
lasso <- cv.glmnet(twdfm[training,], tweets$engaging[training], 
	family="binomial", alpha=1, nfolds=5, parallel=TRUE, intercept=TRUE,
	type.measure="class")

```

```{r}
# computing predicted values
preds <- predict(lasso, twdfm[test,], type="class")
# confusion matrix
table(preds, tweets$engaging[test])
# performance metrics (slightly better!)
accuracy(preds, tweets$engaging[test])
precision(preds, tweets$engaging[test])
recall(preds, tweets$engaging[test])
```

```{r}
best.lambda <- which(lasso$lambda==lasso$lambda.1se)
beta <- lasso$glmnet.fit$beta[,best.lambda]
head(beta)

## identifying predictive features
df <- data.frame(coef = as.numeric(beta),
				word = names(beta), stringsAsFactors=F)

df <- df[order(df$coef),]
head(df[,c("coef", "word")], n=30)

df <- df[order(df$coef, decreasing=TRUE),]
head(df[,c("coef", "word")], n=30)
```

We now see that the coefficients for some features actually became zero.

```{r}
enet <- cv.glmnet(twdfm[training,], tweets$engaging[training], 
	family="binomial", alpha=0.50, nfolds=5, parallel=TRUE, intercept=TRUE,
	type.measure="class")
# NOTE: this will not cross-validate across values of alpha

# computing predicted values
preds <- predict(enet, twdfm[test,], type="class")
# confusion matrix
table(preds, tweets$engaging[test])
# performance metrics
accuracy(preds, tweets$engaging[test])
precision(preds, tweets$engaging[test])
recall(preds, tweets$engaging[test])

best.lambda <- which(enet$lambda==enet$lambda.1se)
beta <- enet$glmnet.fit$beta[,best.lambda]
head(beta)

## identifying predictive features
df <- data.frame(coef = as.numeric(beta),
				word = names(beta), stringsAsFactors=F)

df <- df[order(df$coef),]
head(df[,c("coef", "word")], n=30)

df <- df[order(df$coef, decreasing=TRUE),]
head(df[,c("coef", "word")], n=30)

```


### Xgboost

If we really want the best performance at a low computational cost, the cutting-edge method many people are using is Distributed Gradient Boosting, based on the same ideas as boosted trees / random forests, implemented as `xgboost`. You can read more about the history of this package [here](https://homes.cs.washington.edu/~tqchen/2016/03/10/story-and-lessons-behind-the-evolution-of-xgboost.html).

First, let's prepare the data...

```{r}
library(quanteda)
tweets <- read.csv("../data/UK-tweets.csv", stringsAsFactors=F)
tweets$engaging <- ifelse(tweets$communication=="engaging", 1, 0)
tweets <- tweets[!is.na(tweets$engaging),]
# clean text and create DFM
tweets$text <- gsub('@[0-9_A-Za-z]+', '@', tweets$text)
twcorpus <- corpus(tweets$text)
twdfm <- dfm(twcorpus, remove=stopwords("english"), remove_url=TRUE, 
             ngrams=1:2, verbose=TRUE)
twdfm <- dfm_trim(twdfm, min_docfreq = 2, verbose=TRUE)
# training and test sets
set.seed(123)
training <- sample(1:nrow(tweets), floor(.80 * nrow(tweets)))
test <- (1:nrow(tweets))[1:nrow(tweets) %in% training == FALSE]
```

Now we can train the model:

```{r}
library(xgboost)
# converting matrix object
X <- as(twdfm, "dgCMatrix")
# parameters to explore
tryEta <- c(1,2)
tryDepths <- c(1,2,4)
# placeholders for now
bestEta=NA
bestDepth=NA
bestAcc=0

for(eta in tryEta){
  for(dp in tryDepths){	
    bst <- xgb.cv(data = X[training,], 
		    label =  tweets$engaging[training], 
		    max.depth = dp,
	      eta = eta, 
	      nthread = 4,
	      nround = 500,
	      nfold=5,
	      print_every_n = 100L,
	      objective = "binary:logistic")
    # cross-validated accuracy
    acc <- 1-mean(tail(bst$evaluation_log$test_error_mean))
		cat("Results for eta=",eta," and depth=", dp, " : ",
				acc," accuracy.\n",sep="")
		if(acc>bestAcc){
				bestEta=eta
				bestAcc=acc
				bestDepth=dp
		}
	}
}

cat("Best model has eta=",bestEta," and depth=", bestDepth, " : ",
	bestAcc," accuracy.\n",sep="")
   
```

How well does it perform out-of-sample?

```{r}
# running best model
rf <- xgboost(data = X[training,], 
    label = tweets$engaging[training], 
		max.depth = bestDepth,
    eta = bestEta, 
    nthread = 4,
    nround = 1000,
		print_every_n=100L,
    objective = "binary:logistic")

# out-of-sample accuracy
preds <- predict(rf, X[test,])

## function to compute accuracy
accuracy <- function(ypred, y){
	tab <- table(ypred, y)
	return(sum(diag(tab))/sum(tab))
}
# function to compute precision
precision <- function(ypred, y){
	tab <- table(ypred, y)
	return((tab[2,2])/(tab[2,1]+tab[2,2]))
}
# function to compute recall
recall <- function(ypred, y){
	tab <- table(ypred, y)
	return(tab[2,2]/(tab[1,2]+tab[2,2]))
}

cat("\nAccuracy on test set=", round(accuracy(preds>.50, tweets$engaging[test]),3))
cat("\nPrecision on test set=", round(precision(preds>.50, tweets$engaging[test]),3))
cat("\nRecall on test set=", round(recall(preds>.50, tweets$engaging[test]),3))
```

What we sacrifice is interpretability (yet again!). We can check feature importance, but it's often hard to tell what's going on exactly. Why? We only what features "matter", but not why!


```{r}
# feature importance
labels <- dimnames(X)[[2]]
importance <- xgb.importance(labels, model = rf, data=X, label=tweets$engaging)
importance <- importance[order(importance$Gain, decreasing=TRUE),]
head(importance, n=20)

# adding sign
sums <- list()
for (v in 0:1){
	sums[[v+1]] <- colSums(X[tweets[,"engaging"]==v,])
}
sums <- do.call(cbind, sums)
sign <- apply(sums, 1, which.max)
	
df <- data.frame(
	Feature = labels, 
	sign = sign-1,
	stringsAsFactors=F)
importance <- merge(importance, df, by="Feature")
	
## best predictors
for (v in 0:1){
	cat("\n\n")
	cat("value==", v)
	importance <- importance[order(importance$Gain, decreasing=TRUE),]
	print(head(importance[importance$sign==v,], n=50))
	cat("\n")
	cat(paste(unique(head(importance$Feature[importance$sign==v], n=50)), collapse=", "))
}

```



