TubeSpam

The data comes from a freely available source here: http://www.dt.fee.unicamp.br/~tiago//youtubespamcollection/

Alberto, T. C., Lochter, J. V., & Almeida, T. A. 2015. TubeSpam: Comment Spam Filtering on YouTube. 2015 IEEE 14th International Conference on Machine Learning and Applications (ICMLA), 138–143. ieeexplore.ieee.org.

We want to build a model that looks at the linguistic properties of the data, and uses it for making predictions. class: 1= spam 0=ham

rm(list = ls())
gc()

##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells  608608 32.6    1392490 74.4   692780 37.0
## Vcells 1150298  8.8    8388608 64.0  1929374 14.8

library(devtools)

## Loading required package: usethis

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidytext)
library(quanteda)

## Package version: 4.2.0
## Unicode version: 13.0
## ICU version: 66.1
## Parallel computing: disabled
## See https://quanteda.io for tutorials and examples.

if (!require("quanteda.dictionaries", character.only = TRUE)) {
  if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
  remotes::install_github("kbenoit/quanteda.dictionaries")
}

## Loading required package: quanteda.dictionaries

library(quanteda.dictionaries)

Import the data

Read in and compile each of the individual files that start with ‘Youtube’. Using lapply to bind them all together. Make column headers lower case. Remove duplicates and missing values. Inspect the data.

Create tokens.

Regular expressions cheatsheets: https://hypebright.nl/index.php/en/2020/05/25/ultimate-cheatsheet-for-regex-in-r-2/ and https://gist.github.com/shakahl/43443dcba2c8105218350bf14df2c610 are both nice.

The command \w matches word characters, things that are used in writing typical English words. unnest_tokens will create tokens (unigrams). New column is called ‘word’, and will contain content from ‘content’

Text cleaning step!

data("stop_words")

yt %>%  unnest_tokens(word, content) %>% #give new column a name - we called it 'word'
  count(word, sort=TRUE)

## # A tibble: 4,138 × 2
##    word      n
##    <chr> <int>
##  1 i       636
##  2 this    606
##  3 the     570
##  4 to      498
##  5 a       458
##  6 and     451
##  7 my      412
##  8 you     406
##  9 out     367
## 10 check   364
## # ℹ 4,128 more rows

yt_tokens <- yt %>% 
  unnest_tokens(word, content) %>%      # create the tokens
  mutate(word = tolower(word)) %>%      # make sure everything is lowercase
  anti_join(stop_words) %>%             # remove the stop words
  filter(nchar(word) >= 2) %>%          # remove anything that's short
  filter(grepl("\\w+", word))           # delete any punctuation or symbols (words only)

## Joining with `by = join_by(word)`

x <- paste(sort(unique(yt_tokens$word)), collapse = ' ')  # unique tokens only sorted - just to have a look

Most common words?

yt_tokens %>% 
  count(word) %>% 
  arrange(-n)

## # A tibble: 3,687 × 2
##    word          n
##    <chr>     <int>
##  1 check       364
##  2 br          250
##  3 video       238
##  4 song        233
##  5 youtube     198
##  6 love        172
##  7 subscribe   161
##  8 39          156
##  9 http        143
## 10 amp         129
## # ℹ 3,677 more rows

yt_tokens %>%  filter(word == "39")

## # A tibble: 156 × 6
##    class file                             id    author comment_date        word 
##    <dbl> <chr>                            <chr> <chr>  <dttm>              <chr>
##  1     0 Youtube 01-comments Psy.csv      z135… Spenc… 2014-11-08 05:29:26 39   
##  2     0 Youtube 04-comments KatyPerry.c… z13m… micha… 2014-10-30 05:20:18 39   
##  3     1 Youtube 07-comments LMFAO.csv    z13t… LaS M… 2015-05-28 19:23:35 39   
##  4     1 Youtube 07-comments LMFAO.csv    z13t… LaS M… 2015-05-28 19:23:35 39   
##  5     0 Youtube 07-comments LMFAO.csv    z13y… LoL G… 2015-05-27 13:37:19 39   
##  6     0 Youtube 07-comments LMFAO.csv    z12s… breat… 2015-05-26 22:21:28 39   
##  7     0 Youtube 07-comments LMFAO.csv    z12l… goofy… 2015-05-26 20:03:52 39   
##  8     0 Youtube 07-comments LMFAO.csv    z12s… Like👑… 2015-05-26 16:56:19 39   
##  9     0 Youtube 07-comments LMFAO.csv    z120… Shock… 2015-05-26 02:03:46 39   
## 10     0 Youtube 07-comments LMFAO.csv    z120… Shock… 2015-05-26 02:03:46 39   
## # ℹ 146 more rows

yt %>% filter(id == "z12sulyo4qmyjvjvm23tznvzfnboetfuj") %>% pull(content)

## [1] "I fuckin love this song!<br /><br /><br />After, I&#39;m sexy and I know it "

Some words are used frequently in every comment, we want to know terms that are used in specific comments.

Compute the TF-IDF for the terms.

word_tfidf <- yt_tokens %>% 
  count(id, word, sort = T) %>% 
  bind_tf_idf(word, id, n)

word_tfidf %>% 
  arrange(-tf_idf)

## # A tibble: 10,461 × 6
##    id                                          word         n    tf   idf tf_idf
##    <chr>                                       <chr>    <int> <dbl> <dbl>  <dbl>
##  1 _2viQ_Qnc68fX3dYsfYuM-m4ELMJvxOQBmBOFHqGOk0 beutiful     1     1  7.43   7.43
##  2 z125v3ozoqenvthaz04cdtajsmzwgxkwxug0k       superr       1     1  7.43   7.43
##  3 z12agrbxiwabvrb2g22pfpcqys2eyrwbd04         nicei        1     1  7.43   7.43
##  4 z12ai5agdlawczj5x04cfhbr2vuezzvox1s         collabo…     1     1  7.43   7.43
##  5 z12bjffg3yu4ybx4e04cjzfbyxqcdzyxl3s         hardcore     1     1  7.43   7.43
##  6 z12gcfvpsr3au3vg104cdpw5rvbivre5aeg         eminen       1     1  7.43   7.43
##  7 z12jf1h5lznmuhcuw22ruxdyxp2eh32ch           likkee       1     1  7.43   7.43
##  8 z12jjrfyukjlw5wyl04cddbo5pysslphbqs0k       epic         1     1  7.43   7.43
##  9 z12kstz4owbggv12p230vpcwntjbyrqim04         fit          1     1  7.43   7.43
## 10 z12lcb0rxw3qtlsqm04cebpbjkqjtdza1xs0k       telepho…     1     1  7.43   7.43
## # ℹ 10,451 more rows

Number of unique words.

n_distinct(yt_tokens$word)

## [1] 3687

count(yt_tokens, word, sort = T) %>% head(10) # common words

## # A tibble: 10 × 2
##    word          n
##    <chr>     <int>
##  1 check       364
##  2 br          250
##  3 video       238
##  4 song        233
##  5 youtube     198
##  6 love        172
##  7 subscribe   161
##  8 39          156
##  9 http        143
## 10 amp         129

count(yt_tokens, word, sort = T) %>% tail(10) # least common words

## # A tibble: 10 × 2
##    word                              n
##    <chr>                         <int>
##  1 ｃｏｍｍｅｎｔ                    1
##  2 ｄａｍｎ                          1
##  3 ｅｂａｙ                          1
##  4 ｆａｎｃy                         1
##  5 ｈｔｔｐ                          1
##  6 ｉｓ                              1
##  7 ｓｈｏｅｃｏｌｌｅｃｔｏｒ314     1
##  8 ｔｈｉｓ                          1
##  9 ｕｓｒ                            1
## 10 ｗｗｗ                            1

Create test and training sets the same way that they did in the paper. They did it by using the first 70% of each set of comments as the training, and tested on the last 30% of the comments. We can’t do it across all comments since - well, look:

ggplot(yt, aes(x = comment_date, y = file)) +
  geom_jitter(alpha = 0.2) +      # try comparing with geom_point and/or removing alpha
  theme_minimal()                 # try removing this line (and the +)

The comment dates are very different depending on the video. Look at training and test sets - test set is most recent comments

yt <- yt %>% 
  group_by(file) %>% 
  mutate(date_decile = ntile(comment_date, 10)) %>%              # break up and divide into 10 even parts based on date
  mutate(train = ifelse(date_decile <= 7, 'train', 'test')) %>%  # train-test split on each file
  ungroup()
ggplot(yt, aes(x = comment_date, y = file, color = train)) +     # visualise the split!
  geom_jitter(alpha = 0.2) +
  theme_minimal()

Need to cast it into a matrix for these next several methods. This is a document term matrix - we talked about them in the lecture - a frequency table of words for each observation.

Its very long - we will make it very wide

https://www.tidytextmining.com/dtm.html

yt_m <- word_tfidf %>% 
  cast_dtm(id, word, tf_idf) %>%      # rows, columns, information inside it
  as.matrix

The DTM massive! Don’t try to look at it yet. Lets’ look at the top left corner of this tf-idf term document matrix looks like. (You could calculate a distance matrix from this, and then do some clustering or pca on it…but we won’t)

yt_m[1:5, 1:5]

##                                              Terms
## Docs                                                 amp    48051       http
##   z12jenlhyre0eheyx04ch1aquxfdsvgpd44         0.55340759 0.000000 0.06238548
##   z131idupvn3yhf3mv23dwzhi4pqixvwuw           0.00000000 1.456885 0.54434388
##   z132cvvy1ob3ht2er23dundqdtertjmlg           0.05584272 0.000000 0.04406593
##   _2viQ_Qnc69MEEHHJxZ427KX8MlljJPnUC2YBbvbWwY 0.00000000 0.000000 0.00000000
##   _2viQ_Qnc6_RKHVetk9kLzx8ZC62_J7y73FWFSBTe8Q 0.00000000 0.000000 0.00000000
##                                              Terms
## Docs                                          image2you       ru
##   z12jenlhyre0eheyx04ch1aquxfdsvgpd44          0.000000 0.000000
##   z131idupvn3yhf3mv23dwzhi4pqixvwuw            1.456885 1.456885
##   z132cvvy1ob3ht2er23dundqdtertjmlg            0.000000 0.000000
##   _2viQ_Qnc69MEEHHJxZ427KX8MlljJPnUC2YBbvbWwY  0.000000 0.000000
##   _2viQ_Qnc6_RKHVetk9kLzx8ZC62_J7y73FWFSBTe8Q  0.000000 0.000000

Tokens as Predictors

Create the test and train datasets

We had already split the dataset by marking rows as “train” or “test.”

Then subset the main matrix (yt_m) into training (xs_train) and testing (xs_test) based on the IDs.

Finally, match the labels (class column) to the corresponding IDs, creating factor variables for ‘ham’ vs. ‘spam’.

This gives a standard X_train, y_train, X_test, y_test setup for modeling and evaluation.

id_train <- yt %>% filter(train == 'train') %>% pull(id)   # grab id's if train
id_test <- yt %>% filter(train == 'test') %>% pull(id)     # grab id's for test

xs_train <- yt_m[which(rownames(yt_m) %in% id_train),]     # filter matrix for training 
xs_test <- yt_m[which(rownames(yt_m) %in% id_test),]       # filter matrix for testing 

ys_train <- yt %>% filter(id %in% rownames(xs_train)) %>% 
  pull(class) %>%                                          # filter class for training id (labels)
  factor(levels = c(0, 1), labels = c('ham', 'spam'))
ys_test<- yt %>% filter(id %in% rownames(xs_test)) %>% 
  pull(class) %>%                                          # filter class for testing id (labels)
  factor(levels = c(0, 1), labels = c('ham', 'spam'))

NOTE!!! If you want to run this, do it later. xs_train (our matrix for prediction) is very large. It will take ages to run as it has a lot of variables (over 3700) - you can check with dim(xs_train) You can just load the model that has been already fit.

library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(randomForest)

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:ggplot2':
## 
##     margin

# rf1 <- randomForest(x = xs_train, y = ys_train, do.trace = T,
#                     ntree = 90)
# save.image(file = 'tubespam_rf1.rda')
load('./tubespam_rf1.rda')

OMG, this model is horribly overfit. How do we know?

#train
pred <- tibble(spam_p = predict(rf1, newdata = xs_train, type = 'prob')[,'spam'],
               spam = ys_train)
yt_roc_train <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = ham, case = spam

## Setting direction: controls < cases

#test
pred <- tibble(spam_p = predict(rf1, newdata = xs_test, type = 'prob')[,'spam'],
               spam = ys_test)
yt_roc_test <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = ham, case = spam
## Setting direction: controls < cases

# using base plotting
plot(yt_roc_train, print.auc = T)
lines(yt_roc_test, lty = 2, col = 'blue', print.auc = T)
text(round(auc(yt_roc_test), 3), x = 0.5, y = 0.4, col = "blue", pos = 4)

Might be worth reducing the mtry (number of variables to use)

In randomForest package, mtry controls how many features (predictor variables) are randomly sampled at each split in a decision tree.

If your dataset has a very large number of features, lowering mtry can sometimes improve performance by reducing overfitting and promoting greater diversity across individual trees. However, setting it too low may hamper the forest’s ability to find strong predictors.

A feature importance plot helps you see which variables are driving the classification decisions. Importance can be measured via metrics like like the Mean Decrease in Gini (how much splitting on a variable improves node purity).

If you notice that only a few variables dominate the importance chart, or your model takes very long to train, you might try lowering mtry.

This can:

• Prevent a handful of strong features from always being chosen.

• Encourage your model to consider a broader range of features at each split, potentially improving generalisation.

varImpPlot(rf1)

Linguistic Features

Create generalized linguistic features. First I have to load up a linguistic dictionary called NRC Emotion Lexicon. It is available here: https://saifmohammad.com/WebPages/NRC-Emotion-Lexicon.htm And is free for use for research or education, with acknowledgement. Thank you Dr. Saif M. Mohammad at the National Research Council Canada.

I have to read it in as a text file, then convert it to a different format that we can use. You can open up the text file and check out what Dr Mohammad has done to classify words as emotions.

# Split each line into components and convert to a data frame
nrc <- readLines("NRC-Emotion-Lexicon-Wordlevel-v0.92.txt")

data <- do.call(rbind, strsplit(nrc, "\t"))
df <- as.data.frame(data, stringsAsFactors = FALSE)
colnames(df) <- c("Word", "Emotion", "Association")

# Convert Association column to numeric
df$Association <- as.numeric(df$Association)

# Filter out rows where Association is 0, as we're only interested in associations marked as 1
df <- df[df$Association == 1, ]

# Create an empty list to store emotions as keys and associated words as values
emotion_words_dict <- list()

# Populate the list
for (row in 1:nrow(df)) {
  emotion <- df$Emotion[row]
  word <- df$Word[row]
  
  # Check if the emotion already exists as a key in the list
  if (!is.null(emotion_words_dict[[emotion]])) {
    # Append the word to the existing list of words for this emotion
    emotion_words_dict[[emotion]] <- c(emotion_words_dict[[emotion]], word)
  } else {
    # Create a new entry in the list with this emotion as the key and the word as the value
    emotion_words_dict[[emotion]] <- c(word)
  }
}

emotion_dict <- dictionary(emotion_words_dict)

# Check the structure of the created dictionary
print(emotion_dict)

## Dictionary object with 10 key entries.
## - [trust]:
##   - abacus, abbot, absolution, abundance, academic, accolade, accompaniment, accord, account, accountability, accountable, accountant, accounts, accredited, accurate, achieve, achievement, acrobat, adhering, administrative [ ... and 1,210 more ]
## - [fear]:
##   - abandon, abandoned, abandonment, abduction, abhor, abhorrent, abominable, abomination, abortion, absence, abuse, abyss, accident, accidental, accursed, accused, accuser, accusing, acrobat, adder [ ... and 1,454 more ]
## - [negative]:
##   - abandon, abandoned, abandonment, abduction, aberrant, aberration, abhor, abhorrent, abject, abnormal, abolish, abolition, abominable, abomination, abort, abortion, abortive, abrasion, abrogate, abscess [ ... and 3,296 more ]
## - [sadness]:
##   - abandon, abandoned, abandonment, abduction, abortion, abortive, abscess, absence, absent, absentee, abuse, abysmal, abyss, accident, accursed, ache, aching, adder, adrift, adultery [ ... and 1,167 more ]
## - [anger]:
##   - abandoned, abandonment, abhor, abhorrent, abolish, abomination, abuse, accursed, accusation, accused, accuser, accusing, actionable, adder, adversary, adverse, adversity, advocacy, affront, aftermath [ ... and 1,225 more ]
## - [surprise]:
##   - abandonment, abduction, abrupt, accident, accidental, accidentally, accolade, advance, affront, aghast, alarm, alarming, alertness, alerts, allure, amaze, amazingly, ambush, angel, anomaly [ ... and 512 more ]
## [ reached max_nkey ... 4 more keys ]

I have now used it by applying it to our raw content.

dic_features <- yt %>%
  with(liwcalike(corpus(content, docnames = id),
                 dictionary = emotion_dict,tolower = T)) %>% 
  rename(id = docname)

We want to find if there are URL’s in the comment. Google provided an answer:

https://stackoverflow.com/questions/3809401/what-is-a-good-regular-expression-to-match-a-url

# website_regex <- "(www\\.)?([a-zA-Z0-9]+(-?[a-zA-Z0-9])*\\.)+[\\w]{2,}(\\/\\S*)?"
website_regex <- "[-a-zA-Z0-9@:%._\\+~#=]{1,}\\.[a-zA-Z0-9()]{1,6}\\b([-a-zA-Z0-9()@:%_\\+.~#?&//=]*)"
yt_url <- yt %>% 
  ungroup() %>% 
  mutate(content = tolower(content)) %>% 
  mutate(url = grepl(website_regex, content)) %>%      #transform to lowercase and use grepl pattern match
  mutate(url = ifelse(url, 1, 0)) %>% 
  select(id, url)

Now use left_join to join our url findings to our dic_features (joining by id) re-join with original youtube data so can get original info

yt_ling_features <- left_join(dic_features, yt_url) %>% 
  left_join(select(yt, id, class, train, content)) %>% 
  mutate(class = as.factor(class))

## Joining with `by = join_by(id)`
## Joining with `by = join_by(id)`

glimpse(yt_ling_features)

## Rows: 1,710
## Columns: 32
## $ id           <chr> "LZQPQhLyRh80UYxNuaDWhIGQYNQ96IuCg-AYWqNPjpU", "LZQPQhLyR…
## $ Segment      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ WPS          <dbl> 13.00, 19.50, 8.00, 12.00, 4.50, 9.50, 4.00, 21.00, 10.00…
## $ WC           <int> 13, 39, 8, 12, 9, 19, 4, 21, 10, 15, 4, 37, 7, 3, 15, 71,…
## $ Sixltr       <dbl> 15.38, 10.26, 12.50, 16.67, 11.11, 10.53, 50.00, 0.00, 10…
## $ Dic          <dbl> 0.00, 10.26, 0.00, 58.33, 22.22, 0.00, 25.00, 0.00, 0.00,…
## $ trust        <dbl> 0.00, 0.00, 0.00, 8.33, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ fear         <dbl> 0.00, 0.00, 0.00, 8.33, 11.11, 0.00, 0.00, 0.00, 0.00, 0.…
## $ negative     <dbl> 0.00, 2.56, 0.00, 16.67, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ sadness      <dbl> 0.00, 2.56, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ anger        <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ surprise     <dbl> 0.00, 2.56, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ positive     <dbl> 0.00, 0.00, 0.00, 8.33, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ disgust      <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ joy          <dbl> 0.00, 0.00, 0.00, 8.33, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ anticipation <dbl> 0.00, 2.56, 0.00, 8.33, 11.11, 0.00, 25.00, 0.00, 0.00, 0…
## $ AllPunc      <dbl> 30.77, 23.08, 12.50, 8.33, 22.22, 26.32, 0.00, 14.29, 20.…
## $ Period       <dbl> 0.00, 0.00, 12.50, 0.00, 11.11, 10.53, 0.00, 14.29, 0.00,…
## $ Comma        <dbl> 7.69, 2.56, 0.00, 0.00, 0.00, 5.26, 0.00, 0.00, 0.00, 0.0…
## $ Colon        <dbl> 7.69, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ SemiC        <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ QMark        <dbl> 0.00, 0.00, 0.00, 0.00, 11.11, 0.00, 0.00, 0.00, 0.00, 0.…
## $ Exclam       <dbl> 0.00, 17.95, 0.00, 0.00, 0.00, 10.53, 0.00, 0.00, 20.00, …
## $ Dash         <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ Quote        <dbl> 0.00, 2.56, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ Apostro      <dbl> 0.00, 2.56, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.0…
## $ Parenth      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ OtherP       <dbl> 15.38, 23.08, 12.50, 0.00, 22.22, 26.32, 0.00, 14.29, 20.…
## $ url          <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, …
## $ class        <fct> 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, …
## $ train        <chr> "train", "train", "train", "train", "train", "train", "tr…
## $ content      <chr> "Huh, anyway check out this you[tube] channel: kobyoshi02…

filter into test and train and remove features we don’t want

yt_ling_train <- yt_ling_features %>% 
  filter(train == 'train') %>% 
  select(-id, -train, -content, -Segment)         #remove variables we don't want
yt_ling_test <- yt_ling_features %>% 
  filter(train == 'test') %>% 
  select(-id, -train, -content, -Segment)

save(yt_ling_features, yt_ling_train, yt_ling_test, file = 'youtube_spam_train_test.rda')

Linguistic Features as Predictors

rf2 <- randomForest(class ~ ., 
                    data = yt_ling_train,
                    ntree = 150)

It’s still overfit, but its a lot better for the test set. Now, instead of using specific words to predict, we are using generalisable emotions to predict

#train
pred <- tibble(spam_p = predict(rf2, newdata = yt_ling_train, type = 'prob')[,'1'],
               spam = yt_ling_train$class)
yt_roc_train <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

#test
pred <- tibble(spam_p = predict(rf2, newdata = yt_ling_test, type = 'prob')[,'1'],
               spam = yt_ling_test$class)
yt_roc_test <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

# using base plotting
plot(yt_roc_train, print.auc = T)
lines(yt_roc_test, lty = 2, col = 'blue', print.auc = T)
text(round(auc(yt_roc_test), 3), x = 0.5, y = 0.4, col = "blue", pos = 4)

We don’t know the causal direction, so like your assignment, we will create some fake data.

varImpPlot(rf2)

Explaining the RF2 model

First using fake data with 1000 values!

#set.seed(1223)
N <- 1000
g <- seq(min(yt_ling_train$Sixltr),
         max(yt_ling_train$Sixltr),
         length.out = N)
#d <- sample_n(yt_ling_train, 1)
d <- yt_ling_features[38,]

#quick check of a prediction (its a true positive)
predict(rf2, newdata = d)

## 38 
##  1 
## Levels: 0 1

dn <- d %>%
  mutate(n = N) %>%
  uncount(n) %>%
  mutate(Sixltr = g)
dn$pred <- predict(rf2, newdata = dn, type = "prob")[,2]

#Lets play with the most important variable a bit before we look

predict(rf2, newdata =d, type="prob")

##      0   1
## 38 0.1 0.9
## attr(,"class")
## [1] "matrix" "array"  "votes"

d_change <- d %>% mutate(Sixltr = 50)
predict(rf2, newdata =d_change, type="prob")

##      0   1
## 38 0.1 0.9
## attr(,"class")
## [1] "matrix" "array"  "votes"

d_change2 <- d %>% mutate(Sixltr = 2)
predict(rf2, newdata =d_change2, type="prob")

##            0         1
## 38 0.3733333 0.6266667
## attr(,"class")
## [1] "matrix" "array"  "votes"

ggplot(dn, aes(x = Sixltr, y = pred)) +
  geom_line() +
  coord_cartesian(xlim = c(0, 100)) +
  xlab("Percent of words that are six letters or longer") +
  ylab("P(spam)")

library(iml)
d <- sample_n(yt_ling_train, 200)
yt_pred <- Predictor$new(model = rf2, data = d,
                         y = d$class, type = "prob")
yt_effect <- FeatureEffect$new(yt_pred,center.at = 0.5,
                               grid.size = 100,
                               method = "pdp",
                               feature = "Sixltr")
yt_effect$plot() +
  xlab("Percent of words that are six letters or longer") +
  ylab("P(spam)")

—- We can stop here for today! —–

Topic model

Prep the term document matrix for counts instead of TD-IDF We don’t need to clean, will use all of the words

library(topicmodels) 
library(LDAvis)
library(gt)

yt_dtm <- yt_tokens %>% 
  count(id, word) %>% 
  cast_dtm(id, word, n)

Choose the number of topics and run the LDA.

k <- 10
seed <- 2015
yt_lda <- LDA(yt_dtm, k=k, method = 'Gibbs', control=list(seed=seed, burnin=1000, thin=100, iter=1000))

Top 10 terms from each topic.

terms(yt_lda, 10) %>% 
  as_tibble() %>% 
  gt

Topic 1	Topic 2	Topic 3	Topic 4	Topic 5	Topic 6	Topic 7	Topic 8	Topic 9	Topic 10
br	guys	https	amp	check	video	views	39	subscribe	song
quot	free	im	lt	http	youtube	billion	music	channel	love
youtube	money	www.facebook.com	nice	mixtape	check	shakira	comment	watch	katy
share	visit	watching	gt	awesome	playlist	people	time	videos	perry
waka	day	ref	hear	href	channel	pray	world	follow	songs
fucking	facebook	cool	3	48051	it	song	thumbs	shit	girl
click	website	2015	style	image2you.ru	party	wow	chance	hate	beautiful
rapper	paid	guy	gangnam	vote	megan	subscribers	money	plz	omg
birthday	twitter	class	http	support	hey	mother	covers	funny	remember
miley	online	dance	songs	rand	rock	d	hey	hey	roar

LDA visualization tool.

lda_json <- createJSON(phi = posterior(yt_lda)$terms,
           theta = posterior(yt_lda)$topics,
           vocab = posterior(yt_lda)$terms %>% colnames,
           doc.length = rowSums(as.matrix(yt_dtm)),
           term.frequency = colSums(as.matrix(yt_dtm)))
serVis(lda_json)

## Loading required namespace: servr

Predicting with topics.

You end up with a single dataset (yt_topics) that contains:

All original columns from yt (id, class, etc.) Topic proportion columns (Topic_1, Topic_2, …, Topic_k) for each document. This allows you to use the topics as features for further analysis, such as classification, clustering, or just exploring how topics vary with the class label.

posterior(yt_lda)$topics[1:5, 1:5]

##                                                      1          2          3
## LZQPQhLyRh80UYxNuaDWhIGQYNQ96IuCg-AYWqNPjpU 0.09090909 0.09090909 0.10909091
## LZQPQhLyRh9-wNRtlZDM90f1k0BrdVdJyN_YsaSwfxc 0.08771930 0.14035088 0.10526316
## LZQPQhLyRh9EXArr4ZnVcDonSbvSMHKYOT24e_qR6fE 0.09615385 0.09615385 0.09615385
## LZQPQhLyRh9MSZYnf8djyk0gEF9BHDPYrrK-qCczIY8 0.11538462 0.09615385 0.09615385
## LZQPQhLyRh9U7Lv_DKpJ7lawpBCotxfgHzBy93Tk028 0.09615385 0.09615385 0.09615385
##                                                      4          5
## LZQPQhLyRh80UYxNuaDWhIGQYNQ96IuCg-AYWqNPjpU 0.10909091 0.09090909
## LZQPQhLyRh9-wNRtlZDM90f1k0BrdVdJyN_YsaSwfxc 0.08771930 0.08771930
## LZQPQhLyRh9EXArr4ZnVcDonSbvSMHKYOT24e_qR6fE 0.09615385 0.11538462
## LZQPQhLyRh9MSZYnf8djyk0gEF9BHDPYrrK-qCczIY8 0.09615385 0.09615385
## LZQPQhLyRh9U7Lv_DKpJ7lawpBCotxfgHzBy93Tk028 0.09615385 0.09615385

yt_topics <- posterior(yt_lda)$topics
colnames(yt_topics) <- paste0("Topic_", colnames(yt_topics))

yt_topics <- as_tibble(yt_topics) %>% 
  mutate(id = rownames(yt_topics))

yt_topics <- left_join(yt, yt_topics) %>% 
  mutate(class = as_factor(class)) %>% 
  na.omit()

## Joining with `by = join_by(id)`

Comments are weighted to every topic - we can use those to make predictions. Create our test/train sets

yt_topic_train <- yt_topics %>% 
  filter(train == "train") %>%    
  select(class, Topic_1:Topic_10)
yt_topic_test <- yt_topics %>% 
  filter(train == "test") %>% 
  select(class, Topic_1:Topic_10)

rf3 <- randomForest(class ~ ., data = yt_topic_train,
                    ntree = 300)

Not as good as the linguistic features, but better than raw words, but you could use both!

#train
pred <- tibble(spam_p = predict(rf3, newdata = yt_topic_train, type = 'prob')[,'1'],
               spam = yt_topic_train$class)
yt_roc_train <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

#test
pred <- tibble(spam_p = predict(rf3, newdata = yt_topic_test, type = 'prob')[,'1'],
               spam = yt_topic_test$class)
yt_roc_test <- roc(data = pred, response = spam, predictor = spam_p)

## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

# using base plotting
plot(yt_roc_train, print.auc = T)
lines(yt_roc_test, lty = 2, col = 'blue', print.auc = T)
text(round(auc(yt_roc_test), 3), x = 0.5, y = 0.4, col = "blue", pos = 4)

Check most important variables

varImpPlot(rf3)

What does Topic 10 do?

d <- sample_n(yt_topic_train, 200)
yt_pred <- Predictor$new(model = rf3, data=d, y=d$class, type="prob")
yt_effect <- FeatureEffect$new(yt_pred, center.at = 0.5,
                               grid.size = 100,
                               method = "pdp+ice",
                               feature = "Topic_10")

yt_effect$plot() +
  #coord_cartesian(xlim = c(0,100))+
  xlab("Topic 10")+
  ylab("P(spam)")+
  scale_x_log10()