06-techniques.qmd

# Techniques for Machine Learning Applications {#sec-06-techniques}

```{r}
#| echo: false
library(ggplot2)
book_theme <- theme_minimal() + 
  theme(plot.title=element_text(face="bold"))
ggplot2::theme_set(book_theme)
```

:::::: solutionbox
:::: solutionbox-header
::: solutionbox-icon
:::

Learning Objectives
::::

::: solutionbox-body
-   Apply feature engineering techniques to prepare data for modelling
-   Evaluate and select the most appropriate machine learning model
    based on data characteristics
-   Understand the fundamentals of key machine learning algorithms and
    their applications
:::
::::::

\
Selecting the most suitable machine learning\index{Machine learning}
model involves understanding **the goals of the analysis**, the **nature
of the data**, and the **statistical and machine learning methods** that
best suit the tasks. In @sec-05-machine_learning, we learned about what
machine learning models are, provided examples for building a model
framework, and selected common metrics for model performance calibration
and evaluation.

In this chapter, we focus on the strategies for selecting appropriate
models by leveraging the strengths of different techniques, specifically
for health metrics\index{Health metrics} and for infectious
diseases\index{Infectious diseases}. We will explore various
considerations involved in addressing potential biases, and discuss
actions to prevent them.

## Goals of the Analysis and Nature of Data

The identification of the primary **goal of the analysis** is
fundamental. Whether it involves trend analysis, investigating the
relationships between response and predictor variables, or strictly
forecasting to predict future outcomes, the strategy for model selection
varies accordingly.

**Health Metrics Data:**

-   **Composite Measures:**\index{Composite measures} Health metrics
    like DALYs\index{Disability Adjusted Life Years (DALYs)} are
    composite measures that include both mortality\index{Mortality} and
    morbidity\index{Morbidity} data, often requiring sophisticated
    regression models capable of handling continuous variables and
    multiple predictors. By examining the components of DALYs (e.g.,
    Years of Life Lost (YLLs)\index{Years of Life Lost (YLLs)} and Years
    Lived with Disability
    (YLDs)\index{Years Lived with Disability (YLDs)}), we can identify
    the key drivers such as mortality rates, disease prevalence, and
    risk factors.

-   **Regression Models:** Regression models\index{Regression models !},
    including **linear regression**\index{Regression models ! linear
    regression}, **Ridge
    regression**\index{Regression models ! Ridge regression}, and
    **Lasso regression**\index{Regression models ! Lasso regression},
    are commonly used to handle these continuous variables and address
    challenges like correlation and
    multicollinearity\index{Multicollinearity} with appropriate
    techniques such as regularisation.

**Infectious Disease Data:**

-   **Categorical and Continuous Data:** Infectious disease data can be
    categorical\index{Categorical} (e.g., disease presence or absence)
    or continuous\index{Continuous} (e.g., incidence rates).
    Classification models are suitable for categorical outcomes, while
    regression models are appropriate for continuous data.

-   **Disease Dynamics:**\index{Disease dynamics !} Understanding the
    dynamics of infectious diseases, such as transmission
    rates\index{Disease dynamics ! transmission rates}, incubation
    periods\index{Disease dynamics ! incubation periods}, and
    immunity\index{Disease dynamics ! immunity}, informs the selection
    of models. Common models include compartmental
    models\index{Compartmental model} (e.g., SIR, SEIR) and agent-based
    models.

**Common considerations for health metrics and infectious diseases data
type:**

-   **Seasonality and Trends:** The data may exhibit
    seasonality\index{Seasonality} or trends\index{Trend}, necessitating
    the use of time series\index{Time series} analysis models like
    **ARIMA**\index{AutoRegressive Integrated Moving Average (ARIMA)} or
    seasonal decomposition\index{Decomposition} to capture these
    patterns.

-   **Simulation Models:** These models can predict the impact of
    interventions on DALYs\index{Disability Adjusted Life Years (DALYs)}
    and infectious diseases\index{Infectious diseases}, estimating the
    effectiveness of different interventions and guiding policy
    decisions. Examples of these types of models are: **SIR** models,
    and **Agent-based models**\index{Agent-based models}. In addition,
    **confidence intervals**\index{Confidence intervals} and
    **sensitivity analyses**\index{Sensitivity analyses} help assess the
    uncertainty associated with these predictions.

-   **Bayesian Models:**\index{Bayesian models} These models can
    estimate parameters and make predictions based on prior knowledge
    and observed data, incorporating uncertainty and variability.

-   **Predictive modelling:**\index{Predictive modelling} Such as
    **decision trees**\index{Decision trees}, **support vector machines
    (SVM)**\index{Support Vector Machines (SVM)}, and **Long Short-Term
    Memory (LSTM) neural
    networks**\index{Long Short-Term Memory (LSTM)}, can predict disease
    outbreaks, identify high-risk populations, and optimise resource
    allocation.

## Statistical and Machine Learning Methods

The choice of model depends on the type of data, the relationships
between variables, and the goals of the analysis. Once we have these
factors well identified, we are a step forward in restricting the range
of applicable models.

The next step involves conducting a thorough **exploratory data analysis
(EDA)**\index{Exploratory Data Analysis (EDA)}. This initial exploration
helps to uncover the underlying structure of the data, the relationships
between variables, and the way the response variable—which may also be
referred to as the outcome variable—depends on predictors. This phase is
critical as it informs the necessity of subsequent data adjustments and
transformations.

The importance of data preparation and exploratory data analysis in
machine learning are the building blocks in the preparation of machine
learning digestible data. **Feature
engineering**\index{Feature engineering} is a technique that involves
creating new features from existing ones based on domain knowledge or
transformation of data to improve the model's ability to discern
patterns. For example, creating features like **moving
averages**\index{Moving averages} or differences between consecutive
days can reveal trends and cycles\index{Cycles} that are not immediately
apparent from raw data.

Another example is the **standardisation**\index{Standardisation}
process, which is crucial when dealing with variables measured in
different units. It involves rescaling the features so they have a mean
of zero and a standard deviation of one. This process is particularly
important when variables span several orders of magnitude; without
standardisation, a model might incorrectly interpret the scale of a
feature as a proxy for importance.

Furthermore, the application of transformations, such as **logarithmic**
scaling or the application of **spline** functions can help in managing
skewed data or enhancing model ability to capture non-linear
relationships, which result particularly useful in complex data
modelling. In addition, tailored adjustments, and more sophisticated
manipulations have been implemented over time to allow estimation of
missing values in order to obtain customised, flexible, and more
homogeneous data. For more information on feature engineering, see
[@butcher2020] useful for effective machine learning strategy
application, covering various techniques and appropriate use cases,
focusing on practical understanding and implementation.

## Model Selection Strategies

In developing predictive models for health metrics and infectious
diseases, selecting the appropriate model is critical to ensure accurate
and reliable forecasts. Here are outlined sample strategies employed in
the model selection process, we introduce the **Rabies**\index{Rabies}
dataset used for our discussion and demonstrate the selection of a
suitable model for analysing its impact. Rabies, although nearly 100%
fatal once symptoms appear, presents a unique challenge due to the
relative rarity of cases and limited availability of comprehensive data.
This scarcity complicates efforts to model the disease accurately and
develop effective public health strategies.

To address these challenges, we explore advanced modelling techniques
that can enhance the robustness of our analyses despite data
limitations, which involves evaluating multiple models based on their
performance and selecting the best-fitting models to achieve the most
accurate predictions.

## Example: Rabies {#sec-06-rabies}

The **rabies**\index{Rabies} dataset from the `{hmsidwR}` package
contains information on death rates and disability-adjusted life years
(DALYs) per 100,000 inhabitants due to rabies and all causes of
mortality in Asia and for the Global region from 1990 to 2019. Rabies
([@cdc2024]) is a fatal viral infection, and it is also classified as an
infectious disease that can infect all mammals causing acute
encephalitis. Caused by the rabies virus, which belongs to the
Lyssavirus genus, it is transmitted to humans through the bite of an
infected animal such as bats, raccoons, skunks, foxes, and obviously
dogs, which are the main source of human rabies deaths [@hampson2015].
Rabies defined as **neglected tropical disease
(NTD)**\index{Neglected Tropical Disease (NTD)} predominantly affects
already marginalised, poor and vulnerable populations. Although
effective human vaccines and immunoglobulins exist for rabies, these are
often not readily available or accessible to everyone [@rabies].

In this example we consider the number of DALYs per 100,000 inhabitants
due to rabies in Asia and the Global region, as our response variable,
the dataset is made available in the `{hmsidwR}` package. It is composed
of 240 observations and 7 variables: `measure`, `location`, `cause`,
`year`, `val`, `upper`, `lower`.

```{r}
library(tidyverse)
hmsidwR::rabies %>%
  filter(year >= 1990 & year <= 2019) %>%
  select(-upper, -lower) %>%
  head()
```

Selecting only the `cause == Rabies` , the first thing to notice is that
deaths rates and DALYs are on different units, rates and years
respectively.

```{r}
library(tidyverse)
rabies <- hmsidwR::rabies %>%
  filter(year >= 1990 & year <= 2019) %>%
  select(-upper, -lower) %>%
  pivot_wider(names_from = measure, values_from = val) %>%
  filter(cause == "Rabies") %>%
  rename(dx_rabies = Deaths, dalys_rabies = DALYs) %>%
  select(-cause)

rabies %>% head()
```

```{r}
#| echo: false
#| eval: false
breaks <- seq(min(rabies$year), max(rabies$year), by = 10)

rabies %>%
  mutate(year_10 = cut(year,
                       breaks = breaks,
                       right = FALSE,
                       labels = paste(breaks[-length(breaks)],
                                      breaks[-1] - 1,
                                      sep = "-")),
         year_10 = case_when(is.na(year_10) ~ "2019",
                             TRUE ~ as.factor(year_10))) %>%
  group_by(location, year_10) %>%
  summarize(avg_10yrs_dalys = mean(dalys_rabies))
```

It can be seen that the number of deaths due to rabies is much lower
than the number of DALYs. This difference in scale can affect the
model's ability to learn from the data. To address this issue, we can
scale and centre the numeric variables to make them more comparable.

```{r}
p1 <- rabies %>%
  ggplot(aes(x = year, group = location, linetype = location)) +
  geom_line(aes(y = dx_rabies),
            linewidth = 1) +
  geom_line(aes(y = dalys_rabies))

p2 <- rabies %>%
  # apply a scale transformation to the numeric variables
  mutate(year = as.integer(year),
         across(where(is.double), scale)) %>%
  ggplot(aes(x = year, group = location, linetype = location)) +
  geom_line(aes(y = dx_rabies),
            linewidth = 1) +
  geom_line(aes(y = dalys_rabies))
```

```{r}
#| layout-ncol: 2
#| label: fig-rabies-deaths-dalys
#| fig-cap: "Not Scaled and Scaled and Centred"
#| fig-subcap:
#|   - "Not scaled"
#|   - "Scaled and centred"
#| fig-alt: "DALYs due to Rabies - Not scaled and Scaled and centred"
#| echo: false

legend_grob <- cowplot::get_legend(
  ggplot(data = rabies %>%
           pivot_longer(cols = c("dx_rabies", "dalys_rabies")) %>%
           mutate(name = ifelse(name == "dx_rabies", "Deaths", "DALYs")),
         aes(x = year, group = name, linewidth = name)) +
    geom_line(aes(y = value),
              color = "grey") +
    labs(linewidth = ""))

p1 +
  labs(
    title = "Dalys and Deaths due to Rabies - Not scaled",
    subtitle = "in Asia and Global Region",
    color = "Location",
    y = "Value", x = "Time(Year)") +
  scale_color_manual(values = c("orange", "navy"))

p2 +
  labs(
    title = "Dalys and Deaths due to Rabies - Scaled and centred",
    subtitle = "in Asia and Global Region",
    color = "Location",
    y = "Value", x = "Time(Year)") +
  scale_color_manual(values = c("orange", "navy")) +
  annotation_custom(legend_grob,
    xmin = 2010, xmax = 2015,
    ymin = 1, ymax = 2)
```

Creating new features from existing ones provide additional predictive
power. Then, combine the cause vector in a way to obtain two vectors for
death rates due to rabies and all causes, scale and centre the numeric
variables to obtain homogeneous data to use in the model.

```{r}
all_causes <- hmsidwR::rabies %>%
  filter(year >= 1990 & year <= 2019) %>%
  select(-upper, -lower) %>%
  pivot_wider(names_from = measure, values_from = val) %>%
  filter(!cause == "Rabies") %>%
  rename(dx_allcauses = Deaths, dalys_allcauses = DALYs) %>%
  select(-cause)

dat <- rabies %>%
  full_join(all_causes, by = c("location", "year"))

dat %>% head()
```

To be able to visualise the magnitude of difference between death rates
and DALYs for both rabies and all causes, it is necessary to scale or
standardise the data as shown above.

```{r}
p3 <- dat %>%
  select(-year, -location) %>%
  scale() %>%
  cbind(dat %>% select(year, location)) %>%
  ggplot(aes(x = year, 
             group = location, 
             linetype = location)) +
  geom_line(aes(y = dx_rabies),
            linewidth = 1) +
  geom_line(aes(y = dx_allcauses))

p4 <- dat %>%
  select(-year, -location) %>%
  scale() %>%
  cbind(dat %>% select(year, location)) %>%
  ggplot(aes(x = year,
             group = location,
             linetype = location)) +
  geom_line(aes(y = dalys_rabies),
            linewidth = 1) +
  geom_line(aes(y = dalys_allcauses))
```

```{r}
#| layout-ncol: 2
#| label: fig-deaths-dalys-std
#| fig-cap: "Scaled and centred"
#| fig-subcap:
#|  - "Deaths due to Rabies and All Causes"
#|  - "Dalys due to Rabies and All Causes"
#| fig-alt: "Deaths and DALYs due to Rabies and All Causes - Scaled and centred"
#| echo: false

legend_grob2 <- cowplot::get_legend(
  ggplot(data = dat %>%
           select(-year, -location) %>%
           scale() %>%
           cbind(dat %>% select(year, location)) %>%
           pivot_longer(cols = c("dx_rabies", "dx_allcauses")) %>%
           mutate(name = ifelse(name == "dx_rabies", "Rabies", "All Causes")),
         aes(x = year, group = name, linewidth = name)) +
    geom_line(aes(y = value), color = "grey") +
    labs(linewidth = ""))

p3 +
  labs(
    title = "Deaths due to Rabies and All Causes",
    subtitle = "in Asia and Global Region",
    color = "Location",
    y = "Value", x = "Time(Year)") +
  scale_color_manual(values = c("orange", "navy")) +
  annotation_custom(legend_grob2,
    xmin = 2010, xmax = 2015,
    ymin = 1, ymax = 2)

p4 +
  labs(
    title = "Dalys due to Rabies and All Causes",
    subtitle = "in Asia and Global Region",
    color = "Location",
    y = "Value", x = "Time(Year)") +
  scale_color_manual(values = c("orange", "navy")) +
  annotation_custom(legend_grob2,
    xmin = 2010, xmax = 2015,
    ymin = 1, ymax = 2)
```

For this task, we will use the `{tidymodels}` meta-package, as it
provides a consistent interface for modelling and machine learning
tasks. In particular, we define and execute modelling workflows, to
create tailored data pre-processing tasks on various modelling
specifications, and evaluate the performance using resampling
techniques, to eventually select the best model. A more detailed
explanation of the `{tidymodels}` framework can be found in the book
[@silge].

### Training Data and Resampling

Splitting data into training and test allows the model to train a
subsection of the data and then test its performance on the remaining
part of the data, the test set. In this case, we will use the
`initial_split()` function to split the data into training and test
sets. The proportion assigned to trains can vary but it is usually
assigned to be 80%, also a stratification option can be set.

```{r}
library(tidymodels)

set.seed(11012024)
split <- initial_split(dat, prop = 0.8, strata = location)
training <- training(split)
test <- testing(split)
```

After that, it is important to create a set of folds, which means a set
of subgroups of the original data by grouping following specific
directions based on the type of **resampling
technique**\index{Resampling !}. Resampling techniques are used to
evaluate the model's performance and estimate its generalisation error.
There are various types of resampling techniques, it depends on the
specific characteristics of your dataset, and the goals of your
analysis. Some of the most common resampling techniques include:

-   **k-Fold Cross-Validation**\index{k-Fold Cross-Validation} for
    general model evaluation and hyperparameter tuning.
-   **Bootstrap Resampling**\index{Resampling ! bootstrap} to estimate
    the variability of your model and for smaller datasets.
-   **Time Series Cross-Validation**\index{Time series cross validation}
    for time-dependent data to preserve temporal structure.
-   **Spatial Resampling**\index{Resampling ! spatial} for spatially
    correlated data to account for spatial dependencies.
-   **Stratified Resampling**\index{Resampling ! stratified} when
    dealing with imbalanced datasets to ensure proper representation of
    all classes.

In this case, we will use **k-Fold Cross-Validation** to evaluate the
model's performance. The `vfold_cv()` function creates a set of folds
for cross-validation\index{Cross Validation}, which is used to train and
test the model on different subsets of the data.

```{r}
set.seed(11102024)
folds <- vfold_cv(training, v = 10)
```

### Data Preprocessing and Featuring Engineering

As already seen in the exploratory phase, preprocessing data is a
crucial step in machine learning\\index{Machine learning. This process
can include techniques for handling missing values, standardisation of
the data, encoding categorical variables, and removing highly correlated
variables.

In this case, we will use the `{recipes}` package to create a recipe,
with a set of preprocessing steps. The `recipe()` function allows us to
define a model formula and use various `step_<functions>()` for
manipulating data. We are going to set up 3 recipes, the first is a
basic one which includes all variables and does not perform any data
transformation.

```{r}
rec <- recipe(dalys_rabies ~ ., data = training)
```

The second recipe includes some key steps to transform the data into a
way specific models would be able to understand and learn from it.
Models such as k-nearest neighbours, or support vector machines, that
rely on distance metrics, can be sensitive to differences in feature
scales.

For instance, non-standardised year data can dominate the model's
decision-making process, leading to biased results. By scaling and
centring the data, we ensure that all features contribute equally to the
model's predictions.

We can create more complex recipes with more steps, but for this
example, we will use a step for encoding the location variable (Asia,
Global) into a numeric vector, and a second step to normalise (or
standardise) all predictors.

```{r}
rec1 <- recipe(dalys_rabies ~ ., data = training) %>%
  # convert nominal variables to dummy variables
  step_dummy(all_nominal_predictors()) %>%
  # scale the numeric variables
  step_normalize(all_numeric_predictors())
```

Once the recipe is created, we can apply it to the data using the
`prep()` function, which estimates the necessary parameters for the
transformations and applies them to the data. Then, to check the results
we can use the `juice()` function to extract the processed data.

```{r}
rec1 %>%
  prep() %>%
  juice() %>%
  select(1, 2, 5) %>%
  head()
```

Trained data can be also tested on new data, in this case we test them
on the `test` set with the `bake()` function.

```{r}
rec1 %>%
  prep() %>%
  bake(new_data = test) %>%
  select(1, 2, 5) %>%
  head()
```

DALYs often aggregate various health impacts, and can have highly skewed
distributions. This skewness arises due to several factors: the presence
of outliers, the nature of the health condition being measured, and the
distribution of the data itself. To handle the skewness of the data, we
can apply:

-   Log Transformation: $log(DALYs+1)$
-   Sqrt Transformation: $\sqrt{DALYs}$
-   Yeo-Johnson Transformation, a generalisation of the Box-Cox
    transformation that can handle both positive and negative values:
    $((DALYs+1)^p-1)/p$.

Let's apply the **Yeo-Johnson
transformations**\index{Yeo-Johnson transformations} to the response
variable (`dalys_rabies)` and see how the density distribution changes
with different values of $\lambda$. This is a step that can be tuned
with a machine learning algorithm.

```{r}
#| layout-ncol: 3
#| label: fig-response-transformation
#| fig-cap: "Response variable transformation"
#| fig-subcap: 
#|  - "Log10"
#|  - "Yeo-Johnson p=-2"
#| fig-alt: "Log10, Yeo-Johnson p=-2, Yeo-Johnson p=2"
#| echo: false

ggplot(dat) +
  geom_density(aes(x = dalys_rabies)) +
  scale_x_log10() +
  labs(title = "Log10 transformation")

ggplot(dat) +
  geom_density(aes(x = dalys_rabies)) +
  scale_x_continuous(transform = scales::transform_yj(p = -2)) +
  labs(title = "Yeo-Johnson transformation p=-2")

ggplot(dat) +
  geom_density(aes(x = dalys_rabies)) +
  scale_x_continuous(transform = scales::transform_yj(p = 2)) +
  labs(title = "Yeo-Johnson transformation p=2")
```

Let's now create a third recipe with the `step_YeoJohnson()` function.

```{r}
rec2 <- rec1 %>%
  # apply Yeo-Johnson transformation to the response variable
  step_YeoJohnson(dalys_rabies)

rec2 %>%
  prep() %>%
  juice() %>%
  select(1, 2, 5) %>%
  head()
```

### Correlation, Multicollinearity and Overfitting

To be noted is that we haven't applied any
**correlation**\index{Correlation} selection step on this data.
Filtering out highly correlated predictors, such as those with a
correlation greater than 80% to avoid
multicollinearity\index{Multicollinearity}, would lead to excluding
crucial variables. On the other hand, including all possible covariates
in a model often yields implausible signs on covariates or unstable
coefficients, as well as overfitting [@foreman2012].

When multiple predictors are correlated, but all are crucial for the
analysis (e.g., deaths due to rabies, total deaths, and total DALYs for
all causes), applying a correlation step that filters out correlated
variables can be problematic. One way to overcome bias arising from it
is using regularisation techniques like **Ridge
Regression**\index{Regression ! ridge} or **Lasso
Regression**\index{Regression ! lasso} is often the best approach to
handle multicollinearity without removing any predictors. Alternatively,
**Principal Components Analysis
(PCA)**\index{Principal Components Analysis (PCA)} can reduce
dimensionality while retaining most of the variance. These methods
ensure all important predictors are included in the model without the
adverse effects of multicollinearity.

### Model Specification

The next step is to outline the model specification. There are various
type of models that can be used. We start with a **random
forest**\index{Random Forest}. This choice is typically done due to the
algorithm's features, which is able to create **multiple bootstrap
samples**\index{Bootstrap} (random samples with replacement) from the
original dataset. Each bootstrap sample is used to train a separate
decision tree.

### Model 1: Random Forest

Rabies death rates may exhibit complex relationships with predictor
variables. Random forests\index{Random Forests} are capable of capturing
non-linear relationships between predictors and the target variable.

Also, it handles multicollinearity, missing data, provides variables
importance and is an ensemble learning method, which means they combine
the predictions of multiple individual decision trees to produce a more
accurate and stable prediction.

In our simplified case this type of model will do random samples with
replacement of data. In `{tidymodels}` we can select different types of
engines, in the case of random forest we could use random forest,
ranger, and others. The difference between these engines derives from
the specific type of calculation used to make the estimation. The Ranger
engine is notably faster than random forest, so let's use that for this
example.

```{r}
rf_mod <- rand_forest(mtry = tune(),
                      trees = tune(),
                      min_n = tune(),
                      mode = "regression",
                      engine = "ranger")

wkf <- workflow(preprocessor = rec, 
                spec = rf_mod)

rf_res <- tune_grid(object = wkf,
                    resamples = folds,
                    grid = 5,
                    control = control_grid(save_pred = TRUE))

show_best(rf_res, metric = "rmse") %>%
  select(-n, -std_err)
```

```{r}
rf_res_tuned <- select_best(rf_res, metric = "rmse")

rf_res_tuned
```

```{r}
rf_fit <- wkf %>%
  finalize_workflow(select_best(rf_res,
                                metric = "rmse")) %>%
  fit(training)

rf_fit %>%
  predict(new_data = test) %>%
  bind_cols(test) %>%
  rmse(truth = dalys_rabies, estimate = .pred)
```

```{r}
#| layout-ncol: 2
#| label: fig-rf-predictions
#| fig-cap: "Predictions vs. Truth"
#| fig-subcap:
#|   - "Predictions vs. Truth"
#|   - "Predictions vs. Truth by Year"
#| fig-alt: "Predictions vs. Truth and Predictions vs. Truth by Year"   
#| echo: false

rf_fit %>%
  predict(new_data = test) %>%
  bind_cols(test) %>%
  ggplot(aes(dalys_rabies, .pred,
             color = location)) +
  geom_abline(color = "brown") +
  geom_point(size = 2) +
  scale_color_manual(values = c("orange", "navy")) +
  coord_fixed() +
  expand_limits(x = 0, y = 0) +
  labs(
    title = "Predictions vs. Truth",
    x = "Truth",
    y = "Predictions")

rf_fit %>%
  predict(new_data = test) %>%
  bind_cols(test) %>%
  ggplot(aes(x = year, 
             color = location)) +
  geom_point(aes(y = dalys_rabies)) +
  geom_line(aes(y = .pred)) +
  scale_color_manual(values = c("orange", "navy")) +
  labs(
    title = "Predictions vs. Truth by Year",
    x = "Time(Year)",
    y = "Values")
```

### Model 2: Generalised Linear Model (GLM)

Generalised Linear Models (GLMs)\index{Generalised Linear Models (GLMs)}
involve statistical estimation rather than the iterative parameter
tuning, common in many machine learning techniques. However, adding a
machine learning feature through parameter calibration can be done using
techniques such as cross-validation and grid search to find the best
model settings.

To introduce a machine learning feature with parameter calibration into
our modelling of the rabies data, we can use a technique like
cross-validation combined with a **regularisation method** or an
algorithm that supports parameter tuning. Here, we can employ a model
from the `glmnet` package, which fits a generalised linear model via
**penalised maximum likelihood**. The regularisation path is computed
for the lasso or elastic-net penalty at a grid of values for the
regularisation parameter lambda.

Adding Machine Learning Features with `{glmnet}` and Cross-Validation

```{r}
if (!require(glmnet)) install.packages("glmnet")
library(glmnet)
```

For `glmnet`, we need to input matrices rather than data frames, and
create matrices for the independent variables (predictors) and the
dependent variable (response).

```{r}
predictors <- model.matrix(dalys_rabies ~ .,
                           data = dat)[, -1] # Remove intercept
response <- dat$dalys_rabies
```

Use cross-validation to find the optimal lambda value, which controls
the strength of the regularisation:

```{r}
# Set seed for reproducibility
set.seed(123)

# Fit the model with cross-validation
cv_model <- cv.glmnet(predictors,
                      response,
                      family = "gaussian")
cv_model
```

Extracting the best model, we can see that $\lambda$ is 0.165.

```{r}
# Get the best lambda value
best_lambda <- cv_model$lambda.min
paste("Best Lambda:", best_lambda)
```

And plot the lambda selection with the `plot()` function.

```{r}
#| label: fig-glmnet-predictions
#| fig-cap: "Cross-Validation Optimal Lambda"
#| fig-alt: "Lambda Selection"
# Plot the lambda selection
plot(cv_model)
```

Then, fitting the final model with the selected best lambda, we can
predict and evaluate the model.

```{r}
final_model <- glmnet(predictors,
                      response,
                      family = "gaussian",
                      lambda = best_lambda)

# Predict using the final model
predictions <- predict(final_model,
                       # values of the penalty parameter lambda
                       s = best_lambda,
                       # matrix of new values for x
                       newx = predictors
                       )

# Calculate Mean Squared Error
rmse <- sqrt(mean((response - predictions)^2))
paste("Root Mean Squared Error:", rmse)
```

By incorporating `glmnet` and using lambda selection via
cross-validation, we introduce a **machine learning feature—parameter
calibration** into our analysis. This approach not only helps in
minimising overfitting but also enhances model performance by selecting
the most effective regularisation parameter. The cross-validation
process used here is crucial for confirming that our model's parameters
are optimally tuned for the given data, embodying a key aspect of
machine learning methodologies.

```{r}
#| echo: false
#| layout-ncol: 2
#| label: fig-glmnet-predictions2
#| fig-cap: "Predictions vs. Truth"
#| fig-subcap: 
#|  - "Predictions vs. Truth"
#|  - "Predictions vs. Truth by Year"
#| fig-alt: "Predictions vs. Truth and Predictions vs. Truth by Year"

# plot the predictions
data.frame(dalys_rabies = response, predictors, predictions) %>%
  ggplot(aes(x = dalys_rabies, y = predictions)) +
  geom_point() +
  geom_abline(color = "brown") +
  labs(
    title = "Predictions vs Truth",
    x = "Truth",
    y = "Predictions")

data.frame(dalys_rabies = response, predictors, predictions) %>%
  ggplot(aes(x = year, 
             group = factor(locationGlobal),
             color = factor(locationGlobal))) +
  geom_point(aes(y = dalys_rabies)) +
  geom_line(aes(y = predictions), color = "brown") +
  scale_color_manual(values = c("orange", "navy"),
                     labels = c("Asia", "Global")) +
  labs(
    title = "DALYs due to Rabies Predictions vs Truth",
    subtitle = "by Year",
    x = "Truth",
    y = "Predictions")
```

### Testing Multiple Models

In the example above, we used two models to predict DALYs due to rabies,
a random forest with `{tidymodels}` and a generalised linear model with
`{glmnet}` with a Root Mean-Square
Error\index{Root Mean-Square Error (RMSE)} of 0.448 and 0.257
respectively. The Random Forest model has a higher RMSE, which means it
has a higher prediction error compared to the GLM model. However, we
haven't applied any of the preprocessing steps, and there are many other
models that could be used to predict DALYs such as:

1.  **Support Vector Machines
    (SVM)**\index{Support Vector Machines (SVM)}: SVMs are a powerful
    machine learning algorithm that can be used for both classification
    and regression tasks. They work by finding the hyperplane that best
    separates the data into different classes or groups.

2.  **Extreme Gradient Boosting
    (XGBoost)**\index{Extreme Gradient Boosting (XGBoost)}: Known for
    its high performance in various prediction tasks, XGBoost can handle
    missing values and is effective for large datasets.

3.  **K-Nearest Neighbours (KNN)**\index{K-Nearest Neighbours (KNN)}
    models are a type of instance-based learning algorithm that stores
    all available cases and classifies new cases based on a similarity
    measure.

4.  **Long Short-Term Memory (LSTM)
    Networks**\index{Long Short-Term Memory (LSTM)}: For temporal or
    sequential health data, LSTM networks can capture dependencies over
    time, making them suitable for time-series prediction of disease
    progression and outcomes.

Each of these models has its own strengths and weaknesses, and the
choice of model will depend on the specific characteristics of the data
and the goals of the analysis. By testing multiple models and comparing
their performance, we can identify the best model for the given data and
task.

Let's use the `{parsnip}` package and the `workflow_set()` function to
fit a set of models to the rabies data. We will fit a **Support Vector
Machine (SVM)**, and a **K-Nearest neighbours (KNN) model** to the data
and compare their performance.

```{r}
linear_reg_spec <- 
  linear_reg(penalty = tune(),
             mixture = tune()) %>%
  set_engine("glmnet")

svm_p_spec <-
  svm_poly(cost = tune(),
           degree = tune()) %>%
  set_engine("kernlab") %>%
  set_mode("regression")

knn_spec <-
  nearest_neighbor(neighbor = tune(),
                   dist_power = tune(),
                   weight_func = tune()) %>%
  set_engine("kknn") %>%
  set_mode("regression")
```

```{r}
#| eval: false
library(rules)
library(baguette)
# Combine workflows into a workflow set
workflow_set <- workflow_set(preproc = list(scaled = rec1,
                                            yeo_johnson = rec2),
                             models = list(linear_reg = linear_reg_spec,
                                           svm = svm_p_spec,
                                           knn = knn_spec))

grid_ctrl <-control_grid(save_pred = TRUE,
                         parallel_over = "everything",
                         save_workflow = TRUE)

# Fit and evaluate the models with hyperparameter tuning
grid_results <- workflow_set %>%
  workflow_map(seed = 1503,
               resamples = folds,
               grid = 5,
               control = grid_ctrl)
```

```{r}
#| eval: true
#| echo: false
# saveRDS(grid_results, file = "data/model-data/grid_results.rds")
grid_results <- readRDS("data/model-data/grid_results.rds")
```

```{r}
# Show the results
grid_results %>%
  collect_metrics() %>%
  arrange(mean) %>%
  select(1, 5, 7, 9) %>%
  head()
```

```{r}
#| fig-cap: "Model Performance"
#| fig-alt: "Model Performance"
#| label: fig-model-performance
autoplot(grid_results,
         rank_metric = "rmse",
         metric = "rmse",
         select_best = TRUE) +
  geom_text(aes(y = mean - 0.1, 
                label = wflow_id),
            angle = 90, 
            hjust = 1, 
            color = "black",
            size = 3.5) +
  lims(y = c(-1.5, 0.9)) +
  theme(legend.position = "none")
```

## Summary

The integration of machine learning techniques into public
health\index{Public health} data analysis can significantly enhance the
predictive power and robustness of models. By leveraging the
capabilities of machine learning algorithms, we can extract valuable
insights from complex health data, enabling more informed
decision-making and policy formulation in public health contexts. The
examples provided in this chapter illustrate the application of machine
learning techniques to health metrics data, demonstrating the importance
of feature engineering\index{Feature engineering}, model selection, and
parameter calibration\index{Calibration} in enhancing the predictive
accuracy and relevance of models. By following best practices in machine
learning, public health researchers and practitioners can harness the
power of data-driven insights to address critical health challenges and
improve population health outcomes.

**Best Practices for Machine Learning in Public Health:**

-   Conduct exploratory data analysis to understand the underlying
    structure of the data and relationships between variables.
-   Apply feature engineering techniques to create new variables and
    enhance the model's predictive power.
-   Select machine learning models that are contextually appropriate and
    robust for public health data analysis. Such as Random Forest,
    Generalised Linear Models, and others.
-   Use parameter calibration techniques such as cross-validation,
    regularisation, monte carlo, and grid search to optimise model
    performance.
-   Evaluate model performance using appropriate metrics and
    visualisation tools to assess predictive accuracy and relevance.

The integration of machine learning methodologies into public health
data analysis represents a significant opportunity to advance the field
of public health and enhance our understanding of health metrics and
disease dynamics.