15-malaria.qmd

# The Case of Malaria {#sec-15-malaria}

```{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
-   Understand the key characteristics and impact of Malaria
-   Explore how Malaria spreads through populations and environments
-   Learn to map and visualize Malaria outbreaks using spatial data
:::
::::::

\
In this chapter, we explore the dynamics of Malaria\index{Malaria}
transmission in more detail. We examine the results of various model’s
applications that simulate the spread of the virus.

## Epidemiology

**Malaria** is a mosquito-borne\index{Mosquito borne} infectious disease
that affects humans and other animals. It is caused by parasitic
protozoans belonging to the Plasmodium type. The disease is
**transmitted through the bites of Anopheles mosquitoes**. The symptoms
of malaria typically include fever, fatigue, vomiting, and headaches. If
left untreated, malaria can be fatal. Malaria is a major public
health\index{Public health} concern in many tropical and subtropical
regions, particularly in Africa.

Malaria transmission dynamics\index{Transmission dynamics} are
influenced by various factors, including the prevalence of infected
individuals, the density of mosquito vectors, and environmental
conditions. The transmission of malaria occurs when an infected mosquito
bites a human host, injecting the Plasmodium parasites into the
bloodstream. The parasites then multiply within the host's red blood
cells, leading to the characteristic symptoms of malaria. The parasites
can be transmitted back to mosquitoes when they feed on infected
individuals, completing the transmission cycle.

## Mapping Malaria Outbreaks

Mapping malaria outbreaks\index{Outbreak} is essential for locating
high-risk areas and guiding effective public health interventions.
Geographic Information Systems (GIS)\index{Geographic Information
Systems (GIS)} enable the visualisation of malaria’s spatial
distribution, revealing transmission
patterns\index{Transmission patterns} and hotspots. The analysis of
malaria cases, alongside environmental factors, such as temperature,
humidity, and vegetation cover, allows the identification of conditions
that facilitate transmission and support the development of targeted
control strategies.

In this example, we will use the `malariaAtlas` package to download
malaria data for Nigeria and visualise the distribution of malaria cases
in the country. We will plot the malaria hotspots on a map of Nigeria to
identify regions with the highest incidence of the disease.

```{r}
#| message: false
#| warning: false
# Load necessary libraries
library(malariaAtlas)
library(tidyverse)
library(sf)
library(rnaturalearth)
```

```{r}
#| output: false
tidyverse_conflicts()
```

To download malaria data we can use the `getPR()` function, it releases
the data from the Malaria Atlas Project API. The function requires the
country and species of Plasmodium to be specified. In this example, we
will download data for Nigeria and the *Plasmodium falciparum species*.

```{r}
#| eval: false
nigeria_data <- getPR(country = "Nigeria", 
                      species = "Pf")
```

```{r}
#| eval: false
#| echo: false
# saveRDS(nigeria_data, "data/inst/extdata/15-nigeria_data.rds")
```

```{r}
#| echo: false
nigeria_data <- readRDS("data/inst/extdata/15-nigeria_data.rds")
```

Data contains cases for 23 years spanning from 1985 to 2018. We can
extract the relevant information (year_start (of the survey), longitude,
latitude, and number of positive cases) and convert the data to a
spatial object using the `st_as_sf()` function from the `sf` package.

```{r}
nigeria_data <- nigeria_data %>%
  arrange(year_start) %>%
  select(year_start,longitude,latitude,positive) %>%
  filter(!is.na(longitude) & !is.na(year_start)) 

nigeria_data_sf <- nigeria_data %>%
  st_as_sf(coords = c("longitude", "latitude"), 
                       crs = 4326)

head(nigeria_data_sf)
```

To get the administrative boundaries of Nigeria we use the
`{rnaturlaearth}` package:

```{r}
nigeria_sf <- ne_countries(country = "Nigeria", 
                           returnclass = "sf")
```

And finally plot the malaria hotspots on a map of Nigeria to visualise
the distribution of malaria cases in the country.

```{r}
#| label: fig-malaria-map
#| fig-cap: "The map shows the distribution of malaria cases in Nigeria, with the size and color of the points representing the number of positive cases. The map highlights regions with a high incidence of malaria, providing valuable information for public health interventions."
#| fig-alt: "Map showing the distribution of malaria cases in Nigeria."
# Plot Malaria Hotspots with expanded spatial data range
ggplot() +
  geom_sf(data = nigeria_sf, 
          fill = "gray90", 
          color = "white")+
geom_sf(data = nigeria_data_sf, 
        aes(size=positive, color = positive),
        alpha=0.5)+
  scale_color_viridis_c(option = "plasma", name = "Malaria Cases") +
  guides(size = "none") +
  labs(title = "Malaria Hotspots in Nigeria",
       subtitle = "1985-2018",
       caption = "Source: Malaria Atlas Project | @fgazzelloni") +
  theme(legend.position = "right")
```

Spatial data on climate, population density, and mosquito habitats are
instrumental in predicting high-risk areas for malaria, aiding
preventive measures like bed net distribution and indoor residual
spraying. In Tanzania, for example, spatial models identified villages
with the highest malaria incidence, allowing the government and NGOs to
prioritize these areas for intervention[@selemani2016]. This targeted
approach enhances the cost-effectiveness of malaria control programs,
reducing the disease burden by directing resources to areas with the
greatest need.

```{r}
#| echo: false
#| label: fig-malaria-tanzania
#| fig-cap: "The map shows the distribution of malaria cases in Tanzania, with the size and color of the points representing the number of positive cases. The map highlights regions with a high incidence of malaria, providing valuable information for public health interventions."
#| fig-alt: "Map showing the distribution of malaria cases in Tanzania."  
#tanzania_data <- getPR(country = "Tanzania", species = "Pf")
#saveRDS(tanzania_data, "data/inst/extdata/15-tanzania_data.rds")  
tanzania_data <- readRDS("data/inst/extdata/15-tanzania_data.rds")
tanzania_data_sf <- tanzania_data %>%
  arrange(year_start) %>%
  select(year_start,longitude,latitude,positive) %>%
  filter(!is.na(longitude) & !is.na(year_start)) %>%
  st_as_sf(coords = c("longitude", "latitude"), 
                       crs = 4326)

tanzania_sf <- ne_countries(country = "United Republic of Tanzania", 
                                        returnclass = "sf")
ggplot() +
  geom_sf(data = tanzania_sf, 
          fill = "gray90", 
          color = "white")+
geom_sf(data = tanzania_data_sf, 
        aes(size=positive, color = positive),
        alpha=0.5)+
  scale_color_viridis_c(option = "plasma", name = "Malaria Cases") +
  guides(size = "none") +
  labs(title = "Malaria Hotspots in Tanzania",
       subtitle = "1985-1994",
       caption = "Source: Malaria Atlas Project | @fgazzelloni") +
  theme(legend.position = "right")
```

## Example: Simulating Malaria Transmission Dynamics

In this example, we will use data for Malaria positive cases in Nigeria
from the **Malaria Atlas Project**\index{Malaria Atlas Project} to
evaluate the transmission dynamics\index{Transmission dynamics} using a
simple mathematical model. We will use machine learning to predict
future trends.

To illustrate the modelling and prediction process we will use the
`{caret}` package, which provides a unified interface for training and
evaluating machine learning models. We will train a Random Forest model
to predict future trends in malaria transmission based on historical
data.

```{r}
nigeria_cases <- nigeria_data %>% 
  group_by(year_start) %>%
  reframe(positive = sum(positive)) %>%
  rename(year = year_start) %>%
  drop_na() %>%
  select(year, positive)

head(nigeria_cases)
```

Visualize the dynamics of malaria transmission in Nigeria using a line
plot to show the number of infected cases over time.

```{r}
#| label: fig-malaria-dynamics
#| fig-cap: "The plot shows the dynamics of malaria transmission in Nigeria from 1985 to 2018. The blue line represents the number of infected cases over time, highlighting the fluctuations in malaria incidence during this period."
#| fig-alt: "Line plot showing the dynamics of malaria transmission in Nigeria from 1985 to 2018."
#| fig-width: 6
#| fig-height: 4
nigeria_cases %>%
  ggplot() +
  geom_line(aes(x = year, y = positive), 
            color = "navy", 
            linetype = "solid") +
  labs(title = "Malaria Transmission Dynamics in Nigeria",
       subtitle = "23 years spanning from 1985 to 2018",
       x = "Time(Year)",
       y = "Number of Infected Cases") +
  theme_minimal()
```

### Modelling with caret

To train a machine learning model to predict future trends in malaria
transmission, we will be using 80% of the original data to train the
model and 20% to test it. Then, we will evaluate its performance using
the Root Mean Squared Error
(RMSE)\index{Root Mean Squared Error (RMSE)}.

The `{caret}` package provides a unified interface for training and
evaluating machine learning models in R (@sec-07-packages). It supports
a wide range of algorithms and provides tools for hyperparameter
tuning\index{Hyperparameter tuning},
cross-validation\index{Cross Validation}, and model evaluation.

```{r}
#| eval: true
# Load libraries and check for conflicts
library(caret)
conflicted::conflicts_prefer(dplyr::lag)
```

#### Feature Engineering

Create lagged variables to capture the temporal
dynamics\index{Temporal dynamics} of malaria transmission. For time
series\index{Time series} forecasting, it’s helpful to create variables,
which represent past values of the target variable. This allows the
model to learn trends from previous values.

```{r}
# Create lagged features (e.g., previous year's cases)
nigeria_cases <- nigeria_cases %>%
  arrange(year) %>%
  mutate(lag_1 = lag(positive, 1),
         lag_2 = lag(positive, 2),
         lag_3 = lag(positive, 3)) %>%
  drop_na()

head(nigeria_cases)
```

Create partition for training and testing data using the
`createDataPartition()` function from the `caret` package.

```{r}
set.seed(123)
train_index <- createDataPartition(nigeria_cases$positive, 
                                   p = 0.8, list = FALSE)
train <- nigeria_cases[train_index, ]
test <- nigeria_cases[-train_index, ]
```

#### Model Selection and Training (Parameter Calibration)

Define the machine learning model specific for investigating Malaria
dynamics. In this case, we will use the time and number of lagged cases
to predict the number of infected cases. We will use the `train()`
function from the `{caret}` package to train the model. A suitable model
would be? List the models that can be used for this task:

-   Random Forest\index{Random
    Forest} ("rf")
-   Gradient Boosting\index{Gradient Boosting (GB)} ("gbm")
-   Support Vector Machines ("svm")\index{Support Vector Machines (SVM)}
-   Neural Networks\index{Neural Networks (NNet)} ("nnet") - etc.

In the `method` argument, specify the model to be used (e.g., "rf" for
Random Forest). The `trControl` argument specifies the cross-validation
method for hyperparameter tuning, and the `tuneGrid` argument defines
the hyperparameter grid for tuning.

```{r}
# Define the training control with 5-fold cross-validation
train_control <- trainControl(method = "cv", 
                              number = 5)

# Define a tuning grid for mtry 
# (number of variables sampled at each split)
tuning_grid <- expand.grid(mtry = c(1, 2, 3))

# Train the Random Forest model
set.seed(123)
rf_model <- train(
  positive ~ lag_1 + lag_2 + lag_3,
  data = train,
  method = "rf",
  trControl = train_control,
  tuneGrid = tuning_grid,
  ntree = 500)  # Number of trees 
```

The parameter calibration\index{Parameter calibration} (hyperparameter
tuning\index{Hyperparameter tuning}) is performed by the `train`
function automatically using cross-validation to select the optimal
hyperparameters for the model.

#### Model Evaluation

Testing the model on observed data will show how the model performs in
estimating the number of infected cases.

```{r}
# Predict on the test data
nigeria_cases$pred <- predict(rf_model, 
                              newdata = nigeria_cases)
```

```{r}
#| label: fig-rf-malaria-predictions
#| fig-cap: "The plot shows the observed vs predicted malaria positive cases in Nigeria made with a Random Forest model. The blue line represents the observed cases, while the red line represents the predicted cases. The model's predictions are not closely aligned with the observed values, indicating potential limitations in capturing the dynamics of malaria transmission."
#| fig-alt: "Line plot showing the observed vs predicted malaria positive cases in Nigeria."
#| fig-width: 6
#| fig-height: 4
nigeria_cases %>%
 ggplot(aes(x = year, y = positive)) +
  geom_line(color = "navy", linetype = "solid") +
  geom_line(aes(y = pred), 
            color = "darkred", linetype = "dashed") +
  labs(title = "Malaria Positive Cases in Nigeria",
       subtitle = "Observed vs Predicted - Random Forest",
       x = "Year",
       y = "Number of Positive Cases") +
  theme_minimal()
```

To evaluate the model's ability in predicting future trends, we predict
the number of infected cases on the test data and calculate the Root
Mean Squared Error (RMSE)\index{Root
Mean Squared Error (RMSE)}. Test data are a subset of the original data
that were not used for training the model and provide an independent
evaluation of the model’s performance.

```{r}
# Predict on the test data
test$pred <- predict(rf_model, 
                     newdata = test)

# View predictions alongside actual values
head(test[, c("year", "positive", "pred")])
```

Visualize the predicted vs actual infected cases using a line plot to
compare the model's predictions with the actual data.

```{r}
# Reshape data for plotting
plot_data_rf <- test %>%
  select(year, positive, pred) %>%
  rename(Observed = positive, Predicted = pred) %>%
  pivot_longer(cols = c("Observed", "Predicted"), 
               names_to = "Type", values_to = "Cases")

plot_data_rf %>% head()    
```

```{r}
#| label: fig-rf-malaria-predictions2
#| fig-cap: "The plot shows the observed vs predicted malaria positive cases in Nigeria made with a Random Forest model. The blue line represents the observed cases, while the red line represents the predicted cases. The model's predictions are not closely aligned with the observed values, indicating potential limitations in capturing the dynamics of malaria transmission."
#| fig-alt: "Line plot showing the observed vs predicted malaria positive cases in Nigeria."
# Plot observed vs predicted cases
ggplot(plot_data_rf, aes(x = year, y = Cases, 
                      color = Type, linetype = Type)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_color_manual(values = c("navy", "darkred")) +
  labs(title = "Malaria Positive Cases in Nigeria",
       subtitle = "Observed vs Predicted - Random Forest",
       x = "Year",
       y = "Number of Positive Cases") +
  theme_minimal() +
  theme(legend.position = "top")
```

Evaluate the model's performance using RMSE:

```{r}
# Calculate RMSE
rmse <- sqrt(mean((test$positive - test$pred)^2))
cat("RMSE:", rmse, "\n")
```

A Root Mean Squared Error (RMSE) of 247.3936 indicates the average
difference between the predicted and actual values of infected cases.
Lower RMSE values indicate better model performance.

The output shows that the model’s predictions are not closely aligned
with the observed values, the predicted trend seems to be relatively
flat or declining, while the observed cases show significant
fluctuations. This mismatch suggests that the current model setup may
not be effectively capturing the dynamics of malaria transmission in
this dataset.

## Model Refinement

To improve the model's performance, we can refine the model by adjusting
parameters such as 'mtry' in Random Forest. For example, increase the
number of folds in cross-validation, adjust the tuning grid, or add
additional features to improve the model's performance.

```{r}
#| output: false
# Increase the number of folds in cross-validation
train_control <- trainControl(method = "cv", 
                              number = 10)  

# Scaling data in caret
rf_model2 <- train(
  positive ~ year + lag_1 + lag_2 + lag_3,
  data = train,
  method = "rf",
  trControl = train_control,
  tuneGrid = expand.grid(mtry = 2:4),
  preProcess = c("center", "scale"))
```

Predict future trends using the trained adjusted Random Forest model on
the test data.

```{r}
# Predict on the test data
test$pred_rf2 <- predict(rf_model2, 
                         newdata = test)
```

Visualize the predicted vs actual infected cases using a line plot to
compare the adjusted Random Forest model's predictions with the actual
data.

```{r}
# Reshape data for plotting
plot_data_xgb <- test %>%
  select(year, positive, pred_rf2) %>%
  rename(Observed = positive, Predicted = pred_rf2) %>%
  pivot_longer(cols = c("Observed", "Predicted"), 
               names_to = "Type", values_to = "Cases")

plot_data_xgb %>% head()
```

```{r}
#| label: fig-xgb-malaria-predictions
#| fig-cap: "The plot shows the observed vs predicted malaria positive cases in Nigeria made with an adjusted Random Forest model. The blue line represents the observed cases, while the red line represents the predicted cases. The model's predictions show a closer alignment with the observed values compared to the first Random Forest model, indicating improved performance."
#| fig-alt: "Line plot showing the observed vs predicted malaria positive cases in Nigeria using XGBoost."  
# Plot observed vs predicted cases
ggplot(plot_data_xgb, aes(x = year, y = Cases, 
                      color = Type, linetype = Type)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_color_manual(values = c("navy", "darkred")) +
  labs(title = "Malaria Cases in Nigeria",
       subtitle = "Observed vs Predicted - Adjusted Random Forest (XGBoost)",
       x = "Year",
       y = "Number of Positive Cases") +
  theme_minimal() +
  theme(legend.position = "top")
```

Evaluate the model's performance improvement using RMSE:

```{r}
# Calculate RMSE
rmse <- sqrt(mean((test$positive - test$pred_rf2)^2))
cat("RMSE:", rmse, "\n")
```

To visualise the improvement in model's performance, we can plot the
observed vs predicted for both models:

```{r}
#| label: fig-rf-xgb-malaria-predictions
#| fig-cap: "The plot shows the observed vs predicted malaria positive cases in Nigeria made with both Random Forest and XGBoost models. The blue line represents the observed cases, while the red line represents the predicted cases from the Random Forest and XGBoost models."
#| fig-alt: "Line plot showing the observed vs predicted malaria positive cases in Nigeria using Random Forest and XGBoost."
#| fig-width: 6
#| fig-height: 4
#| echo: false
plot_data_rf %>%
ggplot(aes(x = year, y = Cases, 
       color = Type, linetype = Type)) +
  geom_line() +
  geom_line(data = plot_data_xgb) +
  geom_point(size = 2) +
  scale_color_manual(values = c("navy", "darkred")) +
annotate("text", x = 2007.5, y = 140, size = 3,
         label = "Random Forest", color = "darkred") +
annotate("text", x = 2007.5, y = 250, size = 3,
         label = "XGBoost", color = "darkred") +
# draw a curved arrow from "Random Forest" to the point (2007.5, 140)
geom_curve(aes(x = 2006.5, y = 140, 
               xend = 2007, yend = 190), 
           curvature = -0.2, 
           arrow = arrow(length = unit(0.1, "inches"))) +
  geom_curve(aes(x = 2008, y = 250, 
               xend = 2008, yend = 190), 
           curvature = -0.2, 
           arrow = arrow(length = unit(0.1, "inches"))) +
  labs(title = "Malaria Cases in Nigeria",
       subtitle = "Observed vs Predicted - Random Forest vs XGBoost",
       x = "Year",
       y = "Number of Positive Cases") +
  theme_minimal() +
  theme(legend.position = "top")
```

In conclusion, the model's performance can be further improved by
adjusting hyperparameters, adding additional features, or using more
advanced algorithms. The iterative process of model refinement and
evaluation is essential for developing accurate predictions of malaria
transmission dynamics.

## Summary

In this chapter, we explored the dynamics of malaria transmission and
emphasized the importance of mapping outbreaks to guide public health
interventions. Using the `malariaAtlas` package, we downloaded malaria
data for Nigeria and visualized the spatial distribution of cases across
the country. Additionally, and demonstrated how machine learning
techniques can analyse historical data to predict future trends in
malaria transmission.