Exploratory Analysis

Load packages

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(fpp3)

── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──
✔ tsibble     1.1.3     ✔ fable       0.3.3
✔ tsibbledata 0.4.1     ✔ fabletools  0.3.4
✔ feasts      0.3.1     
── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()    masks base::date()
✖ dplyr::filter()      masks stats::filter()
✖ tsibble::intersect() masks base::intersect()
✖ tsibble::interval()  masks lubridate::interval()
✖ dplyr::lag()         masks stats::lag()
✖ tsibble::setdiff()   masks base::setdiff()
✖ tsibble::union()     masks base::union()

library(ggthemes)
library(here)

here() starts at /Users/samkuhn/Documents/bc/time_series

library(corrplot)

corrplot 0.92 loaded

theme_set(theme_fivethirtyeight())

Load data

train_features <- readr::read_csv(here::here("final_project/data", "dengue_features_train.csv"))

Rows: 1456 Columns: 24
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (1): city
dbl  (22): year, weekofyear, ndvi_ne, ndvi_nw, ndvi_se, ndvi_sw, precipitati...
date  (1): week_start_date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

train_labels <- readr::read_csv(here::here("final_project/data", "dengue_labels_train.csv"))

Rows: 1456 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city
dbl (3): year, weekofyear, total_cases

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Data preprocessing

Combine year and week

train_labels <- train_labels |> 
  mutate(date = lubridate::make_date(year = year) + weeks(weekofyear))

train_features <- train_features |> 
  rename(date = "week_start_date")

Check missing values

check_na <- function(data) {
  data |> 
    summarise(across(everything(), ~ sum(is.na(.))))
}
check_na(train_labels)

# A tibble: 1 × 5
   city  year weekofyear total_cases  date
  <int> <int>      <int>       <int> <int>
1     0     0          0           0     0

check_na(train_features)

# A tibble: 1 × 24
   city  year weekofyear  date ndvi_ne ndvi_nw ndvi_se ndvi_sw
  <int> <int>      <int> <int>   <int>   <int>   <int>   <int>
1     0     0          0     0     194      52      22      22
# ℹ 16 more variables: precipitation_amt_mm <int>, reanalysis_air_temp_k <int>,
#   reanalysis_avg_temp_k <int>, reanalysis_dew_point_temp_k <int>,
#   reanalysis_max_air_temp_k <int>, reanalysis_min_air_temp_k <int>,
#   reanalysis_precip_amt_kg_per_m2 <int>,
#   reanalysis_relative_humidity_percent <int>,
#   reanalysis_sat_precip_amt_mm <int>,
#   reanalysis_specific_humidity_g_per_kg <int>, reanalysis_tdtr_k <int>, …

Impute mean

impute_missing_by_mean <- function(data) {
  # Loop through columns
  for (col in names(data)) {
    # Check if column has missing values
    if (any(is.na(data[[col]]))) {
      # Calculate mean excluding missing values
      col_mean <- mean(data[[col]], na.rm = TRUE)
      
      # Impute missing values with mean
      data[[col]][is.na(data[[col]])] <- col_mean
    }
  }
  
  return(data)
}
train_labels <- impute_missing_by_mean(train_labels)
train_features <- impute_missing_by_mean(train_features)

Viz

Total case count - SJ

train_labels |> 
  filter(city == "sj") |> 
  ggplot(aes(date, total_cases, color = "temp")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  scale_x_date(
    breaks = scales::breaks_width("2 years"),
    labels = scales::label_date_short()
  ) +
  annotate(
    geom = "text",
    x = as.Date("1994-10-15"),
    y = 475,
    label = "Peak case count of 461"
  ) +
  labs(
    title = "Total weekly Dengue Fever case count, SJ"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  )

ggsave(here::here("final_project/figures", "total_case_count_sj.png"))

Saving 7 x 5 in image

Total case count - IQ

train_labels |> 
  filter(
    city == "iq",
    date >= "2000-01-01"
  ) |> 
  ggplot(aes(date, total_cases, color = "temp")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  scale_x_date(
    breaks = scales::breaks_width("2 years"),
    labels = scales::label_date_short()
  ) +
  annotate(
    geom = "text",
    x = as.Date("2005-10-15"),
    y = 125,
    label = "Peak case count of 116"
  ) +
  labs(
    title = "Total weekly Dengue Fever case count, IQ"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  )

ggsave(here::here("final_project/figures", "total_case_count_iq.png"))

Saving 7 x 5 in image

Correlation plot of numeric vectors

numeric_vars <- train_features |> 
  select(-c("city", "date")) |> 
  as.matrix()

corr_matrrix <- cor(numeric_vars)

png(here::here("final_project/figures", "corrplot.png"))
corrplot(corr_matrrix, method = "color", order = "alphabet")
dev.off()

quartz_off_screen 
                2 

Time series of all numeric vectors

plot_time_series <- function(col) {
  train_features |> 
    ggplot(aes(week_start_date, {{col}})) +
    geom_line(linewidth = .75)
}

Train/test splits

train_features <- train_features |> 
  select(-c("year", "weekofyear"))

train_labels <- train_labels |> 
  select(-c("year", "weekofyear"))

train <- train_features |> 
  left_join(train_labels) |> 
  mutate(date = yearweek(date)) |> 
  as_tsibble(index = "date", key = "city")

Joining with `by = join_by(city, date)`

Impute missing dates

train <- train |>
  group_by_key() |> 
  pivot_longer(where(is.numeric), names_to = "features") |> 
  fill_gaps() |> 
  tidyr::fill(value) |> 
  pivot_wider(names_from = "features")

Inspect seasonality

stl_decomp <- train |> 
  model(stl = STL(total_cases ~ trend() + season(52)))

stl_decomp |> 
  filter(city == "sj") |> 
  components() |> 
  autoplot() +
  labs(
    title = "STL decomposition, SJ"
  )

ggsave(here::here("final_project/figures", "stl_decomp_sj.png"))

Saving 7 x 5 in image

stl_decomp |> 
  filter(city == "iq") |> 
  components() |> 
  autoplot() +
  labs(
    title = "STL decomposition, IQ"
  )

ggsave(here::here("final_project/figures", "stl_decomp_iq.png"))

Saving 7 x 5 in image

SJ

sj_df <- train |> 
  filter(city == "sj")

sj_train <- sj_df |> 
  filter_index(.~ "2004-09-09")

sj_test <- sj_df |> 
  filter_index("2004-09-16"~.)

IQ

iq_df <- train |> 
  filter(city == "iq")

iq_train <- iq_df |> 
  filter_index(.~ "2008-06-24")

iq_test <- iq_df |> 
  filter_index("2008-07-01"~.)

Calculate box-cox transformation lambda

lambda <- train |>
  features(features = guerrero) |>
  pull(lambda_guerrero)

Feature variable not specified, automatically selected `.var = ndvi_ne`

Fit models

sj_fit <- sj_train |> 
  model(
    lm = TSLM(total_cases ~ ndvi_ne:station_precip_mm),
    spline_lm = TSLM(total_cases ~ trend(knots = c(1994, 1998)) + ndvi_ne:station_precip_mm),
    nnet = NNETAR(box_cox(total_cases, lambda[1]) ~ AR(p = NULL, period = 52))
  ) |> 
  mutate(ensemble = (lm + spline_lm + nnet) / 3)

iq_fit <- iq_train |> 
  model(
    lm = TSLM(total_cases ~ ndvi_ne:station_precip_mm),
    spline_lm = TSLM(total_cases ~ trend(knots = c(2000, 2004)) + ndvi_ne:station_precip_mm),
    nnet = NNETAR(box_cox(total_cases, lambda[2]) ~ AR(p = NULL, period = 52))
  ) |> 
  mutate(ensemble = (lm + spline_lm + nnet) / 3)

Forecast against test sets

sj_fc <- sj_fit |> 
  forecast(sj_test)

Warning: There were 2 warnings in `mutate()`.
The first warning was:
ℹ In argument: `spline_lm = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleading
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

iq_fc <- iq_fit |> 
  forecast(iq_test)

Plot forecasts against test set

sj_fc |> 
  autoplot(
    bind_rows(sj_train, sj_test),
    level = NULL
  ) +
  labs(
    title = "Forecasted weekly case count for SJ",
    subtitle = "Colored lines are forecasted values from training data"
  )

`mutate_if()` ignored the following grouping variables:
• Column `city`

ggsave(here::here("final_project/figures", "train_fc_sj.png"))

Saving 7 x 5 in image

iq_fc |> 
  autoplot(
    bind_rows(iq_train, iq_test),
    level = NULL
  ) +
  labs(
    title = "Forecasted weekly case count for IQ",
    subtitle = "Colored lines are forecasted values from training data"
  )

`mutate_if()` ignored the following grouping variables:
• Column `city`

ggsave(here::here("final_project/figures", "train_fc_iq.png"))

Saving 7 x 5 in image

Accuracy

accuracy(sj_fc, train) |> arrange(MAE)

# A tibble: 4 × 11
  .model    city  .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
  <chr>     <chr> <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 spline_lm sj    Test   11.4   33.5  18.4  -Inf   Inf 0.464 0.427 0.911
2 ensemble  sj    Test   -5.25  31.5  24.1  -Inf   Inf 0.607 0.402 0.915
3 lm        sj    Test  -10.5   33.1  27.3  -Inf   Inf 0.689 0.422 0.888
4 nnet      sj    Test  -16.5   35.1  30.1  -Inf   Inf 0.760 0.448 0.931

accuracy(iq_fc, train) |> arrange(MAE)

# A tibble: 4 × 11
  .model    city  .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
  <chr>     <chr> <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lm        iq    Test    2.87  11.6  7.10  -Inf   Inf 0.752 0.746 0.847
2 nnet      iq    Test   -4.20  12.6 10.0   -Inf   Inf 1.06  0.808 0.856
3 ensemble  iq    Test   -5.91  15.1 12.4   -Inf   Inf 1.32  0.966 0.894
4 spline_lm iq    Test  -16.5   26.7 22.8   -Inf   Inf 2.42  1.71  0.943

Residual analysis

sj_fc |> 
  filter(.model == "spline_lm") |> 
  gg_tsdisplay() |> 
  labs(
    title = "Residuals, spline regression for SJ"
  )

Plot variable not specified, automatically selected `y = .mean`

Warning: The output of `fortify(<fable>)` has changed to better suit usage with the ggdist package.
If you're using it to extract intervals, consider using `hilo()` to compute intervals, and `unpack_hilo()` to obtain values.

[[1]]

$title
[1] "Residuals, spline regression for SJ"

attr(,"class")
[1] "labels"

iq_fc |> 
  filter(.model == "lm") |> 
  gg_tsdisplay()

Plot variable not specified, automatically selected `y = .mean`

Warning: The output of `fortify(<fable>)` has changed to better suit usage with the ggdist package.
If you're using it to extract intervals, consider using `hilo()` to compute intervals, and `unpack_hilo()` to obtain values.

Forecast with best model

test <- readr::read_csv(here::here("final_project/data", "dengue_features_test.csv"))

Rows: 416 Columns: 24
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (1): city
dbl  (22): year, weekofyear, ndvi_ne, ndvi_nw, ndvi_se, ndvi_sw, precipitati...
date  (1): week_start_date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

test <- test |> 
  select(-c("year", "weekofyear")) |> 
  mutate(date = yearweek(week_start_date)) |> 
  as_tsibble(index = date, key = city)
test <- impute_missing_by_mean(test)

sj_test <- test |> 
  filter(city == "sj")

iq_test <- test |> 
  filter(city == "iq")

sj_best <- sj_train |> 
  model(
    spline_lm = TSLM(total_cases ~ trend(knots = c(1994, 1998)) + season() + ndvi_ne:station_precip_mm)
  )

iq_best <- iq_train |> 
  model(
    lm = TSLM(total_cases ~ ndvi_ne:station_precip_mm)
  )

sj_preds <- sj_best |> 
  forecast(new_data = sj_test)

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `spline_lm = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleading

iq_preds <- iq_best |> 
  forecast(new_data = iq_test)

preds_df <- 
  bind_rows(sj_preds, iq_preds)

Convert numeric to integer case count

preds_df$.mean <- as.integer(preds_df$.mean)

Drop all exogenous variables, convert yearweek to isoweek + year

preds_df <- preds_df |> 
  tibble() |> 
  mutate(
    year = lubridate::isoyear(date),
    weekofyear = lubridate::isoweek(date)
  ) |> 
  select(city, year, weekofyear, total_cases = ".mean") |> 
  arrange(year)

Write to CSV

readr::write_csv(preds_df, here::here("final_project/data", "case_predictions.csv"))