rm(list=ls(all=TRUE))
getwd()
# setwd("C:/Users/mw/Dropbox (VOICES)/TOPIC MODEL")
setwd("C:/Users/luigi/Dropbox/TOPIC MODEL")
getwd()
library(quanteda)
library(readtext)
library(caTools)
library(caret)
library(car)
library(ggplot2)
library(dplyr)
library(reshape2)
library(glmnet)
#####################################################
# FIRST STEP: let's create the DfM for the training-set
#####################################################
# let's focus on MOVIE reviews (either positive or negative)
x <- read.csv("train_review2.csv", stringsAsFactors=FALSE)
str(x)
nrow(x)
table(x$Sentiment)
prop.table(table(x$Sentiment))
myCorpusTwitterTrain <- corpus(x)
Dfm_train <- dfm(myCorpusTwitterTrain , remove = c(stopwords("english")), remove_punct = TRUE, remove_numbers=TRUE,
tolower = TRUE, remove_symbols=TRUE, remove_separators=TRUE, remove_url = TRUE, split_hyphens = TRUE)
# Let's trim the dfm in order to keep only tokens that appear in at least 5% of the reviews
Dfm_train <- dfm_trim(Dfm_train , min_docfreq = 0.05, verbose=TRUE, docfreq_type = "prop")
topfeatures(Dfm_train , 20) # 20 top words
#####################################################
# SECOND STEP: let's create the DfM for the test-set
#####################################################
x10 <- read.csv("test_review2.csv", stringsAsFactors=FALSE)
str(x10)
nrow(x10)
myCorpusTwitterTest <- corpus(x10)
Dfm_test<- dfm(myCorpusTwitterTest , remove = c(stopwords("english")), remove_punct = TRUE, remove_numbers=TRUE,
tolower = TRUE, remove_symbols=TRUE, remove_separators=TRUE, remove_url = TRUE, split_hyphens = TRUE)
Dfm_test<- dfm_trim(Dfm_test, min_docfreq = 0.05, verbose=TRUE, docfreq_type = "prop")
topfeatures(Dfm_test , 20) # 20 top words
#####################################################
# THIRD STEP: Let's make the features identical between train and test-set by passing Dfm_train to dfm_match() as a pattern.
# after this step, we can "predict" by employing only the features included in the training-set!
# this moreover makes the DfM for the test-set smaller if you have a giant test-set compared
# to the training-set (i.e., of course this is not the case!)
#####################################################
setequal(featnames(Dfm_train), featnames(Dfm_test))
length(Dfm_test@Dimnames$features)
length(Dfm_train@Dimnames$features)
test_dfm <- dfm_match(Dfm_test, features = featnames(Dfm_train))
length(test_dfm@Dimnames$features)
setequal(featnames(Dfm_train), featnames(test_dfm ))
#####################################################
# FOURTH STEP: Let's convert the DfM into matrices for the ML algorithms to work
#####################################################
train <- as.matrix(Dfm_train)
test <- as.matrix(test_dfm)
#####################################################
# FIFHT STEP: let's estimate a ML Model
#####################################################
#####################################################
# Let's start with a Regularized regression - RIDGE model
#####################################################
# Ridge regression requires that alpha=0
# As a family here you identify "binomial" cause your output variable has 2 categories (negative, positive)
# if it were a 3 or more categories, you HAVE to use the family "multinomial"
system.time(ROLS <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train, family="binomial", alpha=0, trace.it=1))
# as you can see as lambda increases, the penalty becomes large and forces the coefficients to a value closer (but
# NOT exactly equal to) zero
plot(ROLS, xvar = "lambda")
# let's predict the test-set
system.time(predicted_rr <- predict(ROLS, test,type="class"))
# why you get 100,000 predictions? cause it is 1000 texts in the test-set * 100 Lambda estimated by default!
table(predicted_rr )
# So better always use cv.glmnet and extract the optimal lasso$lambda.min
set.seed(123)
system.time(ridge <- cv.glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train,
family="binomial", alpha=0, nfolds=10, type.measure="class", trace.it=1 ))
plot(ridge)
# The plot output illustrates the 10-fold CV minimum miss-classification error across the lambda values.
# As we constrain our coefficients with log of lambda>2 penalty, the miss-classification error rises considerably.
# The number at the top of the plot (1058) just refers to the number of variables (features in our case) in the model.
# Ridge regression does not force any variables to exactly zero so all features will remain in the model.
# The first and second vertical dashed lines represent the lambda value with the minimum miss-classification error and the largest
# lambra value within one standard error of the minimum miss-classification error
min(ridge $cvm) # minimum miss-classification error
ridge $lambda.min # lambda for this minimum
log(ridge $lambda.min )
ridge $cvm[ridge $lambda == ridge $lambda.1se] # 1 st.error of minimum miss-classification error
ridge $lambda.1se # lambda for this error
log(ridge$lambda.1se)
lse <- as.numeric(ridge$lambda.1se) # let's save the lambda in this latter case
lse
# Once you have identified your preferred lambda, you can simply use the command "predict" to predict the same model on a new data set.
# The only caveat is that you need to supply to predict an s parameter with the preferred models lambda value
# Let's go back to our previous ROLS object
plot(ROLS, xvar = "lambda")
abline(v = log(ridge $lambda.min), col = "red", lty = "dashed")
abline(v = log(ridge $lambda.1se), col = "blue", lty = "dashed")
# these are the results when you use lambda.min
system.time(predicted_glment <- predict(ROLS, test, s = ridge $lambda.min, type="class"))
table(predicted_glment)
prop.table(table(predicted_glment))
# same results you get if you predict with the CV for ridge above
system.time(predicted_rr2 <- predict(ridge , s= ridge $lambda.min, test,type="class"))
table(predicted_rr2 )
prop.table(table(predicted_rr2))
# these are the results when you use lambda.1se
system.time(predicted_glment2 <- predict(ROLS, test, s = lse, type="class"))
table(predicted_glment2)
prop.table(table(predicted_glment2))
# Why should you ever want to use lambda.lse rather than lambda.min?
# The advantage of identifying the lambda with a miss-classification error within one standard error from the minimum
# becomes more obvious with the lasso and elastic net models (see below!)
# CHECK FOR COEFFICIENTS of the model with lambda.min. Why caring about them???
# Often, ML models are considered black boxes due to their complex inner-workings. However, because of their complexity,
# they are typically more accurate for predicting nonlinear, faint, or rare phenomena. Unfortunately, more accuracy often comes at the
# expense of interpretability, and interpretability is crucial. It is not enough to identify a machine learning model that optimizes
# predictive performance; understanding and trusting model results is a hallmark of good science!
# Luckily, several advancements have been made to aid in interpreting ML models over the years.
# Interpreting ML models is an emerging field that has become known as interpretable machine learning (IML).
# Approaches to model interpretability can be broadly categorized as providing global or local explanations.
# Global interpretations help us understand the inputs and their entire modeled relationship with the prediction target.
# Local interpretations help us understand model predictions for a single row of data or a group of similar rows.
# Global interpretability is about understanding how the model makes predictions, based on a holistic view of its features and how they
# influence the underlying model structure. It answers questions regarding which features are relatively influential as well as how these
# features influence the response variable. In this latter case, we assess the magnitude of a variable’s relationship with the response
# as compared to other variables used in the model.
# In text analytics (given that we are dealining with huge DfM wherein the features per-se are not the main focus of our interest)
# we are mainly interested in global interpretations; however if you use ML for other aims, local interpretations become VERY important!!!
# So let's do some Global Interpretation for our RR!
newROLS <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train,
family="binomial", alpha=0, trace.it=1, lambda=ridge $lambda.min)
str(newROLS)
# let's create a dataframe with the features and their coefficients
coeff_df <- as.data.frame(Dfm_train@Dimnames$features)
str(coeff_df)
names(coeff_df)[1] <- "Feature"
coeff_df$coeff <- newROLS$beta@x
str(coeff_df)
# From here, let's replicate what we did for RF earlier
predicted_rr <- predict(newROLS, train, type="class")
table(predicted_rr )
Dfm_train@docvars$predRR <- ifelse(predicted_rr=="neg",0,1)
table(Dfm_train@docvars$predRR)
# NOTE: no perfect prediction of the training-set
table(Dfm_train@docvars$Sentiment, Dfm_train@docvars$predRR)
# adding sign [if 0/1 according to the content of the tweet - negative or positive]
sums <- list()
for (v in 0:1){
sums[[v+1]] <- colSums(train[Dfm_train@docvars[,"predRR"]==v,])
}
sums <- do.call(cbind, sums)
sign <- apply(sums, 1, which.max)
names <- dimnames(train)[[2]]
str(names)
df <- data.frame(
Feature = names,
sign = sign-1,
stringsAsFactors=F)
str(df)
str(coeff_df)
importanceRR <- merge(coeff_df, df, by="Feature")
str(importanceRR)
# let's identify the absolute value for each coefficient
summary(importanceRR$coeff)
importanceRR$coeff <- abs(importanceRR$coeff)
summary(importanceRR$coeff)
## best predictors
for (v in 0:1){
cat("\n\n")
cat("value==", v)
importanceRR<- importanceRR[order(importanceRR$coeff, decreasing=TRUE),]
print(head(importanceRR[importanceRR$sign==v,], n=10))
cat("\n")
cat(paste(unique(head(importanceRR$Features[importanceRR$sign==v], n=10)), collapse=", "))
}
# let's draw a graph with our results focusing on such coefficients
importanceRR$sign2 <- ifelse(importanceRR$sign>0,importanceRR$coeff ,importanceRR$coeff *-1)
# the twop 20 relevant words
importance2RR <- top_n(importanceRR, 20, sign2 )
importance2RR
# the twop 20 irrelevant words
importance3RR <- top_n(importanceRR, -20, sign2 )
importance3RR
importance_newRR <- rbind(importance2RR, importance3RR)
str(importance_newRR)
# reorder the features
importance_newRR <- mutate(importance_newRR, Feature= reorder(Feature, sign2 ))
importance_newRR$sign3<- ifelse(importance_newRR$sign>0,"Positive","Negative")
ggplot(importance_newRR, aes(Feature, sign2 , fill = sign3)) +
geom_col(show.legend = F) +
coord_flip() +
ylab("RR Coefficients (we recoded as negative the values for the negative reviews)") +
scale_fill_manual(values = c("orange", "blue")) +
labs(title = "Movie Reviews",
subtitle = "Negative (-) versus Positive (+) features - RR")
######################################################
######################################################
# Let's now move to a lasso regression model
######################################################
######################################################
# lasso regression requires that alpha=1
set.seed(123)
system.time(lasso <- cv.glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train,
family="binomial", alpha=1, nfolds=10, type.measure="class", trace.it=1 ))
plot(lasso )
min(lasso $cvm) # minimum miss-classification error of the lasso regression
min(ridge $cvm) # minimum miss-classification error of the ridge regression
lasso $lambda.min # lambda for this value
log(lasso $lambda.min)
lasso $lambda.1se
log(lasso $lambda.1se)
# Now the advantage of identifying the lambda with a miss-classification error within one standard error becomes more obvious.
# If we use the lambda that drives the minimum miss-classification error we can reduce our feature set to less than 109.
# However, there will be some variability with this miss-classification error and we can reasonably assume that we can achieve
# a similar miss-classification error with a slightly more constrained model that uses only 44 features.
# If describing and interpreting the predictors is an important outcome of your analysis, this may significantly aid your endeavor.
# This is not necessarily the case while doing text analytics with lots of features...still, that's up to you as a decision
# Note that whereas the ridge regression approach pushes variables to approximately but not equal to zero, the lasso penalty will actually
# push coefficients to zero
system.time(ROLS_lasso <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train, family="binomial", alpha=1, trace.it=1))
# as you can see as lambda increases, the penalty becomes large and forces our coefficients to zero
plot(ROLS_lasso, xvar = "lambda") # with lasso
plot(ROLS, xvar = "lambda") # with ridge
######################################################
######################################################
# Let's now move to elastic nets regression model
######################################################
######################################################
# Any alpha value between 0-1 will perform an elastic net.
# When alpha = 0.5 we perform an equal combination of penalties whereas alpha closer 0
# will have a heavier ridge penalty applied and alpha closer 1 will have a heavier lasso penalty.
# Ideally you could look for the optimal values of both lambda AND alpha. In the below example, however,
# we have fixed alpha = 0.5. Then cv.glmnet will not cross-validate across values of alpha
set.seed(123)
system.time(elastic<- cv.glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train,
family="binomial", alpha=0.5, nfolds=10, trace.it=1,
type.measure="class"))
plot(elastic)
min(elastic$cvm) # minimum miss-classification error of the elastic nets regression model
elastic$lambda.min # lambda for this value
log(elastic$lambda.min)
elastic$lambda.1se
log(elastic$lambda.1se)
# it seems that the ridge model is the one performing better in this case
min(ridge $cvm) # minimum miss-classification error of the ridge regression
min(lasso $cvm) # minimum miss-classification error of the lasso regression
min(elastic$cvm) # minimum miss-classification error of the elastic nets regression model
#### let's check MAE for the test-set (given that in this case we do know the true class-labels for each text in the test-set)
# RIDGE
system.time(ROLS <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train, family="binomial", alpha=0, trace.it=1))
system.time(predicted_glmentR_min <- predict(ROLS, test, s = ridge $lambda.min, type="class"))
# LASSO
system.time(ROLS <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train, family="binomial", alpha=1, trace.it=1))
system.time(predicted_glmentL_min <- predict(ROLS, test, s = lasso $lambda.min, type="class"))
# ELASTIC NETS
system.time(ROLS <- glmnet(y= as.factor(Dfm_train@docvars$Sentiment), x=train, family="binomial", alpha=0.5, trace.it=1))
system.time(predicted_glmentE_min <- predict(ROLS, test, s = elastic $lambda.min, type="class"))
mae<- as.data.frame(rbind(prop.table(table(predicted_glmentR_min)), prop.table(table(predicted_glmentL_min )),
prop.table(table(predicted_glmentE_min)), prop.table(table(x10$Sentiment))))
str(mae)
mae$algorithm <- c("Ridge", "Lasso", "Elastic", "TRUE")
mae_Ridge <- (mean(abs(mae[1,1]-mae[4,1]))+mean(abs(mae[1,2]-mae[4,2])))/2
mae_Lasso <- (mean(abs(mae[2,1]-mae[4,1]))+mean(abs(mae[2,2]-mae[4,2])))/2
mae_Elasitc <- (mean(abs(mae[3,1]-mae[4,1]))+mean(abs(mae[3,2]-mae[4,2])))/2
# Ridge better
mae_Ridge
mae_Lasso
mae_Elasitc
######################################################
######################################################
# Let's run a 5-folds CV with a Ridge Regression (to compare it with other ML algorithms for example)
######################################################
######################################################
# STEP 1: create the 5 folds
library(cvTools)
ttrain <- train # let's change the name of the original train data.frame, given that we are already going to use such name below in the loop
# let's split our training-set in 5 folds
set.seed(123) # set the see for replicability
k <- 5 # the number of folds; it does not matter the number of folds you decide here; the below procedure always will work!
folds <- cvFolds(NROW(ttrain ), K=k)
str(folds)
system.time(for(i in 1:k){
train <- ttrain [folds$subsets[folds$which != i], ] # Set the training set
validation <- ttrain [folds$subsets[folds$which == i], ] # Set the validation set
set.seed(123)
newrf <- cv.glmnet(y= as.factor(Dfm_train[folds$subsets[folds$which != i], ]@docvars$Sentiment), x=train,
family="binomial", alpha=0, nfolds=10, type.measure="class", trace.it=1 ) # Get your new model (just fit on the train data) and ADD the name of the output (in this case "Sentiment")
newpred <- predict(newrf,validation, s = newrf $lambda.min, type="class") # Get the predicitons for the validation set (from the model just fit on the train data)
class_table <- table("Predictions"= newpred, "Actual"=Dfm_train[folds$subsets[folds$which == i], ]@docvars$Sentiment)
print(class_table)
df<-confusionMatrix( class_table, mode = "everything")
df.name<-paste0("conf.mat.rr",i) # create the name for the object that will save the confusion matrix for each loop (=5)
assign(df.name,df)
})
# STEP 3: the metrics
ls()
# we have created 5 objects that have saved the 5 confusion matrices we have created. I can estimate now the performance metrics on such results
# for example:
conf.mat.rr1
RRPredict <- data.frame(col1=vector(), col2=vector(), col3=vector()) # makes a blank data frame with three columns to fill with the predictions;
# why 3 columns? 1 for accuracy; and 2 for the K1 value of the classes in the Sentiment (given you have just two classes!)
for(i in mget(ls(pattern = "conf.mat.rr")) ) {
Accuracy <-(i)$overall[1] # save in the matrix the accuracy value
F1_negative<- (2*(i)$byClass[1]*(i)$byClass[3])/((i)$byClass[1]+(i)$byClass[3]) # save in the matrix the F1 value for negative
F1_positive <- (2*(i)$byClass[2]*(i)$byClass[4])/((i)$byClass[2]+(i)$byClass[4]) # save in the matrix the F1 value for positive
RRPredict <- rbind(RRPredict , cbind(Accuracy , F1_negative, F1_positive ))
}
str(RRPredict )
RRPredict [is.na(RRPredict )] <- 0
# Let's compare the average value for accuracy and f1
acc_rr_avg <- mean(RRPredict[, 1] )
f1_rr_avg <- mean(colMeans(RRPredict[-1] ))
acc_rr_avg
f1_rr_avg