In this notebook we will work through a basic classification problem, using the movie reviews data set. We know the "negative" or "positive" labels for each of the movies. We'll set some of these aside for a test set and train our models on the remainder as a training set, using unigram presence or counts as the features. Then we'll evaluate the predictions quantitatively as well as look at some ways to interpret what the models tell us.

We'll start with Naive Bayes, move to logistic regression and its ridge and LASSO variants, then support vector machines and finally random forests. We'll also combine the models to examine an ensemble prediction.

We'll use these packages (install if necessary):

```{r}
library(dplyr)
library(quanteda)
library(quanteda.textmodels)
library(caret)
```

We'll start with the example given in the `quanteda` documentation. Read in the Pang and Lee dataset of 2000 movie reviews. (This appears to be the same 2000 reviews you used in the dictionary exercise, but in a different order.)

```{r}
corpus <- data_corpus_moviereviews
summary(corpus,5)
```

Shuffle the rows to randomize the order.
```{r}
set.seed(1234)
id_train <- sample(1:2000,1500, replace=F)
head(id_train, 10)
```

Use the 1500 for a training set and the other 500 as your test set. Create dfms for each.
```{r}
docvars(corpus, "id_numeric") <- 1:ndoc(corpus)

dfmat_train <- corpus_subset(corpus, id_numeric %in% id_train) %>% tokens() %>% dfm() #%>% dfm_weight(scheme="boolean")

dfmat_test <- corpus_subset(corpus, !(id_numeric %in% id_train)) %>% tokens %>% dfm() #%>% dfm_weight(scheme="boolean")
```

## Naive Bayes

Naive Bayes is a built in model for quanteda, so it's easy to use:

```{r}
sentmod.nb <- textmodel_nb(dfmat_train, docvars(dfmat_train, "sentiment"), distribution = "Bernoulli")
summary(sentmod.nb)
```

Use the dfm_match command to limit dfmat_test to features (words) that appeared in the training data:
```{r}
dfmat_matched <- dfm_match(dfmat_test, features=featnames(dfmat_train))
```

How did we do? Let's look at a "confusion" matrix.
```{r}
actual_class <- docvars(dfmat_matched, "sentiment")
predicted_class <- predict(sentmod.nb, newdata=dfmat_matched)
tab_class <- table(actual_class,predicted_class)
tab_class
```

Not bad, considering. Let's put some numbers on that:
```{r}
caret::confusionMatrix(tab_class, mode="everything", positive="pos")
```

Given the balance in the data among negatives and positives, "Accuracy" isn't a bad place to start. Here we have Accuracy of 81.6%.

The most commonly reporte alternatives I see are "F1" (the harmonic mean of Precision and Recall, see below) and "AUC" ("area under the curve," referring to the "ROC" or "Receiver Operating Characteristic" curve, discussed in the ROC curve section at the end of this notebook).

If the classes are imbalanced, accuracy becomes a less useful indicator. If, for example, you have 95% positive cases and 5% negative cases in your test set, a model which simply predicts "positive" without even looking at the data would get an accuracy of 0.95. 

Similarly, if you *care* about positive and negative cases differently -- say, your classifier is a Covid test or a jury decision, or negative is a "null hypothesis" -- you will care about precision/recall or specificity/sensitivity or TPR/FPR/TNR/FNR or PDR/NDR. 

"Sensitivity" and "recall" and "true positive rate" (TPR) and "hit rate" and "probability of detection" are all the same thing -- the probability that a positive case is classified as a positive. That is, the extent to which the classifier captures the positive cases it should capture. This is also the same as 1-FNR, the false negative rate. If negative cases are the "null," then FNR is the probability of "Type II error." 

"Specificity" or "selectivity" is the "true negative rate," the probability that negative cases are classified as negative. This is also 1-FPR, where FPR is the "false positive rate," the probability of "false alarm" -- classifying a negative case as positive -- also known as a "Type I error."

"Precision" or "positive detection rate" (PDR) is the probability that a case classified as positive is actually positive. "Negative detection rate" (NDR) is the probability that a case classified as negative is actually negative.

Let's wave our hands at validation by doing some sniff tests. What are the most positive and negative words?

The conditional posterior estimates, the "estimated feature scores" displayed by the summary command above, are stored in the model object's `param` attribute as a $K \times V$ matrix, where $K$ is the number of classes and $V$ is the vocabulary size. Here are the first 10 columns.

```{r}
sentmod.nb$param[,1:10]
```

These are the posterior probability that a word will appear in a document of a given class, $p(w|k)$. We want the posterior probability that a document is of a given class, given it contains that word, $p(k|w)$. This is obtained by Bayes Theorem: $p(\text{pos}|w) = \frac{p(w|\text{pos})}{p(w|\text{pos}) + p(w|\text{neg})}$

```{r}
#Most positive words
sort(sentmod.nb$param["pos",]/colSums(sentmod.nb$param),dec=T)[1:20]
```

There's reasonable stuff there: "outstanding", "seamless", "lovingly", "flawless". There's also some evidence of overfitting: "spielberg's", "winslet", "gattaca", "mulan". We'll see support for the overfitting conclusion below.

```{r}
#Most negative words
sort(sentmod.nb$param["neg",]/colSums(sentmod.nb$param),dec=T)[1:20]
```

Let's get a birds-eye view.
```{r, fig.width=7, fig.height=6}

post.pos <- sentmod.nb$param["pos",]/colSums(sentmod.nb$param)
# Plot weights
plot(colSums(dfmat_train),post.pos, pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Posterior Probabilities, Naive Bayes Classifier, IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances")
text(colSums(dfmat_train),post.pos, colnames(dfmat_train),pos=4,cex=5*abs(.5-post.pos), col=rgb(0,0,0,1.5*abs(.5-post.pos)))
```

Look a little closer at the negative.
```{r, fig.width=7, fig.height=6}
# Plot weights
plot(colSums(dfmat_train),post.pos, pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Posterior Probabilities, Naive Bayes Classifier, IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim=c(10,1000),ylim=c(0,.25))
text(colSums(dfmat_train),post.pos, colnames(dfmat_train),pos=4,cex=5*abs(.5-post.pos), col=rgb(0,0,0,1.5*abs(.5-post.pos)))
```

And a little more closely at the positive words:

```{r, fig.width=7, fig.height=6}
# Plot weights
plot(colSums(dfmat_train),post.pos, pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Posterior Probabilities, Naive Bayes Classifier, IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim=c(10,1000),ylim=c(0.75,1.0))
text(colSums(dfmat_train),post.pos, colnames(dfmat_train),pos=4,cex=5*abs(.5-post.pos), col=rgb(0,0,0,1.5*abs(.5-post.pos)))
```


Let's look a little more closely at the document predictions.

```{r}
predicted_prob <- predict(sentmod.nb, newdata=dfmat_matched, type="probability")

dim(predicted_prob)
head(predicted_prob)
summary(predicted_prob)
```

```{r, fig.height=3, fig.width=4}
plot(density(predicted_prob[,"pos"]), main = "Predicted Probabilities from Naive Bayes classifier", xlab="", ylab = "")
rug(predicted_prob[,"pos"])
```

You can see there one problem with the "naive" part of naive Bayes. By taking all of the features (words) as independent, it thinks it has seen far more information than it really has, and is therefore far more confident about its predictions than is warranted.

What's the most positive review in the test set according to this?

```{r}
# sort by *least negative* since near zero aren't rounded
sort.list(predicted_prob[,"neg"], dec=F)[1]
```

```{r}
id_test <- !((1:2000) %in% id_train)
writeLines(as.character(corpus)[id_test][440])
```
Looks like ``Amistad.'' A genuinely positive review, but note how many times "spielberg" is mentioned. The prediction is biased toward positive just because Spielberg had positive reviews in the training set. We may not want that behavior.

Note also that this is a very long review.

```{r}
# sort by *least neg* since near zero aren't rounded
sort.list(predicted_prob[,"pos"], dec=F)[1]
```

```{r}
writeLines(as.character(corpus)[id_test][211])
```

Schwarzenegger's ``End of Days.''

It also should be clear enough that more words means more votes, so longer documents are more clearly positive or negative. There's an argument for that. It also would underplay a review that read in it's entirety ``terrible,'' even though the review is 100% clear in its sentiment.

What is it most confused about?

```{r}
sort.list(abs(predicted_prob[,"pos"] - .5), dec=F)[1]
```

```{r}
predicted_prob[400,]
```

So ... the model says 52% chance negative, 48% positive.

```{r}
writeLines(as.character(corpus)[id_test][400])
```

A negative review of "Mafia!" a spoof movie I'd never heard of. Satire, parody, sarcasm, and similar are notoriously difficult to correctly classify, so perhaps that's what happened here.

Let's look at a clear mistake. 
```{r}
sort.list(predicted_prob[1:250,"neg"],dec=F)[1]
```
```{r}
predicted_prob[196,]
```
So ... the model says *DEFINITELY* positive.

```{r}
writeLines(as.character(corpus)[id_test][196])
```

Aha! A clearly negative review of "Saving Private Ryan."

This is at least partly an "overfitting" mistake. It probably learned other "Saving Private Ryan" or "Spielberg movies" words -- it looks like `spielberg's` was number #3 on our list above -- and learned that "reviews that talk about Saving Private Ryan are probably positive."

Below, I'll give brief examples of some other classification models for this data.

## Logistic regression, ridge regression, LASSO, and elasticnet

We'll look at three (well really only two) variants of the relatively straightforward regularized logistic regression model. These add a "regularization penalty" to the usual loss function of logistic regression. There are two basic types of regularization penalty: L1 and L2. 

The L2 penalty imposes loss proportional to the sum of the squared coefficients. That is, it biases the model toward finding smaller coefficients. In Bayesian terms, this is also equivalent to using a Gaussian prior. Logistic regression with an L2 penalty is known as "ridge regression."

The L1 penalty imposes loss proportional to the sum of the absolute values of the coefficients. That is, it biases the model toward finding more coefficients equal exactly to zero. It has a feature selection effect. In Bayesian terms, this is also equivalent to using a Laplace prior. Logistic regression with an L1 penalty is known as the "LASSO" (Least Absolute Shrinkage and Selection Operator).

How much of a penalty is applied is determined by a weight, $\lambda$. This is a "hyperparameter" that can be tuned during the training stage.

The elastic net adds *both* penalties, distributing the total regularization cost according to a second hyperparameter, $\alpha$, and estimates the optimal values for both.


```{r}
library(glmnet)
library(doMC)
```

### Ridge regression (Logistic with L2-regularization)

```{r}
registerDoMC(cores=2) # parallelize to speed up
sentmod.ridge <- cv.glmnet(x=dfmat_train,
                   y=docvars(dfmat_train)$sentiment,
                   family="binomial", 
                   alpha=0,  # alpha = 0: ridge regression
                   nfolds=5, # 5-fold cross-validation
                   parallel=TRUE, 
                   intercept=TRUE,
                   type.measure="class")
plot(sentmod.ridge)
```

This shows classification error as $\lambda$ (the total weight of the regularization penalty) is increased from 0.
The minimum error is at the leftmost dotted line, about $\log(\lambda) \approx 3$. This value is stored in `lambda.min`.

```{r}
# actual_class <- docvars(dfmat_matched, "sentiment")
predicted_value.ridge <- predict(sentmod.ridge, newx=dfmat_matched,s="lambda.min")[,1]
predicted_class.ridge <- rep(NA,length(predicted_value.ridge))
predicted_class.ridge[predicted_value.ridge>0] <- "pos"
predicted_class.ridge[predicted_value.ridge<0] <- "neg"
tab_class.ridge <- table(actual_class,predicted_class.ridge)
tab_class.ridge
```


```{r}
confusionMatrix(tab_class.ridge, mode="everything", positive="pos")
```

Accuracy of .818, similar to the others. The misses are a little more evenly distributed across classes, with it being slightly more successful in identifying positive reviews and slightly less successful in identifying negative reviews.


At first blush, the coefficients should tell us what the model learned:

```{r, fig.width=7, fig.height=6}
plot(colSums(dfmat_train),coef(sentmod.ridge)[-1,1], pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Ridge Regression Coefficients, IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
#text(colSums(dfmat_train),coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=200*abs(coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,.3))
text(colSums(dfmat_train),coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=200*abs(coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,75*abs(coef(sentmod.ridge)[-1,1]))) # inexplicable error when knitting of alpha not in [0,1] range when it absolutely is
```

That's both confusing and misleading since the variance of coefficients is largest with the most obscure terms. (And, for plotting, the 40,000+ features include some very long one-off "tokens" that overlap with more common ones, e.g., "boy-drinks-entire-bottle-of-shampoo-and-may-or-may-not-get-girl-back," "props-strategically-positioned-between-naked-actors-and-camera," and "____________________________________________")


With this model, it would be more informative to look at which coefficients have the most impact when making a prediction, by having larger coefficients *and* occurring more, or alternatively to look at which coefficients we are most certain of, downweighting by the inherent error. The impact will be proportional to $log(n_w)$ and the error will be roughly proportional to $1/sqrt(n_w)$.

So, impact:

```{r, fig.width=7, fig.height=6}
plot(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Ridge Regression Coefficients (Impact Weighted), IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
text(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=50*abs(log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,25*abs(log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1])))
#text(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=50*abs(log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,.3))
```

Most positive and negative features by impact:

```{r}
sort(log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1],dec=T)[1:20]
```

```{r}
sort(log(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1],dec=F)[1:20]
```

Regularization and cross-validation bought us a lot more general -- less overfit -- model than we saw with Naive Bayes.

Alternatively, by certainty:

```{r, fig.width=7, fig.height=6}
plot(colSums(dfmat_train),sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Ridge Regression Coefficients (Error Weighted), IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
#text(colSums(dfmat_train),sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=30*abs(sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,.3))
text(colSums(dfmat_train),sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1], colnames(dfmat_train),pos=4,cex=30*abs(sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1]), col=rgb(0,0,0,10*abs(sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1])))
```

Most positive and negative terms:

```{r}
sort(sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1],dec=T)[1:20]
```

```{r}
sort(sqrt(colSums(dfmat_train))*coef(sentmod.ridge)[-1,1],dec=F)[1:20]
```

This view implies some would-be "stop words" are important, and these seem to make sense on inspection. For example, "as" is indicative of phrases in positive reviews comparing movies to well-known and well-liked movies, e.g., "as good as." There's not a parallel "as bad as" that is as common in negative reviews.

### LASSO (Logistic with L1-regularization)

Ridge regression gives you a coefficient for every feature. At the other extreme, we can use the LASSO to get some feature selection.

```{r}
registerDoMC(cores=2) # parallelize to speed up
sentmod.lasso <- cv.glmnet(x=dfmat_train,
                   y=docvars(dfmat_train)$sentiment,
                   family="binomial", 
                   alpha=1,  # alpha = 1: LASSO
                   nfolds=5, # 5-fold cross-validation
                   parallel=TRUE, 
                   intercept=TRUE,
                   type.measure="class")
plot(sentmod.lasso)
```


```{r}
# actual_class <- docvars(dfmat_matched, "Sentiment")
predicted_value.lasso <- predict(sentmod.lasso, newx=dfmat_matched,s="lambda.min")[,1]
predicted_class.lasso <- rep(NA,length(predicted_value.lasso))
predicted_class.lasso[predicted_value.lasso>0] <- "pos"
predicted_class.lasso[predicted_value.lasso<0] <- "neg"
tab_class.lasso <- table(actual_class,predicted_class.lasso)
tab_class.lasso
```

```{r}
confusionMatrix(tab_class.lasso, mode="everything", positive="pos")
```

Slightly better accuracy of .82. The pattern of misses goes further in the other direction from Naive Bayes, overpredicting positive reviews.


```{r, fig.width=7, fig.height=6}
plot(colSums(dfmat_train),coef(sentmod.lasso)[-1,1], pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="LASSO Coefficients, IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
#text(colSums(dfmat_train),coef(sentmod.lasso)[-1,1], colnames(dfmat_train),pos=4,cex=2*abs(coef(sentmod.lasso)[-1,1]), col=rgb(0,0,0,.3))
text(colSums(dfmat_train),coef(sentmod.lasso)[-1,1], colnames(dfmat_train),pos=4,cex=2*abs(coef(sentmod.lasso)[-1,1]), col=rgb(0,0,0,1*abs(coef(sentmod.lasso)[-1,1])))
```

As we want when we run the LASSO, the vast majority of our coefficients are zero ... most features have no influence on the predictions.

It's less necessary but let's look at impact:

```{r, fig.width=7, fig.height=6}
plot(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1], pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="LASSO Coefficients (Impact Weighted), IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
#text(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1], colnames(dfmat_train),pos=4,cex=.8*abs(log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1]), col=rgb(0,0,0,.3))
text(colSums(dfmat_train),log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1], colnames(dfmat_train),pos=4,cex=.8*abs(log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1]), col=rgb(0,0,0,.25*abs(log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1])))
```

Most positive and negative features by impact:

```{r}
sort(log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1],dec=T)[1:20]
```

Interestingly, there are a few there that would be negative indicators in most sentiment dictionaries, like "flaws" and "war".

```{r}
sort(log(colSums(dfmat_train))*coef(sentmod.lasso)[-1,1],dec=F)[1:20]
```

Both lists also have words that indicate a transition from a particular negative or positive aspect, followed by a holistic sentiment in the opposite direction. "The pace dragged at times, but *overall* it is an astonishing act of filmmaking." "The performances are tremendous, but *unfortunately* the poor writing makes this movie fall flat." 

### Elastic net

The elastic net estimates not just $\lambda$ (the overall amount of regularization) but also $\alpha$ (the relative weight of the L1 loss relative to the L2 loss). In R, this can also be done with the `glmnet` package. "I leave that as an exercise."

## A first ensemble

We've got three sets of predictions now, so why don't we try a simple ensemble in which our prediction for each review is based on a majority vote of the three. Sort of like a Rotten Tomatoes rating. They each learned slightly different things, so perhaps the whole is better than its parts.

```{r}
predicted_class.ensemble3 <- rep("neg",length(actual_class))
num_predicted_pos3 <- 1*(predicted_class=="pos") + 1*(predicted_class.ridge=="pos") + 1*(predicted_class.lasso=="pos")
predicted_class.ensemble3[num_predicted_pos3>1] <- "pos"
tab_class.ensemble3 <- table(actual_class,predicted_class.ensemble3)
tab_class.ensemble3
```


```{r}
confusionMatrix(tab_class.ensemble3, mode="everything", positive="pos")
```

Hey, that is better! Accuracy 83.6%!



## Support vector machine

Without explaining SVM at all, let's try a simple one.

```{r}
library(e1071)
```

```{r}
sentmod.svm <- svm(x=dfmat_train,
                   y=as.factor(docvars(dfmat_train)$sentiment),
                   kernel="linear", 
                   cost=10,  # arbitrary regularization cost
                   probability=TRUE)
```

Ideally, we would tune the cost parameter via cross-validation or similar, as we did with $\lambda$ above.

```{r}
# actual_class <- docvars(dfmat_matched, "Sentiment")
predicted_class.svm <- predict(sentmod.svm, newdata=dfmat_matched)
tab_class.svm <- table(actual_class,predicted_class.svm)
tab_class.svm
```


```{r}
confusionMatrix(tab_class.svm, mode="everything", positive="pos")
```

That's actually a bit better than the others, individually if not combined, with accuracy of .828, and a bias toward overpredicting positives.



For a linear kernel, we can back out interpretable coefficients. This is not true with nonlinear kernels such as the "radial basis function."

```{r}
beta.svm <- drop(t(sentmod.svm$coefs)%*%dfmat_train[sentmod.svm$index,])
```

(Note the signs are reversed from our expected pos-neg.)

```{r,fig.width=7,fig.height=6}
plot(colSums(dfmat_train),-beta.svm, pch=19, col=rgb(0,0,0,.3), cex=.5, log="x", main="Support Vector Machine Coefficients (Linear Kernel), IMDB", ylab="<--- Negative Reviews --- Positive Reviews --->", xlab="Total Appearances", xlim = c(1,50000))
text(colSums(dfmat_train),-beta.svm, colnames(dfmat_train),pos=4,cex=10*abs(beta.svm), col=rgb(0,0,0,5*abs(beta.svm)))
```

```{r}
sort(-beta.svm,dec=T)[1:20]
```

```{r}
sort(-beta.svm,dec=F)[1:20]
```

Looks a bit overfit to me and I would probably increase the regularization cost in further iterations.

## Random Forests

```{r}
library(randomForest)
```

Random forests is a very computationally intensive algorithm, so I will cut the number of features way way down just so this can run in a reasonable amount of time.

```{r}
dfmat.rf <- corpus %>%
  tokens() %>%
  dfm() %>%
  dfm_trim(min_docfreq=50,max_docfreq=300,verbose=TRUE)
```

```{r}
dfmatrix.rf <- as.matrix(dfmat.rf)
```

```{r}
set.seed(1234)
sentmod.rf <- randomForest(dfmatrix.rf[id_train,], 
                   y=as.factor(docvars(dfmat.rf)$sentiment)[id_train],
                   xtest=dfmatrix.rf[id_test,],
                   ytest=as.factor(docvars(dfmat.rf)$sentiment)[id_test],
                   importance=TRUE,
                   mtry=20,
                   ntree=100
                   )
#sentmod.rf
```

```{r}
predicted_class.rf <- sentmod.rf$test[['predicted']]
tab_class.rf <- table(actual_class,predicted_class.rf)
confusionMatrix(tab_class.rf, mode="everything", positive="pos")
```

That did a bit worse -- Accuracy .8 -- but we did give it considerably less information.

Getting marginal effects from a random forest model requires more finesse than I'm willing to apply here. We can get the "importance" of the different features, but this alone does not tell us in what direction the feature pushes the predictions.

```{r,fig.height=7,fig.width=6}
varImpPlot(sentmod.rf)
```

Some usual suspects there, but we need our brains to fill in which ones are positive and negative. Some are ambiguous ("town") and some we have seen are subtle ("overall") or likely to be misleading ("war").

## Another ensemble

Now we've got five, so let's ensemble those.

```{r}
  predicted_class.ensemble5 <- rep("neg",length(actual_class))
num_predicted_pos5 <- 1*(predicted_class=="pos") + 1*(predicted_class.ridge=="pos") + 1*(predicted_class.lasso=="pos") + 
  1*(predicted_class.svm=="pos") + 
  1*(predicted_class.rf=="pos")
predicted_class.ensemble5[num_predicted_pos5>2] <- "pos"
tab_class.ensemble5 <- table(actual_class,predicted_class.ensemble5)
tab_class.ensemble5
confusionMatrix(tab_class.ensemble5,mode="everything", positive="pos")
```

And, like magic, now we're up over 85% Accuracy in the test set.

Note that most of these agree with each other on about 80% of cases (Naive Bayes and Ridge Regression agree on 93.4% of cases, and every other pair is between 78.8% and 81.0% agreement). The ensembling works *because* of this general disagreement, not *in spite* of it. If the models made the same mistakes -- agreed with each other -- there would be no gain from ensembling. Ensembling works best when you have different sorts of models that learn different sorts of things.

## Compare to a dictionary

```{r}
library(quanteda.sentiment)
```

```{r}
polarity.lsd <- textstat_polarity(dfmat_test, dictionary = data_dictionary_LSD2015, fun=sent_relpropdiff)
predicted_class.lsd <- rep("neg",nrow(polarity.lsd))
predicted_class.lsd[polarity.lsd$sentiment>0] <- "pos"
tab_class.lsd <- table(actual_class,predicted_class.lsd)
confusionMatrix(tab_class.lsd,mode="everything", positive="pos")
```

Accuracy 67.8%. Yeah. Better than flipping a coin, but ...

## ROC Curves

Another common approach to evaluating classifiers is the ROC curve, which looks in more depth at the underlying scores / probability statements that indicate how sure the classifier is. There are multiple libraries that provide ROC curve functionality. We'll use `ROCit`.

```{r}
library(ROCit)
```

For these, we will use the underlying "score" emitted from the classifiers -- the predicted probabilities, or the dictionary polarity score.

There are multiple versions, but a typical one plots "Sensitivity" (True Positive Rate) vs. "Specificity" (1 - False Positive Rate) as the threshold for predicting class labels moves from the lowest score received to the highest.



```{r}
roc_empirical.lsd <- rocit(score = polarity.lsd$sentiment, class = actual_class,negref = "neg")
roc_empirical.nb <- rocit(score =predict(sentmod.nb, newdata=dfmat_matched, type="probability")[,"pos"] , class = actual_class,negref = "neg")

roc_empirical.ridge <- rocit(score = predicted_value.ridge, class = actual_class,negref = "neg")
roc_empirical.lasso <- rocit(score = predicted_value.lasso, class = actual_class,negref = "neg")
predicted_prob.svm <- attr(predict(sentmod.svm, newdata=dfmat_matched, probability=TRUE),"probabilities")[,"pos"]
roc_empirical.svm <- rocit(score = predicted_prob.svm, class = actual_class,negref = "neg")
roc_empirical.rf <- rocit(score = sentmod.rf$test$votes[,2], class = actual_class,negref = "neg")
```




```{r,fig.height=4,fig.width=4}
plot(c(0,1),c(0,1), type="n",main = "ROC Curves",ylab = "Sensitivity (TPR)", xlab = "1-Specificity (FPR)")
lines(c(0,1),c(0,1), col="grey", lty=2)
lines(roc_empirical.lsd$FPR,roc_empirical.lsd$TPR,col="black",lwd=2)
lines(roc_empirical.nb$FPR,roc_empirical.nb$TPR,col="red",lwd=2)
lines(roc_empirical.ridge$FPR,roc_empirical.ridge$TPR,col="blue",lwd=2)
lines(roc_empirical.lasso$FPR,roc_empirical.lasso$TPR,col="green",lwd=2)
lines(roc_empirical.svm$FPR,roc_empirical.svm$TPR,col="orange",lwd=2)
lines(roc_empirical.rf$FPR,roc_empirical.rf$TPR,col="purple",lwd=2)
legend(.6,.2,legend=c("Lexicoder","Naive Bayes","Ridge Regression","LASSO","Support Vector Machine","Random Forests"), col=c("black","red","blue","green","orange","purple"),lty=1,lwd=2)
```

A good ROC curve steps up steeply from 0, nearing 1 quickly, with its "elbow" as close to the upper left corner as possible.

A random classifier will have the ROC curve shown in gray.

This shows clearly that the analysis using the Lexicoder dictionary is by far the poorest. Among the other five, SVM and LASSO are slightly better, followed by Naive Bayes and Ridge Regression. Random Forest is the poorest of these five, but as previously noted we used much less information in building that model because of computation time.

This can be summarized with the AUC (Area Under Curve) statistic:

```{r}
print(paste("AUC, Lexicoder: ",sprintf("%1.3f",roc_empirical.lsd$AUC)))
print(paste("AUC, Naive Bayes: ",sprintf("%1.3f",roc_empirical.nb$AUC)))
print(paste("AUC, Ridge Regression: ",sprintf("%1.3f",roc_empirical.ridge$AUC)))
print(paste("AUC, LASSO: ",sprintf("%1.3f",roc_empirical.lasso$AUC)))
print(paste("AUC, SVM: ",sprintf("%1.3f",roc_empirical.svm$AUC)))
print(paste("AUC, Random Forests: ",sprintf("%1.3f",roc_empirical.rf$AUC)))
```