The Goal

We have been working with multinomial regression models, which are models we can use when our response variable is categorical with more than 2 levels. We work with such data a lot when we are dealing with text, so multinomial is a model we encounter often.

Our client for today is a cell phone carrier who is trying to determine what customers are saying about different cell phones they carry. They head over to a common source of reviews online. However, a glitch in the system often means phone reviews are labelled under the wrong brand (this happens A LOT in Amazon and other review systems!). In other words, while looking at reviews for Apple phones, there are actually reviews for Samsung mixed in, and vice versa. This makes it very challenging for companies to find the insight that they need when the reviews are not all from the same product!

They provide us a data set called train with 1000 phones that are labelled with their brand. They ask us to use those phones to build a model to (1) predict which brand a phone is and (2) describe any key traits in reviews that might suggest a phone is from a certain brand.

They then also provide us a data set called test with 500 phones to use to assess the predictive abilities of our model.

You can use the following code to load the data.

# Load the libraries
library(dplyr)
library(ggplot2)
library(nnet)
library(broom)

# Load the data
train <- read.csv("https://www.dropbox.com/scl/fi/y9srlcqu36qh4x2gxx1mf/CellPhonetrain.csv?rlkey=qc1bgktk8kdg8xo0uwhx9gbiv&st=bp2sdljz&dl=1")
test  <- read.csv("https://www.dropbox.com/scl/fi/i875r45iwefli11b8bhbq/CellPhonetest.csv?rlkey=ux8hcgfrb7dagm1lirmhszlxz&st=4gb6yu3q&dl=1")

If you look at the train dataset you have loaded, the dataset should have 1000 rows and 119 columns. The first column is our Y variable brand, and the remaining 118 columns are features about the text created from the LIWC.

As a reminder, the LIWC is a tool that allows us access to multiple lexicons that we can use to create features about text data. For instance, how much of the review is written in past tense? How much involves discussions of economics? How long are the sentences, and what is the tone like? The list of possible features is extensive!

If you are curious about what each feature represents, you can look at the list on page 12 of this document. This details each of the features created by the LIWC and gives examples of the types of words in each of the lexicons it uses to create these features.

Model 1: Perception

Our client thinks that one feature that might be important is the percent of words in a review that are related to perception. These are words about how a phone feels, sounds, looks, etc. In other words, these are words about how a person physically interacts with the device.

• Percepetion = the percent of words in a text that have to do with the sensory experiences (hearing, feeling, touching, smelling, etc.)

They ask us to build a model using this as a feature.

Recall that for multinomial regression, we use:

# Load the library
library(nnet)
library(broom)

# Build the model
Model1_Perception <- multinom( Y ~ X , data = dataset, trace = FALSE)

# Print out the coefficients
knitr::kable( tidy(Model1_Perception) ) 

The kind of interpretation work we have done in the previous two questions is useful. However, especially when we are working with a client, it is extremely valuable to be able to visualize the relationships that our model is describing. In other words, we would like to be able to make a plot that shows how the probability of being written about each brand changes as the percent of perception words changes.

We can do this in R, but the code involved is a little tricky, so I built a function that will do it for us. To teach R the function, copy and paste the following into a chunk and press play. DO NOT CHANGE ANYTHING ABOUT THIS FUNCTION.

plot_multnomial <- function(data, model, Y, X, Xname, xlab){
  
  if(length(coefficients(model)[1,])==2){
  
  predictions <- predict( model, type ="probs")
  
  data[, Y] <-as.factor(data[,Y])
  
  J = length(levels(data[,Y]))
  
  holder <- data.frame("group" = NA, "x" = NA, "Probability"=NA)
  a  = 1 
  for(i in 1:nrow(data)){
    for(j in 1:J){
      holder[a,1] <- levels(data[,Y])[j]
      holder[a,2] <- data[i,X]
      holder[a,3] <- predictions[i,j]
      a = a + 1
    }
  }
  
  g1 <- ggplot( holder, aes(x, Probability, col= group)) + 
    geom_line() + 
    labs(x = xlab)
  }
  
  if(length(coefficients(model)[1,])>2){
    
    suppressMessages(library(pdp))
    suppressMessages(library(ggplot2))
    
    data[ , Y] <-as.factor(data[,Y])
    
    toplot1 <- data.frame(partial(model, Xname, which.class = 1, plot = FALSE, prob = TRUE))
    toplot2 <- data.frame(partial(model, Xname, which.class = 2, plot = FALSE, prob = TRUE))
    toplot3 <- data.frame(partial(model, Xname, which.class = 3, plot = FALSE, prob = TRUE))
    
    toplot1$response <- levels(data[,Y])[1]
    toplot2$response <- levels(data[,Y])[2]
    toplot3$response <- levels(data[,Y])[3]
    
    hold <- rbind(toplot1,toplot2,toplot3)
    hold <-data.frame(hold)
    hold$response <- as.factor(hold$response)

    g1 <- ggplot(hold, aes(x = hold[,1], yhat, col = response)) + geom_line() + 
      labs(x = xlab, y = "Probability")
    
  }
  
  g1
  
  
  
}

To use this function, adapt the following:

plot_multnomial(data, model, Y, X, Xname, xlab)

• data: replace with the name of your data set.
• model: replace with the name of your model.
• Y: replace with the number of the column in your data set when \(Y\) is.
• X: replace with the number of the column in your data set when \(X\) is. For perception, this is 96.
• Xname: replace with the name of your \(X\) column. Make sure you put it in quotes, for example: "Perception"
• xlab: replace with the label you want for the x axis label. Make sure you put it in quotes, for example: "The x axis label".

In addition to describing the relationships in the data, if we want to find predicted probabilities from a multinomial regression model for a row \(i\) in the data set, we can use:

fitted(Model1_Perception)[i,]

Model 2: Full Multinomial Model

Recall that the client asked us to build a model that could predict which brand of phone each review is about. So far, we have only been working with one feature, but when we move into prediction it makes a lot more sense to use all the features we have access too!

We have 118 available to us, so let’s use them all. Adding every feature we have access to is called making a full model.

# Train the full model
Model2_Full <- multinom( brand ~ . , data = train, trace = FALSE)

This model will have a LOT of coefficients, so it won’t be super useful to look at a table of coefficients. However, it should be useful for their other goal of predicting which brand each review is about.

To make predictions on the test using multinomial regression, we use the code:

# Make predictions 
predictions <- predict(Model2_Full, type = "class", newdata = test)

# Create a confusion matrix
holder <- table(  predictions, test$brand )
rownames(holder) <- paste("Predicted", rownames(holder))
colnames(holder) <- paste("True", colnames(holder))
knitr::kable( holder )

We can also use the pretty confusion matrix code we did with logistic regression!!

Naive Bayes

Because the client asked us to do prediction, it makes sense to try more than one approach. In practice, we often try multiple predictive approaches. Now that we have tried multinomial, let’s try Naive Bayes!

Let’s start off with just two features to remind ourselves how this algorithm works, and then we will see how to run the whole algorithm using R.

Let \(X_1\) = 1 if the positive sentiment score (emo_pos) is greater than 5, and \(X_2\) = 1 if the review contains at least one exclamation point.

X1 <- ifelse( train$emo_pos > 5, 1, 0 )
X2 <- ifelse( train$Exclam >0, 1, 0)

Let’s use Naive Bayes to find the probability of being about an Apple phone if the positive sentiment score is 10 and there are 3 exclamation points.

Recall that for Naive Bayes, we find our probabilities using:

\[P(y_i = A ~ | ~ x_{i1} ~and~ x_{i2} ) = \frac{P(y_i = A) \times P(x_{i1}~|~A) \times P(x_{i2}~|~A)}{P(x_{i1} ~and~ x_{i2})}\]

I know this is long, but it is important to do a few times to make sure we feel comfortable with what Naive Bayes is actually doing. In practice, it is MUCH more efficient to use R to help with the computations.

To run Naive Bayes in R, we need the naivebayes library:

library(naivebayes)

Remember: If your computer ever tells you a library is not found, this means you need to install the package from Tools.

Since we are using computing now, let’s go back to using all the features rather than just \(X_1\) and \(X_2\). The code we use to make predictions with Naive Bayes is then:

# Train Naive Bayes 
nb <- naive_bayes( Y ~ . , data= ... )
 
# Make predictions 
yhat <- predict(nb , type = "class", newdata = ...)

# Create a confusion matrix
holder <- table(  yhat, test$brand )
rownames(holder) <- paste("Predicted", rownames(holder))
colnames(holder) <- paste("True", colnames(holder))
knitr::kable( holder )

Comparing Models

At this point, we have predictions from multinomial regression:

# Create a confusion matrix
holder <- table(  predictions, test$brand )
rownames(holder) <- paste("Predicted", rownames(holder))
colnames(holder) <- paste("True", colnames(holder))
knitr::kable( holder )

and from Naive Bayes:

# Create a confusion matrix
holder <- table(  yhat, test$brand )
rownames(holder) <- paste("Predicted", rownames(holder))
colnames(holder) <- paste("True", colnames(holder))
knitr::kable( holder )

When it comes time to compare different predictive approaches, it is helpful if we have some metric, meaning some number, we can use to compare them. We could use the accuracy, but this tends to not work as well unless all our options for \(Y\) have roughly same number of rows.

In situations where the classes are not balanced, we reach for a tool called the Macro F1-score.

A standard F1-score is computed as

\[F1 = \frac{TP}{TP+.5(FP + FN)}\]

For example, the F1-score for Apple is computed as

\[F1_{Apple} = \frac{TP_{Apple}}{TP_{Apple}+.5(FP_{Apple} + FN_{Apple})}\]

where

• \(TP_{Apple}\) = the number of Apple reviews that are correctly identified
• \(FP_{Apple}\) = the number of reviews that are predicted to be Apple but are not
• \(FN_{Apple}\) = the number of reviews that are Apple but are predicted to be something else.

It is important to note that F1-scores are between 0 and 1, and the higher, the better. This means an F1-score of .9 indicates a strong predictive performance, while something like .2 does not. Beyond this, we do not intepret F1-scores.

To get the Macro F1-score, all we do is average together the F1-scores for each of the 3 brands. As a note, neither the Macro F1-score nor standard F1-scores are interpretable. All we do is compare them, and the higher, the better!

Like the F1-scores, Macro F1-scores are between 0 and 1, and higher values indicate better prediction.

When we are doing prediction, looking for a higher Macro F1-score can be a tool to determine which model is doing better at prediction overall.

We just saw that computing the Macro F1-score by hand can take a while. Luckily, the caret package has a built in.

library(caret) 

# Make sure R knows brand is categorical
test$brand <- factor(test$brand)

# Build the confusion matrix
CM <- confusionMatrix(data = predicted, reference = actual, mode = "everything")

# Pull out the F1-scores
class_f1_scores <- CM$byClass[, "F1"]

# Compute the Macro F1 Score
mean(class_f1_scores)

Before you submit!

Change the chunk at the top of your document from

knitr::opts_chunk$set(echo = TRUE)

to

knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)