Lab 5 — tree-based methods
Info about the lab
Learning aim
The aim of this lab is to apply tree-based models for natural language processing and to learn how to train and tune tree-based models.
Objectives
By the end of this lab session, students should be able to
Visualize text data.
Extract features from a text data.
Train and plot decision trees and regularized decision trees.
Train random forest and tune the hyperparameters in random forest to minimize the out-of-bag error.
Preliminaries
You should be familiar with libraries caret and
tidyverse.
Mode
Please run the R chunks one by one, look at the output and make sure that you understand how it is produced. There will be questions that either require a short answer - then you type your answer right in this document - or modifying R codes - then you modify the R codes here. In either case, you can discuss your work with the lab instructor.
Data
We will work with the data of Donald Trump’s tweets. Source:
The objective is to predict the number of retweets from the text of each tweet.
Loading data into R
First we will load the data into R. Note that the original dataset is a JSON file that contains a lot of information. For this lab, we only extracted the text and the number of retweets from 5K out of original 40K Trump’s tweets.
library(tidyverse) # for manipulation with data
library(tidytext) # for tidyverse-style text tokenization
library(caret) # for machine learning, including KNN
library(rpart) # for training decision trees
library(rpart.plot) # for plotting decision trees## Warning: package 'rpart.plot' was built under R version 4.6.1library(ranger) # for training random forest
# Note for students: the libraries `cforest` and `randomForest`
# provide better access to guts of random forest models.
# That's why these libraries were used to prepare lecture 5 slides.
# However, cforest and randomForest do not allow tuning the minimum node size.
library(ggwordcloud) # for text visualization
set.seed(8192)
tweets <- read_csv("trump_tweets.csv")
head(tweets)The number of tweets is
nrow(tweets)## [1] 5000Data cleaning
We will remove all non-character symbols except for “@” and “#” because they are special symbols. Then we will change all characters to the lower case. Below is the processed dataset.
tweets <- tweets %>%
mutate(
cleaned_text = text %>%
str_replace_all("[^A-Za-z@#]", " ") %>%
str_squish() %>%
str_to_lower()
)
head(tweets)Question 1
The symbols “@” and “#” have special meaning in twitter and it is probably a good idea to treat them as separate words. Write an R code that replaces all occurrences of “@” with “twacc” and all occurrences of “#” with “htag” in Trump’s tweets.
# Modify the code below
tweets <- tweets %>%
mutate(
cleaned_text = cleaned_text %>%
str_replace_all("@", "twacc ") %>%
str_replace_all("#", "htag ")
)
head(tweets)Document-term matrix
An important tool in natural language processing is document-term-matrix. This is a matrix whose rows represent documents, i.e., tweets, and whose columns represent words. The element \(T_{ij}\) in row \(i\) and column \(j\) is the number of occurrences of word \(j\) in document \(i\). To construct the document-term matrix, the first step is to create a data frame of word counts in each document:
tidy_counts <- tweets %>%
mutate(doc_id = row_number()) %>%
select(doc_id, cleaned_text) %>%
unnest_tokens(word, cleaned_text) %>%
count(doc_id, word, sort = FALSE)
tidy_counts %>% sample_n(10)Question 2
How many different words are there in the vocabulary?
Given dimensions of the term-document matrix, is it possible to use linear regression to predict the number of retweets from the term-document matrix?
Answer The number of words in the vocabulary is
n_distinct(tidy_counts$word)## [1] 12674It is impossible to use linear regression because the number of predictors is greater than the number of variables.
If we included all words that are mentioned at least once in any tweet, we would have too many words and most of them wouldn’t be useful since the appear only in one or two tweets.
To reduce the number of independent variables, we will only keep words that appear in the whole corpus (collection of all tweets), say, at least 30 times. We won’t keep stopwords either since they aren’t meaningful for our purpose.
keep_words <- tidy_counts %>%
count(word, wt = n, name = "freq") %>%
filter(freq >= 30, !word %in% stop_words$word)
keep_words %>% sample_n(10)The next step is to subset tidy_counts to include only
frequent words. There are two ways to do it, the first one is more
efficient (works faster for large datasets), but second one is easier to
read.
# tidy_counts_f <- tidy_counts %>%
# semi_join(keep_words, by = "word")
tidy_counts_f <- tidy_counts %>%
filter(word %in% keep_words$word)Finally, we can create the document-term matrix. Note that it has a
special class, so we will immediately create a data frame out of it. You
are encourated to explore DTM itself to get a feeling what
it looks like.
Again, there are two ways to create the document-term matrix. The first one (commented out below) is more efficient if the number of predictors is large but takes some effort to understand:
# all_docs <- seq_len(nrow(tweets))
# DTM <- tidy_counts_f %>%
# mutate(doc_id = factor(doc_id, levels = all_docs),
# word = factor(word, levels = keep_words$word)) %>%
# cast_sparse(doc_id, word, n)
#
# all_data <- DTM %>%
# as.matrix() %>%
# as_tibble(.name_repair = "minimal") %>%
# mutate(Y = tweets$Y)
all_data <- tibble(doc_id = 1:nrow(tweets)) %>%
left_join(tidy_counts_f) %>%
pivot_wider(id_cols = doc_id,
names_from = "word",
values_from = "n",
values_fill = 0) %>%
select(-any_of(c("doc_id", "NA"))) ## Joining with `by = join_by(doc_id)`all_data %>%
select(amp, iran, golf, hillary, twacc) %>%
slice(1:5)Data visualization
Text data are often visualized with a word cloud. Below is the word cloud of Trump’s tweets.
keep_words %>%
slice_max(freq, n = 100) %>%
ggplot(aes(label = word, size = freq)) +
geom_text_wordcloud() +
scale_size_area(max_size = 15) +
theme_void()
# wordcloud(words = keep_words$word, freq = keep_words$freq, min.freq = 0,
# max.words = 100, random.order=FALSE, rot.per=0.35)Question 3
Plot the word cloud without the tokens “twacc”, “realdonaldtrump”, “htag”, “http”, and “https”.
df_to_word_cloud <- function(word_data) {
word_data %>%
slice_max(freq, n = 100) %>%
ggplot(aes(label = word, size = freq)) +
geom_text_wordcloud() +
scale_size_area(max_size = 15) +
theme_void()
}
keep_words %>%
filter(!word %in% c("twacc", "realdonaldtrump", "htag", "http","https")) %>%
df_to_word_cloud
Question 4
Let us plot the distribution of the number of retweets. As shown in the plot below, it is highly skewed:
ggplot(data = tweets, aes(x = retweet_count)) +
geom_histogram(fill = "orange", color = "black", bins = 20)
It makes more sense to predict the logarithm of the number of retweets than the raw number of retweets then. Compute the logarithm of the number of retweets and store it in the variable “Y” of our dataset. Warning: think of what you will do with tweets with 0 retweet count.
Plot the histogram of the log-retweet_count
# Modify the code below
# You should see a nice bimodal distribution
tweets <- tweets %>%
mutate(Y = log(1 + retweet_count))
all_data <- all_data %>% mutate(Y = tweets$Y)
ggplot(data = tweets, aes(x = Y)) +
geom_histogram(fill = "orange", color = "black", bins = 20) +
xlab("Log-retweet-count")
You can try to figure out why the distribution is bimodal, but it is not necessary.
Training and test sets
We will randomly split the data into 80% training and 20% test sets.
p <- 0.8
ind <- runif(nrow(all_data)) < p
train_data <- all_data %>% filter(ind)
test_data <- all_data %>% filter(!ind)
cat("Dimensions of the training set are", dim(train_data),"\n")## Dimensions of the training set are 3983 178cat("Dimensions of the test set are", dim(test_data),"\n")## Dimensions of the test set are 1017 178Decision tree
We will use the raw interface of the library rpart here
rather than wrapping it to train. This is because it is a
bit simpler and rpart automatically does cross-validation
for choosing the optimal value of the regularization constant.
Training and plotting
Remark As a rule of thumb, it is preferable to use
tidyverse tools such as ggplot() for plotting,
because they are much easier to customize and extend than base R
graphics. In practice, if you need anything more elaborate than a quick
plot(), ggplot() will usually be the better
choice. An exception are simple diagnostic plots that come bundled with
modeling packages (e.g., decision trees with rpart.plot,
cross-validation error curves with plotcp() etc). These are
often implemented in base R and are convenient to use as-is.
First, we construct a decision tree and print it.
mod_tree <- rpart(Y ~ . , data = train_data)
print(mod_tree)## n= 3983
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 3983 33370.59000 5.793440
## 2) https< 0.5 3196 21102.42000 4.965168
## 4) twacc>=0.5 2134 11327.23000 4.102699
## 8) rt< 0.5 2000 8134.59700 3.820299
## 16) interviewed< 0.5 1972 7712.69000 3.769130 *
## 17) interviewed>=0.5 28 53.10444 7.424061 *
## 9) rt>=0.5 134 652.53380 8.317631 *
## 5) twacc< 0.5 1062 4998.10000 6.698227
## 10) http>=0.5 147 218.91660 4.493694 *
## 11) http< 0.5 915 3949.99400 7.052398 *
## 3) https>=0.5 787 1171.62600 9.157042 *And here is the plot:
rpart.plot(mod_tree)
We can get access to the guts of the tree model to be able to extract particular splits and deviations at those splits:
mod_tree$frameQuestion 5
Compute the RSS in two ways: by aggregating deviance of all leaf nodes and as a total training error.
The RSS of the entire dataset is deviance of the root node, Below is the RSS of the decision tree.
sum(mod_tree$frame$dev[mod_tree$frame$var == "<leaf>"])## [1] 13758.86And here is the RSS:
tse <- function(x, y) {
sum((x-y)^2)
}
mod_tree %>%
predict(train_data) %>%
tse(train_data$Y)## [1] 13758.86Tree pruning (regularization)
To prune a decision tree, we first grow it and then we choose a
sub-tree minimizing the regularized loss function \[
\sum_{m=1}^{|T|}\sum_{x^{i}\in R_m}(y^{i}-\hat{y}_{R_i})^2
+\alpha|T|
\] Below is the plot of cross-validation error vs \(\alpha\) (parameter cp):
plotcp(mod_tree)
In this case, we do not need to prune the tree. But for your future reference, here is how we do it. First, we will look at the table
mod_tree$cptable## CP nsplit rel error xerror xstd
## 1 0.33252464 0 1.0000000 1.0003914 0.01369043
## 2 0.14315276 1 0.6674754 0.6678108 0.01382708
## 3 0.07611784 2 0.5243226 0.5247545 0.01248580
## 4 0.02484792 3 0.4482048 0.4487915 0.01095215
## 5 0.01105171 4 0.4233568 0.4240781 0.01072747
## 6 0.01000000 5 0.4123051 0.4160796 0.01068395The optimal value of \(\alpha\) is
opt_cp <- mod_tree$cptable[which.min(mod_tree$cptable[ , 'xerror']) , 'CP']
opt_cp## [1] 0.01And now we prune the tree:
mod_tree_pruned <- prune.rpart(mod_tree, opt_cp)
mod_tree_pruned## n= 3983
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 3983 33370.59000 5.793440
## 2) https< 0.5 3196 21102.42000 4.965168
## 4) twacc>=0.5 2134 11327.23000 4.102699
## 8) rt< 0.5 2000 8134.59700 3.820299
## 16) interviewed< 0.5 1972 7712.69000 3.769130 *
## 17) interviewed>=0.5 28 53.10444 7.424061 *
## 9) rt>=0.5 134 652.53380 8.317631 *
## 5) twacc< 0.5 1062 4998.10000 6.698227
## 10) http>=0.5 147 218.91660 4.493694 *
## 11) http< 0.5 915 3949.99400 7.052398 *
## 3) https>=0.5 787 1171.62600 9.157042 *Variable importance
Variable importance is measured as a total drop in residual sum of
squares due to splits in each variable. Here is how we calculate it with
the function varImp from caret:
varImp(mod_tree_pruned) %>%
arrange(-Overall) %>%
head(n = 20)And here is how we can extract it directly from the model:
mod_tree_pruned$variable.importance## https twacc rt
## 11096.542899 4796.047652 2540.097216
## http interviewed obama
## 829.189620 368.802040 206.917343
## join makeamericagreatagain whitehouse
## 197.397205 183.297405 165.278749
## democrats cont hillary
## 155.964861 141.018643 121.451484
## honor crooked obamacare
## 84.598802 80.967656 80.967656
## amp donald gop
## 18.955949 5.640746 5.640746
## hotel national
## 5.640746 5.640746Test error
Finally, let us construct predictions and report the test root mean
squared error calculated in caret (here, we did not create
our own version of the very simple RMSE function because we
are loading caret anyway).
mod_tree %>%
predict(test_data) %>%
RMSE(test_data$Y)## [1] 1.865645Random forest
We will use the library ranger because it allows tuning
more than just one hyperparameter.
Training
First, we train a random forest model and print the result. Note that
by default, it tries three values of mtry (the number of
predictors allowed at each step), \(2\), \(p/2\) and \(p\) and two values of
splitrule (this is something that we did not cover in class
and you can either ignore it or google what it is) and only one value of
min.node.size, 5.
Here we set train control to oob, i.e., out-of-bag error
(with 5-fold cross validation, training time will be 5 times slower) and
num.trees to 50 (with default value of 500, training time
will be 10 times slower).
set.seed(100)
mod_rf <- train(Y ~ . , data = train_data, method = "ranger",
num.trees = 50,
importance = 'impurity',
trControl = trainControl("oob"))
print(mod_rf)## Random Forest
##
## 3983 samples
## 177 predictor
##
## No pre-processing
## Resampling results across tuning parameters:
##
## mtry splitrule RMSE Rsquared MAE
## 2 variance 2.020853 0.6289601 1.702157
## 2 extratrees 2.382233 0.4836240 2.039927
## 89 variance 1.626012 0.6880502 1.171465
## 89 extratrees 1.614116 0.6919872 1.161927
## 177 variance 1.644185 0.6826860 1.180640
## 177 extratrees 1.617662 0.6915464 1.167122
##
## Tuning parameter 'min.node.size' was held constant at a value of 5
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 89, splitrule = extratrees
## and min.node.size = 5.Predictions
Here we will construct predictions and report the test error
mod_rf %>%
predict(test_data) %>%
RMSE(test_data$Y)## [1] 1.582583Variable importance
Here is variable importance for top 20 most important variables:
varImp(mod_rf)## ranger variable importance
##
## only 20 most important variables shown (out of 177)
##
## Overall
## https 100.000
## twacc 75.702
## realdonaldtrump 25.702
## rt 24.511
## http 19.511
## hillary 5.999
## democrats 4.927
## fake 3.949
## htag 3.649
## trump 3.549
## interviewed 3.465
## cnn 3.395
## america 3.206
## people 2.688
## media 2.408
## luck 2.233
## obama 2.104
## apprenticenbc 2.059
## amp 2.043
## foxnews 1.864If we want all, here is how we can get the full information (we only printed the top 10 most important variables, but you can easily print the whole vector):
var_importance <- mod_rf$finalModel$variable.importance %>%
sort(decreasing = TRUE)
var_importance %>% head(10)## https twacc realdonaldtrump rt http
## 7721.8103 5847.6366 1991.0417 1899.1366 1513.4958
## hillary democrats fake htag trump
## 471.2336 388.6107 313.1711 289.9933 282.3058Here is the built-in plot of variable importance:
varImp(mod_rf) %>%
plot(top = 10)
And here is how you can make a custom plot with
ggplot2:
var_importance <- mod_rf$finalModel$variable.importance %>%
sort(decreasing = TRUE) %>% head(10)
var_importance %>%
enframe(name = "variable", value = "importance") %>%
mutate(variable = fct_reorder(variable, importance)) %>%
ggplot(aes(variable, importance)) +
geom_col() +
coord_flip() 
Tuning random forest
Random forest has a number of hyperparameters that can be tuned with OOB error (faster) or with cross-validation (slower).
Random forest models may take a lot of time to train. For the sake of saving time, we will create a mini-version of our dataset to demonstrate the tuning process.
set.seed(199)
mini_data <- train_data %>% slice_sample(n = 500)
rfGrid <- expand.grid(mtry = c(10, 20, 30, 40, 50, 60),
min.node.size = c(5, 10, 20, 40),
splitrule = "variance")
mod_rf_tune <- train(Y ~ . , data = mini_data, method = "ranger",
num.trees = 100,
importance = 'impurity',
tuneGrid = rfGrid,
trControl = trainControl("oob"))
mod_rf_tune## Random Forest
##
## 500 samples
## 177 predictors
##
## No pre-processing
## Resampling results across tuning parameters:
##
## mtry min.node.size RMSE Rsquared MAE
## 10 5 1.819839 0.6350919 1.411404
## 10 10 1.825282 0.6320670 1.416271
## 10 20 1.822569 0.6342226 1.423065
## 10 40 1.830677 0.6339977 1.432435
## 20 5 1.819014 0.6309308 1.339478
## 20 10 1.796023 0.6397083 1.337772
## 20 20 1.777173 0.6480410 1.315884
## 20 40 1.791350 0.6417560 1.342114
## 30 5 1.790549 0.6450794 1.291179
## 30 10 1.837643 0.6255028 1.339370
## 30 20 1.816696 0.6339217 1.317052
## 30 40 1.829556 0.6276364 1.340923
## 40 5 1.841165 0.6259442 1.329190
## 40 10 1.844718 0.6254696 1.334618
## 40 20 1.816654 0.6355616 1.311711
## 40 40 1.824669 0.6293592 1.343431
## 50 5 1.835547 0.6302223 1.319438
## 50 10 1.831429 0.6320274 1.326252
## 50 20 1.815427 0.6375645 1.294562
## 50 40 1.839853 0.6263216 1.332915
## 60 5 1.851863 0.6273492 1.329219
## 60 10 1.840637 0.6294377 1.314511
## 60 20 1.828681 0.6335251 1.315683
## 60 40 1.854598 0.6205367 1.342260
##
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 20, splitrule = variance
## and min.node.size = 20.The tuning process can be plotted to get a better picture of what is going on:
plot(mod_rf_tune)
The optimal values of the hyperparameters are
mod_rf_tune$bestTuneNow we will retrain the model on the whole training data with these values
mod_rf_tuned <- train(Y ~ . , data = train_data, method = "ranger",
num.trees = 100,
importance = 'impurity',
tuneGrid = expand.grid(mod_rf_tune$bestTune),
trControl = trainControl("oob"))
mod_rf_tuned## Random Forest
##
## 3983 samples
## 177 predictor
##
## No pre-processing
## Resampling results:
##
## RMSE Rsquared MAE
## 1.560631 0.7093683 1.142877
##
## Tuning parameter 'mtry' was held constant at a value of 20
## Tuning
## parameter 'splitrule' was held constant at a value of variance
##
## Tuning parameter 'min.node.size' was held constant at a value of 20Here is the test error. It should be a bit better than the first version of random forest.
mod_rf_tuned %>%
predict(test_data) %>%
RMSE(test_data$Y)## [1] 1.548515Question 6
Extract only 70 most important variables and train a new random forest with 50 trees on just those 70 variables. Tune hyperparameters using the following values:
mtry= 3, 5, 8, 15min.node.size= 2, 4, 8, 16, 32
Retrain random forest with the optimal values of hyperparameters but with 500 trees instead of 50 trees.
Report the test error of the final model.
top_variables <- mod_rf_tuned$finalModel$variable.importance %>%
sort(decreasing = TRUE) %>% head(70) %>% names
t_data <- train_data %>%
select(all_of(top_variables), Y)
v_data <- test_data %>%
select(all_of(top_variables), Y)
rfGrid_2 <- expand.grid(mtry = c(3, 5, 8, 15),
min.node.size = c(2, 4, 8, 16, 32),
splitrule = "variance")
mod_rf_topvar <- train(Y ~ . , data = t_data, method = "ranger",
num.trees = 50,
importance = 'impurity',
tuneGrid = expand.grid(rfGrid_2),
trControl = trainControl("oob"))
mod_rf_topvar## Random Forest
##
## 3983 samples
## 70 predictor
##
## No pre-processing
## Resampling results across tuning parameters:
##
## mtry min.node.size RMSE Rsquared MAE
## 3 2 1.700490 0.6824198 1.353845
## 3 4 1.716438 0.6799758 1.380109
## 3 8 1.717438 0.6736306 1.382346
## 3 16 1.730252 0.6744493 1.391610
## 3 32 1.695673 0.6846646 1.362388
## 5 2 1.595151 0.7019814 1.221286
## 5 4 1.598809 0.7009338 1.216815
## 5 8 1.592570 0.7023447 1.209848
## 5 16 1.596552 0.7014013 1.215199
## 5 32 1.594411 0.7030912 1.225699
## 8 2 1.561509 0.7094065 1.154402
## 8 4 1.555794 0.7116010 1.149088
## 8 8 1.568004 0.7069067 1.158538
## 8 16 1.556805 0.7112922 1.154281
## 8 32 1.565099 0.7081617 1.162558
## 15 2 1.572930 0.7050435 1.148112
## 15 4 1.562084 0.7090474 1.136686
## 15 8 1.566665 0.7073324 1.140338
## 15 16 1.561109 0.7092797 1.140003
## 15 32 1.571879 0.7051903 1.150333
##
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 8, splitrule = variance
## and min.node.size = 4.The optimal values of hyperparameters are
mod_rf_topvar$bestTunemod_rf_final <- train(Y ~ . , data = t_data, method = "ranger",
num.trees = 500,
importance = 'impurity',
tuneGrid = expand.grid(mod_rf_topvar$bestTune),
trControl = trainControl("oob"))And the error of the final model is
mod_rf_final %>%
predict(test_data) %>%
RMSE(test_data$Y)## [1] 1.54252Answers
Here are the answers:
- INSERT LINK
Declaration of AI Usage
I used Chat GPT to update codes from 2024 to make them more
tidyverse-styled.