9_predict.qmd

---
execute:
  echo: true
---

::: {.content-visible when-format="pdf"}
```{=latex}
\setDOI{10.4324/9781003393764.9}
\thispagestyle{chapterfirstpage}
```
:::

# Predict {#sec-predict-chapter}

```{r}
#| label: setup-options
#| child: "../_common.qmd"
#| cache: false
```

::: {.callout}
**{{< fa regular list-alt >}} Outcomes**

- Identify the research goals of predictive data analysis
- Describe the workflow for predictive data analysis
- Recognize quantitative and qualitative methods for evaluating predictive models
:::

In this chapter, I introduce supervised learning as an approach to text analysis. Supervised learning aims to establish a relationship between a target (or outcome) variable and a set of feature variables derived from text data. By leveraging this relationship, statistical generalizations (models) can be created to accurately predict values of the target variable based on the values of the feature variables. Throughout the chapter, we explore practical tasks and theoretical applications of statistical learning in text analysis.

::: {.callout}
**{{< fa terminal >}} Lessons**

**What**: Advanced Visualization\
**How**: In an R console, load {swirl}, run `swirl()`, and follow prompts to select the lesson.\
**Why**: To dive deeper into {ggplot2} to enhance visual summaries and provide an introduction to {factoextra} and {ggfortify} that extend {ggplot2} capabilities to model objects.
:::

## Orientation {#sec-predict-orientation}

Predictive data analysis (PDA)\index{predictive data analysis (PDA)} is a powerful analysis method for making predictions about new or future data based on patterns in existing data. PDA is a type of supervised learning\index{supervised learning}, which means that it involves training a model on a labeled dataset where the input data and desired output are both provided. The model is able to make predictions\index{text regression} or classifications\index{text classification} based on the input data by learning the relationships between the input and output data. Supervised machine learning is an important tool for linguists studying language and communication, as it allows us to analyze language data to identify patterns or trends in language use, assess hypotheses, and prescribe actions.\index{machine learning}

The approach to conducting predictive analysis shares some commonalities with exploratory data analysis (@sec-explore-orientation) (as well as inferential analysis @sec-infer-chapter), but there are also some key differences. Consider the workflow in @tbl-predict-workflow.

<!-- Workflow -->

::: {#tbl-predict-workflow tbl-colwidths="[5, 15, 80]"}

| Step | Name | Description |
|------|-------------|-------------|
| 1 | Identify | Consider the research question and aim and identify relevant variables |
| 2 |  | Split the data into representative training and testing sets |
| 3 |  | Apply variable selection and engineering procedures |
| 4 | Inspect | Inspect the data to ensure that it is in the correct format and that the training and testing sets are representative of the data |
| 5 | Interrogate | Train and evaluate the model on the training set, adjusting models or hyperparameters as needed, to produce a final model |
| 6 | (Optional) Iterate | Repeat steps 3-5 to select new variables, models, hyperparameters |
| 7 | Interpret | Interpret the results of the final model in light of the research question or hypothesis |

Workflow for predictive data analysis
:::

Focusing on the overlap with other analysis methods, we can see some fundamental steps such as identifying relevant variables, inspecting the data, interrogating the data, and interpreting the results. And if our research aim is exploratory in nature, iteration may also be a part of the workflow.\index{research aim}

There are two main differences, however, between the PDA and the EDA workflow we discussed in @sec-explore-chapter\index{exploratory data analysis (EDA)}. The first is that PDA requires partitioning the data into training and testing sets. The training set is used to develop the model, and the testing set is used to evaluate the model's performance. This strategy is used to ensure that the model is robust and generalizes well to new data. It is well known, and makes intuitive sense, that using the same data to develop and evaluate a model likely will not produce a model that generalizes well to new data. This is because the model will have potentially conflated the nuances of the data ('the noise') with any real trends ('the signal') and therefore will not be able to generalize well to new data. This is called **overfitting**\index{overfitting} and by holding out a portion of the data for testing, we can evaluate the model's performance on data that it has not seen before and therefore get a more accurate estimate of the generalizable trends in the data.

::: {.callout .halfsize}
**{{< fa medal >}} Dive deeper**

Prediction modeling is a hot topic. Given the potential to make actionable predictions about future outcomes, it attracts a lot of attention from organizations which aim to leverage data to make informed decisions. It's use in research is also growing beyond the development of better models and using predictive models to address research questions and hypotheses.

We will apply predictive modeling in the context of language data as a semi-inductive method. However, it is also increasingly used in hypothesis testing scenarios, see @Gries2014, @Deshors2016, and @Baayen2011 for examples\index{Gries}\index{Baayen}.
:::

Another procedure to avoid the perils of overfitting, is to use resampling methods as part of the model evaluation on the training set. Resampling\index{resampling} is the process of repeatedly drawing samples from the training set and evaluating the model on each sample. The two most common resampling methods are **bootstrapping**\index{bootstrapping} (resampling with replacement) and **cross-validation**\index{cross-validation} (resampling without replacement). The performance of these multiple models is summarized and the error between them is assessed. The goal is to minimize the performance differences between the models while maximizing the overall performance. These measures go a long way to avoiding overfitting and therefore maximizing the chance that the training phase will produce a model which is robust at the testing phase.

The second difference, not reflected in the workflow but inherent in predictive analysis, is that PDA requires a fixed outcome variable\index{outcome variable}. This means that the outcome variable must be defined from the outset and cannot be changed during the analysis. Furthermore, the informational nature of the outcome variable will dictate what type of algorithm we choose to interrogate the data and how we will evaluate the model's performance.

If the outcome is categorical in nature, we will use a **classification algorithm**\index{classification algorithms} (*e.g.* logistic regression, naive Bayes, *etc.*). Classification evaluation metrics include accuracy, precision, recall, and F1 scores (a metric which balances precision and recall) which can be derived from and visualized in a cross-tabulation of the predicted and actual outcome values\index{research interpretation}.

If the outcome is numeric in nature, we will use a **regression algorithm**\index{regression algorithms} (*e.g.* linear regression, support vector regression, *etc.*). Since the difference between prediction and actual values is numeric, metrics that quantify numerical differences, such as root mean square error (RMSE) or $R^2$, are used to evaluate the model's performance\index{research interpretation}.

The evaluation of the model is quantitative on the one hand, but it is also qualitative in that we need to consider the implications of the model's performance in light of the research question or hypothesis. Furthermore, depending on our research question we may be interested in exploring the features that are most important to the model's performance. This is called **feature importance**\index{feature importance} and can be derived from the model's coefficients or weights. Notably, however, some of the most powerful models in use today, such as deep neural networks, are not easily interpretable and therefore feature importance is not easily derived. This is something to keep in mind when considering the research question and the type of model that will be used to address it\index{research interpretation}.

## Analysis {#sec-predict-analysis}

<!-- Goals of this section -->

In this section, we now turn to the practical application of predictive data analysis. The discussion will be separated into classification and regression tasks, as model selection and evaluation procedures differ between the two. For each task, we will frame a research goal and work through the process of building a predictive model to address that goal. Along the way we will cover concepts and methods that are common to both classification and regression tasks and specific to each.

<!-- Research questions -->

To frame our analyses, we will posit research aimed at identifying language usage patterns in second language use, one for a classification task and one for a regression task. Our first research question will be to assess whether Spanish language use can be used to predict natives and L1 English learners (categorical). Our second research question will be to gauge the extent to which the L1 English learners' Spanish language placement test scores (numeric) can be predicted based on their language use.

```{r}
#| label: predict-class-curate-data-run
#| eval: false
#| echo: false

# Read in the datasets
learners_tbl <-
  read_csv("data/cedel2_learners.csv") |> # read in the learners dataset outcome variable |>
  mutate(subcorpus = "Learner") |>
  select(subcorpus, place_score = placement_test_score_percent, proficiency, text) # select variables

natives_tbl <-
  read_csv("data/cedel2_natives.csv") |> # read in the natives dataset
  mutate(place_score = NA) |> # create an outcome variable
  mutate(proficiency = "Native") |> # create a proficiency variable
  select(subcorpus, place_score, proficiency, text) # select the text column

# Combine the datasets by row
cedel_tbl <-
  bind_rows( # combine the datasets by row
    learners_tbl,
    natives_tbl
  ) |>
  mutate(doc_id = row_number()) |> # create a document id
  select(doc_id, subcorpus, place_score, proficiency, text) # select variables

# Write to disk
cedel_tbl |>
  write_csv("data/cedel2_transformed.csv")

# Create data dictionary
create_data_dictionary(
  data = cedel_tbl,
  file_path = "data/cedel2_transformed_dd.csv",
  model = "gpt-3.5-turbo",
  grouping = "subcorpus"
)
```

```{r}
#| label: tbl-predict-cedel2-data-dictionary
#| tbl-cap: "Data dictionary for the CEDEL2 corpus"
#| tbl-colwidths: [15, 15, 15, 55]
#| echo: false

# Read in the data dictionary for the CEDEL2 corpus
read_csv("data/cedel2_transformed_dd.csv") |>
  tt(width = 1)
```

We will use data from the CEDEL2 corpus [@Lozano2022]\index{Corpus Escrito del Español como L2 (CEDEL)}. We will include a subset of the variables from this data that are relevant to our research questions\index{research question}. The data dictionary for this dataset is seen in @tbl-predict-cedel2-data-dictionary.

Let's go ahead and read the transformed dataset and preview it in @exm-predict-cedel-read.

::: {#exm-predict-cedel-read}
```r
# Read in the dataset
cedel_tbl<-
  read_csv("../data/cedel2/cedel2_transformed.csv")

# Preview
glimpse(cedel_tbl)
```

```{r}
#| label: predict-cedel-read
#| echo: false

# Read in the dataset
cedel_tbl<-
  read_csv("data/cedel2_transformed.csv")

# Preview
glimpse(cedel_tbl)
```
\index{R packages!readr}
\cindex{read_csv()}\cindex{glimpse()}
:::

The output of @exm-predict-cedel-read provides some structural information about the dataset, number of rows and columns as well as variable types.

```{r}
#| label: predict-cedel-diagnostics
#| include: false

# Load packages
library(skimr)
library(janitor)

# Diagnostics
cedel_tbl |>
  skim()

# There appear to be 1051 NA values for place_score. This is likely due to the fact that the natives dataset does not have a place_score variable, let's take a look.

cedel_tbl |>
  filter(is.na(place_score)) |>
  select(subcorpus, place_score, proficiency) |>
  head()

# Yes, placement score is NA for the natives dataset.

# Next, let's look at how the proficiency variable is distributed.

cedel_tbl |>
  tabyl(proficiency) |>
  adorn_pct_formatting(digits = 1)

# Lower beginner and upper advanced have very few observations. But since this will not be a class, we can proceed with the analysis.

# Next, let's look at the distribution of the place_score variable.

cedel_tbl |>
  ggplot(aes(x = place_score)) +
  geom_histogram(aes(y = after_stat(density))) +
  geom_density(alpha = 0.5) +
  labs(title = "Distribution of placement scores",
       x = "Placement score",
       y = "Frequency")

# Some left skew, with an apparent bimodal distribution. Since this is only a learner dataset, we can assume that the bimodal distribution is due to the fact that the learners are at different levels of proficiency. Let's explore this further.

cedel_tbl |>
  ggplot(aes(x = place_score, fill = proficiency)) +
  # geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.5, binwidth = 2) +
  geom_density(alpha = 0.5) +
  labs(title = "Distribution of placement scores by proficiency",
       x = "Placement score",
       y = "Frequency")

# Looks like our two peaks are due to the lower and upper proficiency levels. This is good to know and will be useful for our analysis.

# Let's gauge the variablity of the place_score variable by proficiency. A boxplot will be useful for this.

# Before we do that, let's make proficiency and subscorpus a factor variable in the order we want.

cedel_tbl <-
  cedel_tbl |>
  mutate(proficiency =
    factor(proficiency,
          levels =
            c("Lower beginner",
              "Upper beginner",
              "Lower intermediate",
              "Upper intermediate",
              "Lower advanced",
              "Upper advanced",
              "Native")
            ),
          subcorpus = factor(subcorpus, levels = c("Learner", "Native")
          )
        )

cedel_tbl |>
  ggplot(aes(x = proficiency, y = place_score, fill = proficiency)) +
  geom_boxplot() +
  labs(title = "Variability of placement scores by proficiency",
       x = "Proficiency",
       y = "Placement score")

# Finally, let's look at the distribution of the subcorpus variable.

cedel_tbl |>
  tabyl(subcorpus) |>
  adorn_pct_formatting(digits = 1)
```

After I performed some diagnostics and made some adjustments based on a descriptive assessment\index{descriptive assessment}, the dataset is in good order to proceed with the analysis. I updated the variables `subcorpus` and `proficiency` as factor variables and ordered them in a way that makes sense for the analysis. The `place_score` variable is distributed well across the proficiency levels. The `subcorpus` variable is less balanced, with around 65% of the texts being from learners. This is not a problem, but it is something to keep in mind when building and interpreting the predictive models.

We will be using the Tidymodels framework\index{Tidymodels} in R to perform these analyses. {tidymodels} is a meta-package, much like {tidyverse}, that provides a consistent interface for machine learning modeling. Some key packages unique to {tidymodels} are {recipes}, {parsnip}, {workflows}, and {tune}. {recipes} includes functions for pre-processing and engineering features. {parsnip} provides a consistent interface for specifying modeling algorithms. {worflows} allows us to combine recipes and models into a single pipeline. Finally, {tune} give us the ability to evaluate and adjust, or 'tune', the parameters of models.

Since we are using text data, we will also be using {textrecipes} which makes various functions available for pre-processing text including extracting and engineering features.

Let's go ahead and do the setup, loading the necessary packages, seen in @exm-predict-packages-data.

::: {#exm-predict-packages-data}
```r
# Load packages
library(tidymodels)   # modeling metapackage
library(textrecipes)  # text pre-processing

# Prefer tidymodels functions
tidymodels_prefer()
```

```{r}
#| label: predict-packages-data-run
#| echo: false

# Load packages
library(tidymodels)   # modeling metapackage
library(textrecipes)  # text pre-processing
library(janitor)      # data inspection

# Prefer tidymodels functions
tidymodels_prefer()
```
\index{R packages!tidymodels}\index{R packages!textrecipes}
\cindex{tidymodels_prefer()}\cindex{library()}
:::

### Text classification {#sec-predict-text-classification}

<!-- Research question: outcome and features -->

The goal of text classification analysis is to develop a model that can accurately label text samples as either native or learner. This is a binary classification problem. We will approach this problem from an exploratory perspective, and therefore our aim is to identify features from the text that best distinguish between the two classes and explore the features that are most important to the model's performance.\index{text classification}

Let's modify the data frame to include only the variables we need for this analysis, assigning it to `cls_tbl`. In the process, we will rename the `subcorpus` variable to `outcome` to reflect that it is the outcome variable. This is seen in @exm-predict-class-data.

::: {#exm-predict-class-data}
```{r}
#| label: predict-class-data

# Rename subcorpus to outcome
cls_tbl <-
  cedel_tbl |>
  select(outcome = subcorpus, proficiency, text)
```
:::

<!-- Step 1: identify features -->

Let's begin the workflow from @tbl-predict-workflow by identifying the features that we will use to classify the texts. There may be many features that we could use. These could be features derived from raw text (*e.g.* characters, words, ngrams, *etc.*), feature vectors (*e.g.* word embeddings), or meta-linguistic features (*e.g.* part-of-speech tags, syntactic parses, or semantic features) that have been derived from these through manual or automatic annotation.\index{feature selection}

If as part of our research question\index{research question} the types of features are included, then we should proceed toward deriving those features. If not, a simple approach is to use words as the predictor features\index{predictor features}. This will serve as a baseline for more complex models, if necessary.

This provides us the linguistic unit we will use, but we still need to decide how to operationalize what we mean by 'use' in our research statement. Do we use raw token counts? Do we use normalized frequencies? Do we use some type of weighting scheme? These are questions that we need to consider as we embark on this analysis. Since we are exploring, we can use trial-and-error or consider the implications of each approach and choose the one that best fits our research question ---or both.

Let's approach this with a bit more nuance as we already have some domain knowledge about word use. First, we know that the frequency distribution of words is highly skewed\index{skewed distribution}, meaning that a few words occur very frequently and most words occur very infrequently. Second, we know that the most frequent words in a language are often function words\index{function words} (*e.g.* 'the', 'and', 'of', *etc.*) and that these words are not very informative for distinguishing between classes of texts. Third, we know that comparing raw counts across texts conflates the influence text class lengths.

With these considerations in mind, we will tokenize by words and apply a metric known as the **term-frequency inverse-document frequency** ($tf$-$idf$)\index{term-frequency inverse-document frequency (tf-idf)}. The $tf$-$idf$ measure, as the name suggests, is the product of $tf$ and $idf$ for each term. In effect, it produces a weighting scheme, which will downweight words that are common across all documents and upweight words that are unique to a document. It also mitigates the varying lengths of the documents. This is a common approach in text classification and is a good starting point for our analysis.\index{feature engineering}

<!-- Step 2: Initial split -->

With our features and engineering approach identified, we can move on to step 2 of our workflow and split the data into training and testing sets. We make the splits to our data at this point to draw a line in the sand between the data we will use to train the model and the data we will use to test the model. A typical approach in supervised machine learning\index{supervised learning} is to allocate around 75-80% of the data to the training set and the remaining 20-25% to the testing set, depending on the number of observations\index{observations}. We have `r nrow(cls_tbl)` observations in our dataset, so we can allocate 80% of the data to the training set and 20% of the data to the testing set.

In @exm-predict-class-split, we will use the `initial_split()` function from {rsample} to split the data into training and testing sets. The `initial_split()` function takes a data frame and a proportion and returns a `split` object which contains the training and testing sets. We will use the `strata` argument to stratify the data by the `outcome` variable. This will ensure that the training and testing sets have the same proportion of native and learner texts.\index{outcome variable}

\pagebreak

::: {#exm-predict-class-split}
```{r}
#| label: predict-class-split


set.seed(123) # for reproducibility

# Split the data into training and testing sets
cls_split <-
  initial_split(
    data = cls_tbl,
    prop = 0.8,
    strata = outcome
  )

# Create training set
cls_train <- training(cls_split)  # 80% of data

# Create testing set
cls_test <- testing(cls_split)    # 20% of data
```
\index{R packages!rsample}
\cindex{initial_split()}\cindex{training()}\cindex{testing()}\cindex{set.seed()}
:::

A confirmation of the distribution of the data across the training and testing sets as well as a breakdown of the outcome variable, created by {janitor}'s `tabyl()` function, can be seen in @exm-predict-class-split-tabyl.

::: {#exm-predict-class-split-tabyl}
```{r}
#| label: predict-class-split-tabyl
#| results: hold

# View the distribution
# Training set
cls_train |>
  tabyl(outcome) |>
  adorn_totals("row") |>
  adorn_pct_formatting(digits = 1)

# Testing set
cls_test |>
  tabyl(outcome) |>
  adorn_totals("row") |>
  adorn_pct_formatting(digits = 1)
```
\index{R packages!janitor}
\cindex{tabyl()}\cindex{adorn_totals()}\cindex{adorn_pct_formatting()}
:::

We can see that the split was successful. The training and testing sets have very similar proportions of native and learner texts.

<!-- Step 3: Integrate: plan to select and engineer features -->

We are now ready to create a 'recipe', step 3 in our analysis. A recipe is Tidymodels terminology for a set of instructions or blueprint which specifies the outcome variable and the predictor variable and determines how to pre-process and engineer the feature variables.

We will use the `recipe()` function from {recipes} to create the recipe. The `recipe()` function minimally takes a formula and a data frame and returns a `recipe` object. **R formulas**\index{R formula} provide a way to specify relationships between variables and are used extensively in R data modeling. Formulas specify the outcome variable ($y$) and the predictor variable(s)\index{predictor variables} ($x_1 .. x_n$). For example, `y ~ x` can be read as "y as a function of x". In our particular case, we will use the formula `outcome ~ text` to specify that the outcome variable is the `outcome` variable and the predictor variable is the `text` variable. The code is seen in @exm-predict-class-recipe.\cindex{\textasciitilde}

::: {#exm-predict-class-recipe}
```{r}
#| label: predict-class-recipe

# Create a recipe
base_rec <-
  recipe(
    formula = outcome ~ text, # formula
    data = cls_train
    )

# Preview
base_rec
```
\index{R packages!recipes}
\cindex{recipe()}

```
── Recipe ─────────────────────────────────────────

── Inputs
Number of variables by role
outcome:   1
predictor: 1
```

:::

The recipe object at this moment contains just one instruction, what the variables are and what their relationship is.

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

R formulas are a powerful way to specify relationships between variables and are used extensively in data modeling including exploratory, predictive, and inferential analysis. The basic formula syntax is `y ~ x` where `y` is the outcome variable and `x` is the feature variable. The formula syntax can be extended to include multiple feature variables, interactions, and transformations. For more information on R formulas, see R for Data Science [@Wickham2017]
:::

{recipes} provides a wide range of `step_*()` functions which can be applied to the recipe to specify how to engineer the variables in our recipe call. These include functions to scale (*e.g* `step_center()`, `step_scale()`, *etc.*) and transform (*e.g.* `step_log()`, `step_pca()`, *etc.*) numeric variables, and functions to encode (*e.g.* `step_dummy()`, `step_labelencode()`, *etc.*) categorical variables.\index{feature engineering}

These step functions are great when we have selected the variables we want to use in our model and we want to engineer them in a particular way. In our case, however, we need to derive features from the text in the `text` column of datasets before we engineer them.

To ease this process, {textrecipes} provides a number of step functions for pre-processing text data. These include functions to tokenize\index{tokens} (*e.g.* `step_tokenize()`), remove stop words\index{stopwords} (*e.g.* `step_stopwords()`), and to derive meta-features\index{metadata} (*e.g.* `step_lemma()`, `step_stem()`, *etc.*)[^1]. Furthermore, there are functions to engineer features in ways that are particularly relevant to text data, such as feature frequencies and weights (*e.g.* `step_tf()`, `step_tfidf()`, *etc.*) and token filtering (*e.g.* `step_tokenfilter()`).\index{standardization}

::: {.callout .halfsize}
**{{< fa medal >}} Dive deeper**
For other tokenization strategies and feature engineering methods, see {textrecipes} documentation [@R-textrecipes]. There are, however, packages which provide integration with `textrecipes` for other languages, for example, {washoku} for Japanese text processing [@R-washoku].
:::

[^1]: Note that functions for meta-features require more sophisticated text analysis software to be installed on the computing environment (e.g. {spacyr} for `step_lemma()`, `step_pos()`, *etc.*). See {textrecipes} documentation for more information.

So let's build on our basic recipe `cls_rec` by adding steps relevant to our task. To extract our features, we will use the `step_tokenize()` function to tokenize the text into words. The default behavior of the `step_tokenize()` function is to tokenize the text into words, but other token units can be derived and various options can be added to the function call (as {tokenizers} is used under the hood). Adding the `step_tokenize()` function to our recipe is seen in @exm-predict-class-recipe-tokenize.

::: {#exm-predict-class-recipe-tokenize}
```{r}
#| label: predict-class-recipe-tokenize

# Add step to tokenize the text
cls_rec <-
  base_rec |>
  step_tokenize(text) # tokenize

# Preview
cls_rec
```
\index{R packages!tokenizers}\index{R packages!textrecipes}
\cindex{step_tokenize()}

::: {.content-visible when-format="pdf"}
\vspace{1em}
:::

```
── Recipe ─────────────────────────────────────────

── Inputs
Number of variables by role
outcome:   1
predictor: 1

── Operations
• Tokenization for: text
```
:::

The recipe object `cls_rec` now contains two instructions, one for the outcome variable and one for the feature variable. The feature variable instruction specifies that the text should be tokenized into words.

We now need to consider how to engineer the word features. If we add `step_tf()` we will get a matrix of token counts by default\index{matrix}, with the option to specify other weights. The step function `step_tfidf()` creates a matrix of term frequencies\index{frequency} weighted by inverse document frequency\index{inverse document frequency}.

We decided in step 1 that we will start with $tf$-$idf$\index{term-frequency inverse-document frequency (tf-idf)}, so we will add `step_tfidf()` to our recipe. This is seen in @exm-predict-class-recipe-tfidf.

::: {#exm-predict-class-recipe-tfidf}
```{r}
#| label: predict-class-recipe-tfidf

# Add step to tokenize the text
cls_rec <-
  cls_rec |>
  step_tfidf(text, smooth_idf = FALSE)

# Preview
cls_rec
```
\index{R packages!textrecipes}
\cindex{step_tfidf()}
```
── Recipe ─────────────────────────────────────────
Number of variables by role
outcome:   1
predictor: 1

── Operations
• Tokenization for: text
• Term frequency-inverse document frequency with: text
```
:::

<!-- Step 4: Inspect:  -->

<!-- [ ] tmp fmt -->
\pagebreak

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

The `step_tfidf()` function by default adds a *smoothing term* to the inverse document frequency ($idf$) calculation\index{inverse document frequency}. This setting has the effect of reducing the influence of the $idf$ calculation. Thus, terms that appear in many (or all) documents will not be downweighted as much as they would be if the smoothing term was not added. For our purposes, we want to downweight or eliminate the influence of the most frequent terms, so we will set `smooth_idf = FALSE`.
:::

To make sure things are in order and that the recipe performs as expected, we can use the functions `prep()` and `bake()` to inspect the recipe. The `prep()` function takes a recipe object and a data frame and returns a `prep` object. The `prep` object contains the recipe and the data frame with the feature variables engineered according to the recipe. The `bake()` function takes a `prep` object and an optional new dataset to apply the recipe to. If we only want to see the application to the training set, we can use the `new_data = NULL` argument.

In @exm-predict-class-recipe-prep, we use the `prep()` and `bake()` functions to create a data frame with the feature variables. We can then inspect the data frame to see if the recipe performs as expected.

::: {#exm-predict-class-recipe-prep}
```{r}
#| label: predict-class-recipe-prep

# Prep and bake
cls_bake <-
  cls_rec |>
  prep() |> # create a prep object
  bake(new_data = NULL) # apply to training set

# Preview
dim(cls_bake)
```
\index{R packages!recipes}
\cindex{prep()}\cindex{bake()}\cindex{dim()}
:::

The resulting engineered features\index{feature engineering} data frame\index{data frame} has `r format(nrow(cls_bake), big.mark = ",")` observations and `r format(ncol(cls_bake), big.mark = ",")` variables. That is a lot of features! Given the fact that for each writing sample, only a small subset of them will actually appear, most of our cells will be filled with zeros. This is what is known as a **sparse matrix**\index{sparse matrix}\index{matrix}.

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

When applying tokenization and feature engineering steps to text data the result is often contained in a matrix object. Using {recipes} a data frame with a matrix-like structure is returned. Remember, a matrix is a data frame where all the vector types\index{vector types} are the same.

Furthermore, the features are prefixed with the variable name and transformation step labels. In @exm-predict-class-recipe-prep we applied $tf$-$idf$ to the `text` variable. Therefore the features are prefixed with `tfidf_text_`.
:::

But we should pause. This is an unwieldy number of features, on for every single word, for a model and it is likely that many of these features are not useful for our classification task. Furthermore, the more features we have, the more chance these features will capture the nuances of these particular writing samples increasing the likelihood we overfit the model\index{overfitting}. All in all, we need to reduce the number of features.

We can filter out features by stopword\index{stopwords} list or by frequency\index{frequency} of occurrence. Let's start by frequency of occurrence. We can set the maximum number of the top features with an arbitrary threshold to start. The `step_tokenfilter()` function can filters out features on a number of criteria. Let's use the `max_tokens` argument to set the maximum number of features to 100.

This particular step needs to be applied before the `step_tfidf()` step, so we will add it to our recipe before the `step_tfidf()` step. This is seen in @exm-predict-class-recipe-tokenfilter.

::: {#exm-predict-class-recipe-tokenfilter}
```{r}
#| label: predict-class-recipe-tokenfilter
#| results: hold

# Rebuild recipe with tokenfilter step
cls_rec <-
  base_rec |>
  step_tokenize(text) |>
  step_tokenfilter(text, max_tokens = 100) |>
  step_tfidf(text, smooth_idf = FALSE)

# Prep and bake
cls_bake <-
  cls_rec |>
  prep() |>
  bake(new_data = NULL)

# Preview
dim(cls_bake)

cls_bake[1:5, 1:5]
```
\index{R packages!recipes}\index{R packages!textrecipes}
\cindex{step_tokenfilter()}\cindex{step_tokenize()}\cindex{step_tfidf()}
\cindex{prep()}\cindex{bake()}\cindex{dim()}
:::

We now have a manageable set of features, and fewer of which will have a as many zeros. Only during the interrogation step will we know if they are useful.

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

The `prep()` and `bake()` functions are useful for inspecting the recipe and the engineered features, but they are not required to build a recipe. When a recipe is added to a workflow, the `prep()` and `bake()` functions are called automatically as part of the process.
:::

<!-- Step 5: Interrogate -->

We are now ready to turn our attention to step 5 of our workflow, interrogating the data. In this step, we will first select a classification algorithm\index{classification algorithms}, then add this algorithm and our recipe to a workflow object. We will then use the workflow object to train and assess the resulting models, adjusting them until we believe we have a robust final model to apply on the testing set for our final evaluation.

<!-- Select a classification algorithm -->

There are many classification algorithms to choose from with their own strengths and shortcomings. In @tbl-predict-class-algorithms, we list some of the most common classification algorithms and their characteristics.

::: {#tbl-predict-class-algorithms tbl-colwidths="[15, 28, 28, 29]"}

| Algorithm | Strengths | Shortcomings | Tuning Recommendation |
|-----------|-----------|--------------|-----------------------|
| Logistic regression | Interpretable, fast, high-dimensional data | Linear relationship, not for complex tasks | Cross-validate regularization strength |
| Naive Bayes | Interpretable, fast, high-dimensional data, multi-class | Assumes feature (naive) independence, poor with small data | None |
| Decision trees | Nonlinear, interpretable, numerical/ categorical data | Overfitting, high variance | Cross-validate maximum tree depth |
| Random forest | Nonlinear, numerical/ categorical data, less overfitting | Less interpretable, poor with high-dimensional data | Cross-validate number of trees |
| Support vector machines | Nonlinear, high-dimensional data, numerical/ categorical | Requires parameter tuning, memory intensive | Cross-validate regularization parameter |
| Neural networks | Nonlinear, large data, auto feature learning | Overfitting, difficult to interpret, expensive | Cross-validate learning rate |

Common classification algorithms
:::

<!-- Models: simple to complex  -->

In the process of selecting an algorithm, simple, computationally efficient, and interpretable models are preferred over complex, computationally expensive, and uninterpretable models, all things being equal. Only if the performance of the simple model is not good enough should we move on to a more complex model.

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

{parsnip} provides a consistent interface to many different models, `r nrow(parsnip::model_db)` at the time of writing. You can peruse the list of models by running `parsnip::model_db`.

You can also retrieve the list of potential engines for a given model specification with the `show_engines()` function. For example, `show_engines("logistic_reg")` will return a data frame with the engines available for the logistic regression model specification. Note, the engines represent R packages that need to be installed to use the engine.
:::

With this end mind, we will start with a simple logistic regression\index{logistic regression} model to see how well we can classify the texts in the training set with the features we have engineered. We will use the `logistic_reg()` function from {parsnip} to specify the logistic regression model. We then select the implementation engine (`glmnet` General Linear Model) and the mode of the model (`classification`). The implementation engine is the software that will be used to fit the model. The code to set up the model specification is seen in @exm-predict-class-model-spec.

::: {#exm-predict-class-model-spec}
```{r}
#| label: predict-class-model-spec
#| output: false

# Create a model specification
cls_spec <-
  logistic_reg() |>
  set_engine("glmnet")

# Preview
cls_spec
```

```
Logistic Regression Model Specification (classification)

Computational engine: glmnet
```
\index{R packages!parsnip}\index{R packages!glmnet}
\cindex{logistic_reg()}\cindex{set_engine()}
:::

<!-- WARN: define hyperparameters here -->

Now, different algorithms will have different parameters that can be adjusted which can affect the performance of the model (see @tbl-predict-class-algorithms). As not to confuse these parameters with the features, which are also parameters of the model, these are given the name **hyperparameters**\index{hyperparameters}. The adjustment process is called **hyperparameter tuning**\index{hyperparameter tuning} and involves fitting the model to the training set with different hyperparameters and evaluating the model's performance to determine the best hyperparameter values to use for the model.

<!-- [ ] tmp fmt -->
\pagebreak

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

You can find the hyperparameters for a model-engine by consulting the `parsnip::model_db` object and unnesting the `parameters` column. For example, `parsnip::model_db |> filter(model == "logistic_reg") |> unnest(parameters)` will return a data frame with the hyperparameters for the logistic regression model.

To learn more about the hyperparameters for a specific model, you can consult the documentation for `parsnip` model (*e.g.* `?logistic_reg`).
:::

For example, the logistic regression model using `glmnet` can be tuned to prevent overfitting. The regularization typically applied is the LASSO (L1) penalty[^lasso]. The `logistic_reg()` function takes the arguments `penalty` and `mixture`. We set `mixture = 1`, but we now need to decide what value to use for the strength of the `penalty` argument. Values can range from 0 to 1, where 0 indicates no penalty and 1 indicates a maximum penalty.

[^lasso]: The LASSO (least absolute shrinkage and selection operator) is a type of regularization that penalizes the absolute value of the coefficients. In essence, it smooths the coefficients by shrinking them towards zero to avoid coefficients picking up on particularities of the training data that will not generalize to new data.

Instead of guessing, we will use {tune} to tune the hyperparameters of the model. The `tune()` function serves as a placeholder for the hyperparameters we want to tune. We can add the `tune()` function to our model specification to specify the hyperparameters we want to tune. The code is seen in @exm-predict-class-model-spec-tune.

::: {#exm-predict-class-model-spec-tune}
```{r}
#| label: predict-class-model-spec-tune-lasso
#| output: false

# Create a model specification (with tune)
cls_spec <-
  logistic_reg(penalty = tune(), mixture = 1) |>
  set_engine("glmnet")

# Preview
cls_spec
```
```
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet
```
:::

We can see now that the `cls_spec` model specification now includes the `tune()` function as the value for the `penalty` argument.

<!-- Combine recipe and model specification -->

To tune our model, we will need to combine our recipe and model specification into a workflow object which sequences our feature selection, engineering, and model selection. We will use the `workflow()` function from {workflows} to do this. The code is seen in @exm-predict-class-workflow.

::: {#exm-predict-class-workflow}
```{r}
#| label: predict-class-workflow
#| output: false

# Create a workflow
cls_wf <-
  workflow() |>
  add_recipe(cls_rec) |>
  add_model(cls_spec)

# Preview
cls_wf
```
\index{R packages!workflows}
\cindex{workflow()}\cindex{add_recipe()}\cindex{add_model()}

```
══ Workflow ═══════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ───────────────────────────────────
3 Recipe Steps

• step_tokenize()
• step_tokenfilter()
• step_tfidf()

── Model ──────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet
```
:::

We now have a workflow `cls_wf` that includes our recipe and model specification, including the `tune()` function as a placeholder for a range of values for the penalty hyperparameter. To tune the penalty hyperparameter\index{hyperparameters}, we use the `grid_regular()` function from {dials} to specify a grid of values to try. Let's choose a random set of 10 values, as seen in @exm-predict-class-model-spec-tune-grid-values.

::: {#exm-predict-class-model-spec-tune-grid-values}
```{r}
#| label: predict-class-model-spec-tune-grid-values

# Create a grid of values for the penalty hyperparameter
cls_grid <-
  grid_regular(penalty(), levels = 10)

# Preview
cls_grid
```
\index{R packages!dials}
\cindex{grid_regular()}
:::


The 10 values chosen to be in the grid range from nearly 0 to 1, where 0 indicates no penalty and 1 indicates a strong penalty.

Now to perform the tuning and arrive at an optimal value for `penalty` we need to create a tuning workflow. We do this by calling the `tune_grid()` function using our tuning model specification workflow, a resampling object, and our hyperparameter grid and return a `tune_grid` object.

Resampling\index{resampling} is a strategy that allows us to generate multiple training and testing sets from a single dataset ---in this case the training data we split at the outset. Each generated training-testing pair is called a fold. Which is why this type of resampling is called **k-fold cross-validation**\index{k-fold cross-validation}. The `vfold_cv()` function from {rsample} takes a data frame and a number of folds and returns a `vfold_cv` object. We will apply the `cls_wf` workflow to the 10 folds of the training set with `tune_grid()`. For each fold, each of the 10 values of the penalty hyperparameter will be tried and the model's performance will be evaluated. The code is seen in @exm-predict-class-model-spec-tune-grid-cv.

::: {#exm-predict-class-model-spec-tune-grid-cv}
```{r}
#| label: predict-class-model-spec-tune-grid-cv


set.seed(123) # for reproducibility

# Create a resampling object
cls_vfold <- vfold_cv(cls_train, v = 10)

# Tune the model
cls_tune <-
  tune_grid(
    cls_wf,
    resamples = cls_vfold,
    grid = cls_grid
    )

# Preview
cls_tune
```
\index{R packages!rsample}
\cindex{vfold_cv()}\cindex{tune_grid()}\cindex{set.seed()}
:::

The `cls_tune` object contains the results of the tuning for each fold. We can see the results of the tuning for each fold by calling the `collect_metrics()` function on the `cls_tune` object, as seen in @exm-predict-class-model-spec-tune-grid-collect. Passing the `cls_tune` object to `autoplot()` produces the visualization in @fig-predict-class-model-spec-tune-grid-collect.

::: {#exm-predict-class-model-spec-tune-grid-collect}

```r
# Collect the results of the tuning
cls_tune_metrics <-
  collect_metrics(cls_tune)

# Visualize metrics
autoplot(cls_tune)
```

```{r}
#| label: fig-predict-class-model-spec-tune-grid-collect
#| fig-cap: "Metrics for each fold of the tuning process"
#| fig-alt: "Two line plots one showing accuracy and the other showing ROC-AUC for each fold of the tuning process. The y-axis is the metric value and the x-axis is the penalty hyperparameter value."
#| fig-width: 6
#| fig-asp: 0.5
#| echo: false

# Collect the results of the tuning
cls_tune_metrics <-
  collect_metrics(cls_tune)

# Visualize metrics
cls_tune_metrics |>
  filter(.metric != "brier_class") |>
  ggplot(aes(x = penalty, y = mean)) + 
  geom_line() + 
  geom_point() +
  scale_x_log10() + 
  facet_wrap( ~ .metric, nrow = 2) + 
  labs(x = "Amount of regularization", y = "") +
  theme_qtalr(font_size = 10)
```
\index{R packages!tune}\index{R packages!ggplot2}
\cindex{collect_metrics()}\cindex{autoplot()}
:::

The most common metrics for model performance in classification are accuracy\index{accuracy} and the area under the **receiver operating characteristic area under the curve** (ROC-AUC)\index{receiver operating characteristic area under the curve (ROC-AUC)}. Accuracy is simply the proportion of correct predictions. The ROC-AUC provides a single score which summarizes how well the model can distinguish between classes. The closer to 1 the more discriminative power the model has.

In the plot of the metrics, @fig-predict-class-model-spec-tune-grid-collect, we can see that the many of the penalty values performed similarly, with a drop-off in performance at the higher values. Conveniently, the `show_best()` function from {tune} takes a `tune_grid` object and returns the best performing hyperparameter values\index{hyperparameters}. The code is seen in @exm-predict-class-model-spec-tune-grid-collect-best.

::: {#exm-predict-class-model-spec-tune-grid-collect-best}
```{r}
#| label: predict-class-model-spec-tune-grid-collect-best

# Show the best performing hyperparameter value
cls_tune |>
  show_best(metric = "roc_auc")
```
\index{R packages!tune}
\cindex{show_best()}
:::

We can make this selection programmatically by using the `select_best()` function. This function needs a metric to select by. We will use the ROC-AUC and select the best value for the penalty hyperparameter. The code is seen in @exm-predict-class-model-spec-tune-grid-collect-select.

\pagebreak

::: {#exm-predict-class-model-spec-tune-grid-collect-select}
```{r}
#| label: predict-class-model-spec-tune-grid-collect-select

# Select the best performing hyperparameter value
cls_best <-
  select_best(cls_tune, metric = "roc_auc")

# Preview
cls_best
```
\index{R packages!tune}
\cindex{select_best()}
:::

All of that to tune a hyperparameter! Now we can update the model specification and workflow with the best performing hyperparameter value using the previous `cls_wf_tune` workflow and the `finalize_workflow()` function. The `finalize_workflow()` function takes a workflow and the selected parameters and returns an updated `workflow` object, as seen in @exm-predict-class-tune-hyperparameters-update-workflow.

::: {#exm-predict-class-tune-hyperparameters-update-workflow}
```{r}
#| label: predict-class-tune-hyperparameters-update-workflow
#| output: false

# Update model specification
cls_wf_lasso <-
  cls_wf |>
  finalize_workflow(cls_best)

# Preview
cls_wf_lasso
```
\index{R packages!workflows}
\cindex{finalize_workflow()}

```
══ Workflow ═══════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ───────────────────────────────────

• step_tokenize()
• step_tokenfilter()
• step_tfidf()

── Model ──────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = 0.000464158883361278
  mixture = 1

Computational engine: glmnet
```
:::

Our model specification and the workflow are updated with the tuned hyperparameter.

As a reminder, we are still working in step 5 of our workflow, interrogating the data. So far, we have selected and engineered the features, split the data into training and testing sets, and selected a classification algorithm. We have also tuned the hyperparameters of the model and updated the model specification and workflow with the best performing hyperparameter value.

The next step is to assess the performance of the model on the training set given the features we have engineered, the algorithm we have selected, and the hyperparameters we have tuned. Instead of evaluating the model on the training set directly, we will use cross-validation\index{k-fold cross-validation} on the training set to gauge the variability of the model.

The reason for this is that the model's performance on the entire training set at once is not a reliable indicator of the model's performance on new data ---just imagine if you were to take the same test over and over again, you would get better and better at the test, but that doesn't mean you've learned the material any better. Cross-validation is a technique that allows us to estimate the model's performance on new data by simulating the process of training and testing the model on different subsets of the training data.

Similar to what we did to tune the hyperparameters, we can use cross-validation to gauge the variability of the model. The `fit_resamples()` function takes a workflow and a resampling object and returns metrics for each fold. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-workflow-cv.

::: {#exm-predict-class-tune-hyperparameters-evaluate-workflow-cv}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-workflow-cv

# Cross-validate workflow
cls_lasso_cv <-
  cls_wf_lasso |>
  fit_resamples(
    resamples = cls_vfold,
    # save predictions for confusion matrix
    control = control_resamples(save_pred = TRUE)
  )
```
\index{R packages!tune}
\cindex{fit_resamples()}
:::

We want to aggregate the metrics across the folds to get a sense of the variability of the model. The `collect_metrics()` function takes the results of a cross-validation and returns a data frame with the metrics.

\pagebreak

::: {#exm-predict-class-tune-hyperparameters-evaluate-workflow-cv-collect}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-workflow-cv-collect

# Collect metrics
collect_metrics(cls_lasso_cv)
```
\cindex{collect_metrics()}
:::

From the accuracy and ROC-AUC metrics in @exm-predict-class-tune-hyperparameters-evaluate-workflow-cv-collect it appears we have a decent candidate model, but there is room for potential improvement. A good next step is to evaluate the model errors and see if there are any patterns that can be addressed before considering what approach to take to improve the model.

::: {.callout .halfsize}
**{{< fa regular hand-point-up >}} Tip**

To provide context in terms of what is a good model performance, it is useful to compare the model's performance to a null model. A **null model**\index{null model} (or baseline model) is a simple model that is easy to implement and provides a benchmark for the model's performance. For classification tasks\index{text classification}, a common null model is to predict the most frequent class. In modeling, this is the minimal benchmark we want to beat, if we are doing better than this, we are doing better than chance.
:::

For classification tasks, a good place to start is to visualize a confusion matrix. A **confusion matrix**\index{confusion matrix} is a cross-tabulation of the predicted and actual outcomes. The `conf_mat_resampled()` function takes a `fit_resamples` object (with predictions saved) and returns a table (`tidy = FALSE`) with the confusion matrix for the aggregated folds. We can pass this to the `autoplot()` function to plot as in @exm-predict-class-tune-hyperparameters-evaluate-workflow-cv-confusion.

```{r}
#| label: fig-predict-class-tune-hyperparameters-evaluate-workflow-cv-confusion
#| include: false

# Plot confusion matrix
p <- cls_lasso_cv |>
  conf_mat_resampled(tidy = FALSE) |>
  autoplot(type = "heatmap") +
  theme_qtalr(font_size = 10) +
  theme(legend.position = "none")

ggsave("figures/fig-class-tune-hyperparameters-evaluate-workflow-cv-confusion.png", p, width = 4.5, height = 4, dpi = 300)
```

::: {#exm-predict-class-tune-hyperparameters-evaluate-workflow-cv-confusion}

```r
# Plot confusion matrix
cls_lasso_cv |>
  conf_mat_resampled(tidy = FALSE) |>
  autoplot(type = "heatmap")
```
\index{R packages!tune}\index{R packages!ggplot2}
\cindex{conf_mat_resampled()}\cindex{autoplot()}
:::

::: {#fig-class-tune-hyperparameters-evaluate-workflow-cv-confusion}

![](figures/fig-class-tune-hyperparameters-evaluate-workflow-cv-confusion.png){width="75%" fig-alt="Heatmap of the confusion matrix for the aggregated folds of the cross-validation. Actual classes are on the y-axis and predicted classes are on the x-axis."}

Confusion matrix for the aggregated folds of the cross-validation
:::

The top left to bottom right diagonal contains the true positives and true negatives. These are the correct predictions. The top right to bottom left diagonal contains the false positives and false negatives ---our errors. The convention is to speak of one class being the positive class and the other class being the negative class. In our case, we will consider the positive class to be the 'learner' class and the negative class to be the 'natives' class.

We can see that there are more learners falsely predicted to be natives than the other way around. This may be due to the fact that there are simply more learners than natives in the dataset or this could signal that there are some learners that are more similar to natives than other learners. Clearly this can't be the entire explanation as the model is not perfect, even some natives are classified falsely as learners! But it may be an interesting avenue for further exploration. Perhaps these are learners that are more advanced or have a particular style of writing that is more similar to natives.

::: {.callout .halfsize}
**{{< fa medal >}} Dive deeper**

Another perspective often applied to evaluate a model is the receiver operating characteristic (ROC) curve\index{receiver operating characteristic (ROC) curve}. The ROC curve is a plot of the true positive rate (TPR) against the false positive rate (FPR) for different classification thresholds. This metric, and visualization, can be useful to gauge the model's ability to distinguish between the two classes. {yardstick} provides the `roc_curve()` function to calculate the ROC curve on an `fit_resamples` object.
:::

<!-- Model improvement: tune max tokens -->

To improve supervised learning models, consider: 

1. Engineering the features differently
2. Selecting different (or additional) features
3. Changing the algorithm
4. Tuning the hyperparameters differently

Of these options, adjusting the feature engineering\index{feature engineering} process is the option that diverges least from our current workflow `cls_wf_lasso`. Recall that in our recipe specification we set a token filter to limit the number of features to 100. We can adjust this number to see if it has an effect on the model's performance.

To help select the optimal number of tokens\index{tokens}, we again can use the tuning process we explored for the hyperparameters. This time, however, the `tune()` placeholder will be included as the argument to the `max_tokens` argument in the `step_tokenfilter()` function.

I repeat the recipe with the tuning\index{hyperparameter tuning} placeholder in @exm-predict-class-tune-hyperparameters-tokenfilter.

::: {#exm-predict-class-tune-hyperparameters-tokenfilter}
```{r}
#| label: predict-class-tune-hyperparameters-tokenfilter

# Create a recipe with a token filter step
cls_rec <-
  recipe(
    formula = outcome ~ text,
    data = cls_train
    ) |>
  step_tokenize(text) |>
  step_tokenfilter(text, max_tokens = tune()) |>
  step_tfidf(text)
```
\index{R packages!recipes}\index{R packages!textrecipes}
\cindex{step_tokenfilter()}\cindex{step_tokenize()}\cindex{step_tfidf()}\cindex{\textasciitilde}
:::

With the updated recipe, we can update the `cls_wf_lasso` and tune the `max_tokens` hyperparameter. The code is seen in @exm-predict-class-tune-hyperparameters-update-rec.

::: {#exm-predict-class-tune-hyperparameters-update-rec}
```{r}
#| label: predict-class-tune-hyperparameters-update-rec

# Update workflow with token filter tuning
cls_wf_lasso <-
  cls_wf_lasso |>
  update_recipe(cls_rec)
```
\index{R packages!workflows}
\cindex{update_recipe()}
:::

One thing to note is that we will want to consider what values of `max_tokens` we want to use to tune the hyperparameter. So instead of only specifying the levels in the `grid_regular()` function, we are best off to provide a range of values that we think are reasonable. Let's add a range of values between our current value 100 and 2,000 to start. And let's tell the grid to select five values from this range.

The code is seen in @exm-predict-class-tune-hyperparameters-tokenfilter-grid.

\pagebreak

::: {#exm-predict-class-tune-hyperparameters-tokenfilter-grid}
```{r}
#| label: predict-class-tune-hyperparameters-tokenfilter-grid

# Create a grid of values for the max tokens hyperparameter
cls_grid <-
  grid_regular(max_tokens(range = c(100, 2000)), levels = 5)

# Preview
cls_grid
```
\index{R packages!dials}
\cindex{grid_regular()}
:::

From here, the process is the same as before. We will use the `tune_grid()` function to tune the `max_tokens` hyperparameter, select the best value, and finalize the workflow, as seen from @exm-predict-class-model-spec-tune-grid-cv through @exm-predict-class-tune-hyperparameters-update-workflow.

```{r}
#| label: predict-class-tune-hyperparameters-tokenfilter-tune
#| include: false
#| cache: true


set.seed(123) # for reproducibility

# Tune the model
cls_tune <-
  tune_grid(
    cls_wf_lasso,
    resamples = cls_vfold,
    grid = cls_grid
    )

# Select the best performing hyperparameter value
cls_best <-
  select_best(cls_tune, metric = "roc_auc")

# Update model specification
cls_wf_lasso_tokens <-
  cls_wf_lasso |>
  finalize_workflow(cls_best)
```

After tuning the `max_tokens` hyperparameter\index{hyperparameters}, the best performing value is 1,050. We now used the updated `cls_wf_lasso_tokens` workflow to cross-validate\index{k-fold cross-validation} the model and collect the metrics. The code is seen in @exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-collect.

::: {#exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-collect}
```{r}
#| label: predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-collect

# Cross-validate workflow
cls_lasso_tokens_cv <-
  cls_wf_lasso_tokens |>
  fit_resamples(
    resamples = cls_vfold,
    # save predictions for confusion matrix
    control = control_resamples(save_pred = TRUE)
  )

# Collect metrics
collect_metrics(cls_lasso_tokens_cv)
```
\index{R packages!tune}
\cindex{fit_resamples()}\cindex{collect_metrics()}
:::

The metrics from @exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-collect show that the model's performance has improved for both the accuracy\index{accuracy} and the ROC-AUC\index{receiver operating characteristic area under the curve (ROC-AUC)}. The confusion matrix from @exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion shows that the number of false positives and false negatives has decreased. This is a good sign that the model is more robust.

```{r}
#| label: fig-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion
#| include: false

# Plot confusion matrix 
p <- cls_lasso_tokens_cv |>
  conf_mat_resampled(tidy = FALSE) |>
  autoplot(type = "heatmap") +
  theme_qtalr(font_size = 10) +
  theme(legend.position = "none")

ggsave("figures/fig-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion.png", p, width = 4.5, height = 4, dpi = 300)
```
<!--TODO: Issue with the sizing of the heatmap(s)-->
::: {#exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion}
```r
# Plot confusion matrix
cls_lasso_tokens_cv |>
  conf_mat_resampled(tidy = FALSE) |>
  autoplot(type = "heatmap")
```
\index{R packages!tune}\index{R packages!ggplot2}
\cindex{conf_mat_resampled()}\cindex{autoplot()}
:::

::: {#fig-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion}

![](figures/fig-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion.png){width="75%" fig-alt="Heatmap of the confusion matrix for the aggregated folds of the cross-validation. Actual classes are on the y-axis and predicted classes are on the x-axis."}

Confusion matrix for the aggregated folds of the cross-validation
:::

From @fig-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-confusion, it appears that the model is more robust with the updated `max_tokens` hyperparameter. We could continue to explore other model improvement strategies, but for now we will move on to the next step in our workflow.

<!-- Step 6: Evaluate -->

We are now ready to move on to step 7, evaluating the model on the test set\index{model evaluation}. To do this we need to fit the tuned workflow to the training set, which is the actual training phase. We will use the `last_fit()` function from {workflows} to fit the workflow to the training set.

The `last_fit()` function takes a workflow and a split object and returns a `last_fit` object. The `last_fit` object contains the results of the model fit on the training set and the results of the model evaluation on the test set. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test.

We will use the `last_fit()` function to train the final model and predict the outcome on the test set. The `collect_metrics()` function takes a data frame with the actual and predicted outcomes and returns a data frame with the metrics for the model. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test.

::: {#exm-predict-class-tune-hyperparameters-evaluate-test}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-test

# Fit the model to the training set and evaluate on the test set
cls_final_fit <-
  last_fit(
    cls_wf_lasso_tokens,
    split = cls_split
  )

# Evaluate model on testing set
collect_metrics(cls_final_fit)
```
\index{R packages!workflows}\index{R packages!tune}
\cindex{last_fit()}\cindex{collect_metrics()}
:::

The performance metrics are very close to those we achieved on the training set in @exm-predict-class-tune-hyperparameters-tokenfilter-evaluate-workflow-cv-collect. This is a good sign that the model is robust as it performs well on both training and test sets. We can evaluate the confusion matrix\index{confusion matrix} on the test set as well. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test-confusion and the visualization in @fig-class-tune-hyperparameters-evaluate-workflow-cv-confusion.

```{r}
#| label: fig-predict-class-tune-hyperparameters-evaluate-test-confusion
#| fig-width: 6
#| fig-asp: 0.5
#| include: false

# Plot confusion matrix
p <- cls_final_fit |>
  collect_predictions() |>
  conf_mat(truth = outcome, estimate = .pred_class) |>
  autoplot(type = "heatmap") +
  theme_qtalr(font_size = 10) +
  theme(legend.position = "none")

ggsave("figures/fig-predict-class-tune-hyperparameters-evaluate-test-confusion.png", p, width = 4.5, height = 4, dpi = 300)
```

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-confusion}

```r
# Plot confusion matrix
cls_final_fit |>
  collect_predictions() |>
  conf_mat(truth = outcome, estimate = .pred_class) |>
  autoplot(type = "heatmap")
```
\index{R packages!workflows}\index{R packages!tune}\index{R packages!yardstick}\index{R packages!ggplot2}
\cindex{collect_predictions()}\cindex{conf_mat()}\cindex{autoplot()}
:::

::: {#fig-class-tune-hyperparameters-evaluate-workflow-cv-confusion}
![](figures/fig-predict-class-tune-hyperparameters-evaluate-test-confusion.png){width="75%" fig-alt="Heatmap of the confusion matrix for the test set. Actual classes are on the y-axis and predicted classes are on the x-axis."}

Confusion matrix for the test set
:::

On the test set the false instances are balanced, which is a good sign that the model is robust. Ideally, there would be no errors, but this is not realistic. The model is not perfect, but it is useful.

Now a model that can predict the nativeness of a writer based on their writing sample is a useful tool in itself. You could imagine that this could be a pre-processing step for a language learning application, for example. But for a study that is more interested in learning about what features are most important for predicting the native versus non-native features of a writer, we still have some work to do\index{feature importance}. We can inspect the errors on the test set to gain some insight into what writing samples, and which proficiency levels of the writers, are most difficult to predict. We can also inspect the estimates for the features in the model to gain some insight into what features are most important for predicting the outcomes.

Let's first approach this from a document-proficiency point of view. First, we will want to integrate the predictions with the test set to inspect the errors. We can use the `collect_predictions()` function to collect the predictions from the `last_fit` object and attach them with the test set `cls_test` with `bind_cols`. Note, we can drop the `outcome` variable from `cls_test` as we have this column in our fitted model. The code is seen in @exm-predict-class-tune-hyperparameters-integrate-test.

::: {#exm-predict-class-tune-hyperparameters-integrate-test}
```{r}
#| label: predict-class-tune-hyperparameters-integrate-test

# Collect predictions from the model
cls_lasso_fit_preds_test <-
  cls_final_fit |>
  collect_predictions() |>
  bind_cols(cls_test[, -1])

# Preview
glimpse(cls_lasso_fit_preds_test)
```
\index{R packages!dplyr}\index{R packages!tune}
\cindex{collect_predictions()}\cindex{bind_cols()}\cindex{glimpse()}
:::

I will then select the columns with the actual outcome, the predicted outcome, the proficiency level, and the text and separate the predicted outcome to inspect them separately, as seen in @exm-predict-class-tune-hyperparameters-evaluate-test-errors.

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-errors}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-test-errors

# Inspect errors
cls_lasso_fit_preds_test |>
  filter(outcome != .pred_class) |>
  select(outcome, .pred_class, proficiency, text)

# Inspect learners falsely predicted to be natives
cls_lasso_fit_preds_test |>
  filter(outcome == "Learner", .pred_class == "Native") |>
  select(outcome, .pred_class, proficiency, text) |>
  count(proficiency, sort = TRUE)
```
\index{R packages!dplyr}
\cindex{filter()}\cindex{select()}\cindex{count()}
:::

Interestingly, the majority of misclassified learners are advanced, which could be expected as they are more similar to natives. There are some beginners that are misclassified as natives, but this is not as common. Yes, it is still an open question as to why some natives are classified as learners.

We can inspect the estimates for the features\index{predictor features} in the model to gain some insight into what features are most important for predicting the outcomes. The `extract_fit_parsnip()` function takes a trained model specification `cls_final_fit` and returns a data frame with the estimated coefficients for each feature. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test-estimates.

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-estimates}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-test-estimates

# Extract estimates
cls_final_fit_features <-
  cls_final_fit |>
  extract_fit_parsnip() |>
  tidy()
```
\index{R packages!hardhat}
\cindex{extract_fit_parsnip()}\cindex{tidy()}
:::

The estimates are the log odds of the outcome. In a binary classification task, the log odds of the outcome is the log of the probability of the outcome divided by the probability of the other outcome. In our case, the reference outcome is "Learner", so negative log-odds indicate that the feature is associated with the "Learner" outcome and positive log-odds indicate that the feature is associated with the "Native" outcome.

The estimates are in log-odds, so we need to exponentiate them to get the odds. The odds are the probability of the outcome divided by the probability of the other outcome. The probability of the outcome is the odds divided by the odds plus one. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test-estimates-probability.

\pagebreak

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-estimates-probability}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-test-estimates-probability

# Calculate probability
cls_final_fit_features |>
  mutate(probability = exp(estimate) / (exp(estimate) + 1))
```
\index{R packages!dplyr}
\cindex{mutate()}\cindex{exp()}
:::

So just looking at the snippet of the features returned from @exm-predict-class-tune-hyperparameters-evaluate-test-estimates-probability, we can see that the features 'a' and 'abandonado' are associated with the "Native" outcome, 'abuela' is associated with "Learners", and the other features are neutral (`probability` = 0.5).

A quick way to extract the most important features for predicting each outcome is to use the `vi()` function from {vip}\index{feature importance}. It takes a trained model specification and returns a data frame with the most important features. The code is seen in @exm-predict-class-tune-hyperparameters-evaluate-test-vip.

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-vip}
```{r}
#| label: predict-class-tune-hyperparameters-evaluate-test-vip

# Load package
library(vip)

# Avoid conflicts for function names from other packages
conflicted::conflicts_prefer(vip::vi)

# Extract important features
var_importance_tbl <-
  cls_final_fit |>
  extract_fit_parsnip() |>
  vi()

# Preview
var_importance_tbl
```
\index{R packages!vip}\index{R packages!hardhat}
\cindex{vi()}\cindex{extract_fit_parsnip()}\cindex{library()}
:::

The `Variable` column contains each feature (with the feature type and corresponding variable `tfidf_text_`), `Importance` provides the absolute log-odds value, and the `Sign` column indicates whether the feature is associated with the "NEG" ("Learner") or the "POS" ("Native") outcome. We can recode the `Variable` and `Sign` columns to make them more interpretable and then plot them using `ggplot()`, as in @exm-predict-class-tune-hyperparameters-evaluate-test-vip-plot.

::: {#exm-predict-class-tune-hyperparameters-evaluate-test-vip-plot}
```{r}
#| label: fig-predict-class-tune-hyperparameters-evaluate-test-vip-plot
#| fig-cap: "Most important features for predicting the outcome"
#| fig-alt: "scatterplot of the most important features for predicting the outcome. The y-axis is the importance value and the x-axis is the feature. The features are faceted by the outcome."
#| fig-width: 6
#| fig-asp: 0.65

# Recode variable and sign
var_importance_tbl <-
  var_importance_tbl |>
  mutate(
    Feature = str_remove(Variable, "tfidf_text_"),
    Outcome = case_when(
      Sign == "NEG" ~ "Learner",
      Sign == "POS" ~ "Native"),
    ) |>
  select(Outcome, Feature, Importance)

# Plot
var_importance_tbl |>
  slice_max(Importance, n = 50) |>
  ggplot(aes(x = reorder(Feature, Importance), y = Importance)) +
  geom_point() +
  coord_flip() +
  facet_wrap(~ Outcome, scales = "free_y") +
  labs(x = NULL, y = "Importance", fill = NULL) +
  theme_minimal(base_size = 10)
```
\index{R packages!dplyr}\index{R packages!ggplot2}
\cindex{mutate()}\cindex{select()}\cindex{slice_max()}\cindex{geom_point()}\cindex{coord_flip()}\cindex{facet_wrap()}\cindex{labs()}\cindex{theme_minimal()}\cindex{reorder()}
:::

We can inspect @fig-predict-class-tune-hyperparameters-evaluate-test-vip-plot, and qualitatively assess what these features may be telling us about the differences between the learners and the natives\index{research interpretation}.

In this section, we've build a text classifier using a regularized logistic regression model. We've tuned the hyperparameters to arrive at a robust model that performs well on both the training and test sets. We've also evaluated the model errors and inspected the most important features for predicting the outcome.

### Text regression {#sec-predict-text-regression}

We will now turn our attention to the second task in this section, text regression\index{text regression}. In this task, we will use the same original dataset as in the classification task, but we will predict the placement score based on the learner writing samples. I will make reference to but not repeat the steps we took in the classification task, as many of the steps are the same. This is one of the benefits of using Tidymodels\index{Tidymodels}---the workflow is by-and-large the same for different tasks.

Let's start by extracting the observations (only learners) and the relevant variables from the original dataset. The code is seen in @exm-predict-reg-data.

::: {#exm-predict-reg-data}
```{r}
#| label: predict-reg-data

# Extract observations and relevant variables
reg_tbl <-
  cedel_tbl |>
  filter(subcorpus == "Learner") |>
  select(outcome = place_score, proficiency, text)

# Preview
glimpse(reg_tbl)
```
\index{R packages!dplyr}
\cindex{filter()}\cindex{select()}\cindex{glimpse()}
:::

<!-- Step 1: Identify variables -->

In this task, our outcome variable\index{outcome variable} is numeric and our predictor variable `text` is the same as before\index{predictor variables}. It might be useful to engineer the features differently, but we will start with the same feature engineering process as before, namely the $tf$-$idf$ method for the top 1,050 words\index{feature engineering}.

<!-- Step 2: Initial split -->
<!-- Step 3: Create recipe -->
<!-- Step 4: Verify recipe -->

I create a data split, `reg_split` and the training and testing sets, `reg_train` and `reg_test` and create a `reg_rec` recipe object which contains the starting recipe for the regression task. And, since we are using the same recipe as before, there is no need to validate the recipe. We can skip straight to the model building.

```{r}
#| label: predict-reg-initial-setup
#| include: false

set.seed(123) # for reproducibility

# Split data
reg_split <-
  initial_split(reg_tbl, prop = 0.8)

# Training set
reg_train <-
  training(reg_split)

# Test set
reg_test <-
  testing(reg_split)

# Create a recipe
reg_rec <-
  recipe(outcome ~ text, data = reg_train) |>
  step_tokenize(text) |>
  step_tokenfilter(text, max_tokens = 1050) |>
  step_tfidf(text, smooth_idf = FALSE)
```

<!-- Step 5: Interrogate -->

As before, we will want to start with a simple model and then build up to more complex models. The list in @tbl-predict-class-algorithms, includes algorithms that are commonly used in classification tasks. Interestingly, many of these same algorithms can be applied to regression. One exception is that instead of logistic regression, linear regression is used for numeric outcomes\index{regression algorithms}\index{machine learning}\index{supervised learning}. As with logistic regression, linear regression model\index{linear regression} is one of the simpler models. And just as with logistic regression, we will want to tune the regularization hyperparameter\index{hyperparameter tuning} of the linear regression model. Instead of detailing these steps again, let me summarize the process, in @tbl-workflow, and then we will discuss the results from the regularized linear regression model.

::: {#tbl-workflow tbl-colwidths="[5, 95]"}

| Step | Description |
| --- | --- |
| 1 | Build a model specification with a placeholder to tune the model. |
| 2 | Create a workflow with the recipe and the model specification. |
| 3 | Create a grid of values for the regularization hyperparameter. |
| 4 | Tune the model using cross-validation. |
| 5 | Select the best performing hyperparameter value (based on RMSE). |
| 6 | Update the model specification and workflow with the best performing hyperparameter value. |
| 7 | Fit the model to the training set and evaluate the performance using cross-validation. |

Steps to build and tune a model
:::

```{r}
#| label: predict-reg-linear-model-lasso
#| include: false
#| cache: true

# Create model specification
reg_spec <-
  linear_reg(penalty = tune(), mixture = 1) |>
  set_engine("glmnet")

# Create workflow
reg_wf <-
  workflow() |>
  add_recipe(reg_rec) |>
  add_model(reg_spec)

# Create tuning grid
reg_grid <- grid_regular(penalty(), levels = 10)

# Create the cross-validation folds
reg_vfold <- vfold_cv(reg_train, v = 10)


set.seed(123) # for reproducibility

# Create tuning workflow
reg_tune <-
  tune_grid(
    reg_wf,
    resamples = reg_vfold,
    grid = reg_grid
  )

# Select best parameter
chosen_penalty <- reg_tune |> select_best(metric = "rmse")

# Update workflow
reg_wf_lasso <-
  reg_wf |>
  finalize_workflow(chosen_penalty)

# Cross-validated workflow
reg_lasso_cv <-
  reg_wf_lasso |>
  fit_resamples(
    resamples = reg_vfold,
    control = control_resamples(save_pred = TRUE)
  )
```

Applying the steps 1 through 7, we have cross-validated\index{k-fold cross-validation} results for our model in the `reg_lasso_cv` object. We can collect the relevant metrics, root mean squared error (RMSE)\index{root mean squared error (RMSE)} and R-squared ($R^2$) values\index{R-squared}. Let's aggregate these measures using the code is seen in @exm-predict-reg-metrics-lasso.

::: {#exm-predict-reg-metrics-lasso}
```{r}
#| label: predict-reg-metrics-lasso

# Collect metrics
collect_metrics(reg_lasso_cv)
```
\index{R packages!tune}
\cindex{collect_metrics()}
:::

```{r}
#| label: predict-reg-metrics-lasso-results
#| echo: false

# Collect metrics
estimates <-
  reg_lasso_cv |>
  collect_metrics() |>
  pull(3)
```

Now, the **root mean squared error** (RMSE)\index{root mean squared error (RMSE)} estimate is `r format(estimates[1])`. RMSE is a measure of the difference between the predicted and the actual values expressed in the same units as the outcome variable. In this case, the outcome variable is the placement test score percent. So the RMSE is `r format(estimates[1])` percentage points. **R-squared** ($R^2$)\index{R-squared} is a measure of the proportion of the variance in the outcome variable that is explained by the model. The $R^2$ estimate is `r format(estimates[2])`. This means that the model explains `r scales::percent(estimates[2])` of the variance in the outcome variable. Taken together, this isn't the greatest model.

But how good or bad is it? This is where we can use the null model to compare the model to\index{null model}. The null model is a model that predicts the mean of the outcome variable for each of the outcomes. We can use the `null_model()` function to create a null model and submit it to cross-validation, @exm-predict-reg-model-spec-tune-fit-evaluate-null.

::: {#exm-predict-reg-model-spec-tune-fit-evaluate-null}
```{r}
#| label: predict-reg-model-spec-tune-fit-evaluate-null

# Create null model
null_model <-
  null_model() |>
  set_engine("parsnip") |>
  set_mode("regression")

# Cross-validate null model
null_cv <-
  workflow() |>
  add_recipe(reg_rec) |>
  add_model(null_model) |>
  fit_resamples(
    resamples = vfold_cv(reg_train, v = 10),
    metrics = metric_set(rmse)
  )

# Collect metrics
collect_metrics(null_cv)
```
\index{R packages!parsnip}\index{R packages!workflows}\index{R packages!tune}
\cindex{null_model()}\cindex{set_engine()}\cindex{set_mode()}\cindex{fit_resamples()}\cindex{collect_metrics()}
::: 

::: {.callout .halfsize}
**{{< fa exclamation-triangle >}} Warning**

For model specifications in which the model can be used in a classification or regression task, the model specification must be set to the correct mode before fitting the model. We have not set the mode for the `logistic_reg()` or `linear_reg()` model specifications, as the task is inferred. However, we have set the mode for the `null_model()`, and other model specifications that can be used in both classification and regression tasks.
:::

```{r}
#| label: predict-reg-metrics-null-text-block-1
#| echo: false

null_estimate <- null_cv |> collect_metrics() |> pull(3)
```

Our regression model performs better than the null model (`r round(null_estimate, 1)`), which means that it is picking up on some signal in the data.

Let's visualize the distribution of the predictions and the errors from our model to see if there are any patterns of interest. We can use the `collect_predictions()` function to extract the predictions of the cross-validation and plot the true outcome against the predicted outcome using `ggplot()`, as in @exm-predict-reg-lr-eval-rmse.

<!--TODO: add theme_qtalr-->

::: {#exm-predict-reg-lr-eval-rmse}
```{r}
#| label: fig-predict-reg-lr-eval-rmse
#| fig-cap: "Distribution of the RMSE for the cross-validated linear regression model"
#| fig-alt: "scatterplot of the true outcome against the predicted outcome for the cross-validated linear regression model. The x-axis is the true outcome and the y-axis is the predicted outcome. The shape of the points indicate the fold."
#| fig-width: 6
#| fig-asp: 0.5

# Visualize predictions
reg_lasso_cv |>
  collect_predictions() |>
  ggplot(aes(outcome, .pred, shape = id)) +
  geom_point(alpha = 0.5, position = position_jitter(width = 0.5)) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.5) + # trend for each fold
  labs(
    x = "Truth",
    y = "Predicted score",
    shape = "Fold"
  )
```
\index{R packages!ggplot2}\index{R packages!tune}
\cindex{collect_predictions()}\cindex{geom_point()}\cindex{geom_smooth()}\cindex{labs()}\cindex{position_jitter()}\cindex{ggplot()}\cindex{aes()}
:::

From @fig-predict-reg-lr-eval-rmse, we see data points for each predicted and truth value pair for each of the ten folds. There is a trend line for each fold which shows the linear relationship between the predicted and truth values for each fold. The trend lines are more similar than different, which is a good sign that the model is not wildly overfitting\index{overfitting} the training data. Looking closer, however, we can see the errors. Some are noticeably distant from the linear trend lines, *i.e.* outliers, in particular for test scores in the lower ranges\index{scatterplot}.

If the $R^2$\index{R-squared} value is in the ballpark, this means that somewhere around 40% of the variation is not explained by the frequency of the top 1,050 words. This is not surprising, as there are many other factors that contribute to the proficiency level of a text.

We have a model that is performing better than the null model\index{null model}, but it is not performing well enough to be very useful. We will need to update the model specification and/ or the features to try to improve the model fit. Let's start with the model. There are many different model specifications we could try, but we will likely need to use a more complex model specification to capture the complexity that we observe in the errors from the current linear regression model.

Let's try a decision tree model. **Decision trees**\index{decision trees} are models that are able to model nonlinear relationships and interactions between the features and the outcome and tend to be less influenced by outliers. Furthermore, decision trees are interpretable, which is a nice feature for an exploratory-oriented analysis\index{feature importance}. These are all desirable characteristics. Decision trees, however, can be prone to overfitting. For this reason, we will tune the maximum depth of the tree\index{hyperparameters} to minimize overfitting\index{overfitting}.

To implement a new model in Tidymodels, we need to create a new model specification and a new workflow. We will use the `decision_tree()` function from {parsnip} to create the model specification. The `decision_tree()` function takes a `tree_depth` argument that we want to tune. We create the new model specification with the tuning placeholder in @exm-predict-reg-model-spec-decision-tree.

::: {#exm-predict-reg-model-spec-decision-tree}
```{r}
#| label: predict-reg-model-spec-decision-tree
#| output: false

# Create model specification
reg_spec <-
  decision_tree(tree_depth = tune()) |>
  set_engine("rpart") |>
  set_mode("regression")

# Preview
reg_spec
```
\index{R packages!parsnip}\index{R packages!rpart}
\cindex{decision_tree()}\cindex{set_engine()}\cindex{set_mode()}\cindex{tune()}
```
Decision Tree Model Specification (regression)

Main Arguments:
  tree_depth = tune()

Computational engine: rpart
```
:::

With the model and tuning specification in place, we can now continue through the steps outlined in @tbl-workflow for this decision tree model. To create the grid of values for the tree depth hyperparameter, we will include the `grid_regular()` function with 10 levels, as seen in @exm-predict-reg-decision-tree-depth.

::: {#exm-predict-reg-decision-tree-depth}
Tuning values for the tree depth hyperparameter

```r
reg_grid <-
  grid_regular(tree_depth(), levels = 10)
```
\index{R packages!tune}
\cindex{grid_regular()}\cindex{tree_depth()}
:::

```{r}
#| label: predict-reg-decision-tree-depth
#| cache: true
#| include: false

# Create a workflow
reg_wf <-
  workflow() |>
  add_recipe(reg_rec) |>
  add_model(reg_spec)

# Create a grid of values for the tree depth hyperparameter
reg_grid <-
  grid_regular(tree_depth(), levels = 10)

# Create the cross-validation folds
reg_vfold <- vfold_cv(reg_train, v = 10)


set.seed(123) # for reproducibility

# Create tuning workflow
reg_tune <-
  tune_grid(
    reg_wf,
    resamples = reg_vfold,
    grid = reg_grid
  )

# Select best parameter
chosen_depth <- reg_tune |> select_best(metric = "rmse")

# Update workflow
reg_wf_tree <-
  reg_wf |>
  finalize_workflow(chosen_depth)

# Cross-validated workflow
reg_tree_cv <-
  reg_wf_tree |>
  fit_resamples(
    resamples = reg_vfold,
    control = control_resamples(save_pred = TRUE)
  )
```

We can collect the metrics and inspect the RMSE\index{root mean squared error (RMSE)} and $R^2$\index{R-squared} values. The code is seen in @exm-predict-reg-metrics-tree.

::: {#exm-predict-reg-metrics-tree}
```{r}
#| label: predict-reg-metrics-tree

# Collect metrics
collect_metrics(reg_tree_cv)
```
\index{R packages!tune}
:::

```{r}
#| label: predict-reg-metrics-tree-results
#| echo: false

# Collect metrics
estimates <-
  reg_tree_cv |>
  collect_metrics() |>
  pull(3)
```

The performance for the decision tree is worse than the regularized linear regression model. The RMSE is `r format(estimates[1])` and the $R^2$ is `r format(estimates[2])`. And, if we compare the standard error between the two models, we can see that the decision tree model has a lower standard error. This means that the decision tree model is likely overfitting, despite our efforts to tune tree depth.

Given the sensitivity of the decision tree branching process and random initialization, it is possible that the decision tree model is capturing too much nuance, and not enough generalities. Re-running the model with a different seed may result in a different model. This is a limitation with decision tree models, but it is also a feature, if we consider combining multiple decision trees to make a prediction. This is the basis of **ensemble models**\index{ensemble models}. An ensemble model is a model that combines multiple models with the goal to draw out the strengths of each model and minimize the weaknesses.

A **random forest**\index{random forest} is an ensemble model that combines multiple decision trees to make a prediction. In addition, random forests also perform random feature selection. This helps to reduce the correlation between the decision trees and thus works to reduce the overall variance of the model.

Let's try a random forest model to address our text regression task. We will use the `rand_forest()` function from {parsnip} to create the model specification. The `rand_forest()` function also takes a hyperparameter for the number of trees to be used in the model. We will select the `ranger` engine. Additionally, we will add the `importance` argument to ensure that we can extract feature importance if this model proves to be useful. We create the new model specification in @exm-predict-reg-model-spec-random-forest.

::: {#exm-predict-reg-model-spec-random-forest}
```{r}
#| label: predict-reg-model-spec-random-forest
#| output: false

# Create model specification
reg_spec <-
  rand_forest(trees = tune()) |>
  set_engine("ranger", importance = "impurity") |>
  set_mode("regression")

# Preview
reg_spec
```
```
Random Forest Model Specification (regression)

Main Arguments:
  trees = tune()

Engine-Specific Arguments:
  importance = impurity

Computational engine: ranger
```
\index{R packages!parsnip}\index{R packages!ranger}
\cindex{rand_forest()}\cindex{set_engine()}\cindex{set_mode()}\cindex{tune()}
:::

::: {.callout .halfsize}
**{{< fa regular lightbulb >}} Consider this**

The model building process is iterative and many of the steps are the same. This is a good indication that creating a custom function to build and tune the model would be a good idea.

Consider the following: What would you include in the function? What would you leave out? What required and/ or optional arguments would you include? What would you hard code? What would you return?
:::

Again, we apply the steps in @tbl-workflow to build and tune the random forest model. As part of this process, I will limit the range of the number of trees from 100 to 500 in five levels in the tuning grid, as seen in @exm-predict-reg-random-forest-trees.

::: {#exm-predict-reg-random-forest-trees}
Tuning values for the number of trees hyperparameter

```r
reg_grid <-
  grid_regular(trees(range = c(100, 500)), levels = 5)
```
\cindex{grid_regular()}\cindex{trees()}
:::

```{r}
#| label: predict-reg-random-forest-trees
#| cache: true
#| include: false

# Create a workflow
reg_wf <-
  workflow() |>
  add_recipe(reg_rec) |>
  add_model(reg_spec)

# Create a grid of values for the number of trees hyperparameter
reg_grid <-
  grid_regular(trees(range = c(250, 750)), levels = 5)

# Create the cross-validation folds
reg_vfold <- vfold_cv(reg_train, v = 10)


set.seed(123) # for reproducibility

# Create tuning workflow
reg_tune <-
  tune_grid(
    reg_wf,
    resamples = reg_vfold,
    grid = reg_grid
  )

# Select best parameter
chosen_trees <- reg_tune |> select_best(metric = "rmse")

# Update workflow
reg_wf_rf <-
  reg_wf |>
  finalize_workflow(chosen_trees)

# Cross-validated workflow
reg_rf_cv <-
  reg_wf_rf |>
  fit_resamples(
    resamples = reg_vfold,
    control = control_resamples(save_pred = TRUE)
  )
```

Let's collect the metrics and inspect the RMSE and $R^2$ values. The code is seen in @exm-predict-reg-metrics-rf.

::: {#exm-predict-reg-metrics-rf}
```{r}
#| label: predict-reg-metrics-rf

# Collect metrics
collect_metrics(reg_rf_cv)
```
\index{R packages!tune}
\cindex{collect_metrics()}
:::

```{r}
#| label: predict-reg-metrics-rf-results
#| echo: false

# Collect metrics
estimates <-
  reg_rf_cv |>
  collect_metrics() |>
  pull(3)
```

The random forest model performs better than the decision tree model and the regularized linear regression model. The RMSE is `r format(estimates[1])` and the $R^2$ is `r format(estimates[2])`. We also see that the standard error falls between the models we have tried so far.

Before we settle on this model, let's try one more model. In this case, we will introduce a neural network model. **Neural networks**\index{neural networks} are models that are able to model nonlinear relationships and interactions between the features and the outcome. They are also able to model complex relationships between the features and the outcome. We will use the `mlp()` function from {parsnip} to create the model specification. We will choose the `brulee` engine which allows us to tune the learning rate. The learning rate is a hyperparameter that controls the size of the steps that the model takes to update the weights during training.

::: {.callout .halfsize}
**{{< fa exclamation-triangle >}} Warning**

The `brulee` engine requires that the Torch computing resources are available on the computing environment. To facilitate the installation, {torch} provides the `install_torch()` function. This is a one-time operation. From this point on, R packages which depend on Torch will be able to take advantage of this rich machine learning framework.

:::

 We create the new model specification with the tuning placeholder in @exm-predict-reg-model-spec-mlp.

::: {#exm-predict-reg-model-spec-mlp}
```{r}
#| label: predict-reg-model-spec-mlp
#| output: false

# Create model specification
reg_spec <-
  mlp(learn_rate = tune()) |>
  set_engine("brulee") |>
  set_mode("regression")

# Preview
reg_spec
```
\index{R packages!parsnip}\index{R packages!brulee}
\cindex{mlp()}\cindex{set_engine()}\cindex{set_mode()}\cindex{tune()}

```
Single Layer Neural Network Model Specification (regression)

Main Arguments:
  learn_rate = tune()

Computational engine: brulee

Model fit template:
brulee::brulee_mlp(x = missing_arg(), y = missing_arg(), learn_rate = tune())
```
:::

And include the code in @exm-predict-reg-mlp-learning-rate to create a grid of values for the learning rate hyperparameter, as part of the model building workflow.

::: {#exm-predict-reg-mlp-learning-rate}
Tuning values for the learning rate hyperparameter

```r
reg_grid <-
  grid_regular(learn_rate(), levels = 10)
```
\cindex{grid_regular()}\cindex{learn_rate()}
:::

```{r}
#| label: predict-reg-mlp-learn-rate
#| cache: true
#| include: false

# Create a workflow
reg_wf <-
  workflow() |>
  add_recipe(reg_rec) |>
  add_model(reg_spec)

# Create a grid of values for the learning rate hyperparameter
reg_grid <-
  grid_regular(learn_rate(), levels = 10)

# Create the cross-validation folds
reg_vfold <- vfold_cv(reg_train, v = 10)


set.seed(123) # for reproducibility

# Create tuning workflow
reg_tune <-
  tune_grid(
    reg_wf,
    resamples = reg_vfold,
    grid = reg_grid
  )

# Select best parameter
chosen_learn_rate <- reg_tune |> select_best(metric = "rmse")

# Update workflow
reg_wf_mlp <-
  reg_wf |>
  finalize_workflow(chosen_learn_rate)

# Cross-validated workflow
reg_mlp_cv <-
  reg_wf_mlp |>
  fit_resamples(
    resamples = reg_vfold,
    control = control_resamples(save_pred = TRUE)
  )
```

Let's collect the metrics and inspect the RMSE and $R^2$ values. The code is seen in @exm-predict-reg-metrics-mlp.

::: {#exm-predict-reg-metrics-mlp}
```{r}
#| label: predict-reg-metrics-mlp

# Collect metrics
collect_metrics(reg_mlp_cv)
```
:::

So in summary, we've tried four different model specifications. The regularized linear regression model, the decision tree model, the random forest model, and the neural network model. The random forest model performed the best. For each of these models, however, we have only tried word features measured by $tf$-$idf$\index{term-frequency inverse-document frequency (tf-idf)}. We could imagine that the performance of these models could be improved by varying the features to include bigrams\index{ngrams}, for example. We could also explore different measures of word usage. Furthermore, for some of our models, we could try different engines and/ or hyperparameters (some have more than one!).

We could continue to try to explore these possible combinations, and you likely would in your research. But at this point we have a model that is performing better than the null model and is performing better than the other models we have tried. So we will consider this model to be good enough for our purposes.

Let's take our random forest model\index{random forest}, fit it to our training data, apply it to the testing data, and collect the metrics on the test set. The code is seen in @exm-predict-reg-model-spec-random-forest-workflow-tune-fit.

::: {#exm-predict-reg-model-spec-random-forest-workflow-tune-fit}
```{r}
#| label: predict-reg-model-spec-random-forest-workflow-tune-fit

# Fit the model to the training set and
# evaluate on the test set
reg_final_fit <-
  last_fit(
    reg_wf_rf,
    split = reg_split
  )

# Evaluate model on testing set
collect_metrics(reg_final_fit)
```
\index{R packages!tune}
\cindex{last_fit()}\cindex{collect_metrics()}
:::

OK. The difference between the cross-validated metrics and the metrics for the test set differ ---but only slightly. This suggests that the model is robust and that we have not overfit the data from the training set.

<!-- Step 7: Interpret -->

Now, our likely goal as an academic is to understand something about the features that contribute to the performance of the model. So let's approach extracting feature importance\index{feature importance} from the random forest model we build with the `ranger` engine. Remember, we added an `importance` argument to the `set_engine()` function and set it to 'impurity'. We can now take advantage by using {vip} to extract the feature importance. The code is seen in @exm-predict-reg-final-fit-vip.

::: {#exm-predict-reg-final-fit-vip}
```{r}
#| label: predict-reg-final-fit-vip

# Extract feature importance
reg_vip <-
  reg_final_fit |>
  extract_fit_parsnip() |>
  vi(scale = TRUE)

# Preview
reg_vip |>
  slice_head(n = 10)
```
\index{R packages!vip}\index{R packages!hardhat}\index{R packages!dplyr}
\cindex{vi()}\cindex{extract_fit_parsnip()}\cindex{slice_head()}
:::

We can now visualize the feature importance of the model. The code is seen in @exm-predict-reg-final-fit-vip-visualize.

::: {#exm-predict-reg-final-fit-vip-visualize}
```r
# Extract predictions
reg_vip |>
  mutate(Variable = str_replace(Variable, "^tfidf_text_", "")) |>
  slice_max(Importance, n = 20) |>
  # reorder variables by importance
  ggplot(aes(reorder(Variable, Importance), Importance)) +
  geom_point() +
  coord_flip() +
  labs(
    x = "Feature",
    y = "Importance"
  )
```
:::

```{r}
#| label: fig-predict-reg-final-fit-vip-visualize
#| fig-cap: "Feature importance of the random forest model"
#| fig-alt: "A dot plot which shows the importance of the features in the random forest model. The y-axis is the feature (sorted by importance) and the x-axis is the importance value."
#| fig-width: 6
#| fig-asp: 0.5
#| echo: false

# Extract predictions
reg_vip |>
  mutate(Variable = str_replace(Variable, "^tfidf_text_", "")) |>
  slice_max(Importance, n = 20) |>
  # reorder variables by importance
  ggplot(aes(reorder(Variable, Importance), Importance)) +
  geom_point() +
  coord_flip() +
  labs(
    x = "Feature",
    y = "Importance"
  ) +
  theme_qtalr(font_size = 10)
```
\index{R packages!ggplot2}\index{R packages!dplyr}\index{R packages!stringr}
\cindex{mutate()}\cindex{slice_max()}\cindex{geom_point()}\cindex{coord_flip()}\cindex{labs()}\cindex{ggplot()}\cindex{reorder()}\cindex{str_replace()}

In this section, we've built text regression models focusing on the ability to change algorithms and hyperparameters. We have also seen some of the differences between evaluating model performance between classification and regression tasks. There are many more combinations of model specifications and feature selection and engineering that can be applied. In your research, you will find yourself using these tools to explore the best model for your data.

## Activities {.unnumbered}

In the following activities, we will apply the concepts and techniques we have learned in this chapter. We will use the Tidymodels framework to build and evaluate supervised machine learning models for text classification and regression tasks.

::: {.callout}
**{{< fa regular file-code >}} Recipe**

**What**: Building predictive models\
**How**: Read Recipe 9, complete comprehension check, and prepare for Lab 9.\
**Why**: To continue to build experience building predictive models with the Tidymodels framework.
:::

::: {.callout}
**{{< fa flask >}} Lab**

**What**: Text classification\
**How**: Clone, fork, and complete the steps in Lab 9.\
**Why**: To apply your knowledge of supervised machine learning to a text classification task.
:::

## Summary {.unnumbered}

In this chapter, we outlined the workflow for approaching predictive modeling and the Tidymodels framework. We then applied the workflow to text classification and regression tasks. In the process, we gained experience identifying, selecting, and engineering features on the one hand, and building and tuning models on the other. To evaluate the models, we used cross-validation for performance and finalized our interpretation with techniques to extract feature importance.