1
$\begingroup$

I wanted to forecast Dengue count for the given weather data covariates:

MeanT Mean temperature

MeanRH Mean relative humidity

MeanSM Mean soil moisture

MeanVP Mean vapor pressure

MeanWS Mean wind speed

MeanWG Mean wind gust

DI Dengue incidence

I tried using machine learning and stat models such as LSTM, 1D CNN, Zero-Inflated Negative Binomial, and Hurdle, but they did not work. Now, I am attempting to use GAM and MARS, but it seems like these two violate all the assumptions, such as normality and constant variance. However, they did a pretty good job in forecasting on the test data. In fact, the MARS model had lower RMSE and MAE scores. There is another issue here: both models suffer from overfitting. I am not sure how to address these issues. My goal is to find important variables and make forecast. I appreciate your suggestions!

# List of packages
packages <- c("readxl", "dplyr", "ggplot2", "mgcv", "earth", "forecast", "Metrics", "readr", "curl")

# Install missing packages
for (pkg in packages) {
  if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
  }
}

# Load packages
for (pkg in packages) {
  library(pkg, character.only = TRUE)
}


url <- "https://gist.githubusercontent.com/JishanAhmed2019/4f5e915a4b490e1dfeb1ba4defcc85cf/raw/4bcf2f76d54f532f9287a8d40b46125b5a684a92/seasons_data.xlsx"
dengue <- read.csv(url, sep="\t", header=T)


dengue <- dengue %>% 
  mutate(
    Date = as.Date(Date),
    YearDay = as.numeric(format(Date, "%j")),
    Year = as.numeric(format(Date, "%Y")),
    # Ensure 'season' is a factor with appropriate levels
    season = factor(season, levels = c("Summer", "Rainy season", "Autumn", "Late autumn", "Winter", "Spring"))
  )

# Splitting the data
split_date <- as.Date("2023-08-22")  # Adjust accordingly
dengue_train <- dengue %>% filter(Date < split_date)
dengue_test <- dengue %>% filter(Date >= split_date)

# Fit models
# Note: Ensure 's()' terms are appropriate and 'k' is not set too high for the number of unique observations
gam_mod <- gam(Dhaka_DI ~ s(MeanT) + s(MeanRH) + s(MeanSM) + s(MeanVP) + s(MeanWS) + s(MeanWG) + s(Month, bs = "cc", k = 12)+
                        s(YearDay, k = 12, bs="cc") + season, 
               data = dengue_train, method = "REML")     

mars_mod <- earth(Dhaka_DI ~ MeanT + MeanRH + MeanSM + MeanVP + MeanWS + MeanWG +Month+YearDay + season, 
                  data = dengue_train)

# Diagnostics for GAM - check for patterns or outliers
par(mfrow = c(2, 2))
plot(gam_mod)

# Predictions
dengue_train$Pred_GAM <- predict(gam_mod, newdata = dengue_train)
dengue_test$Pred_GAM <- predict(gam_mod, newdata = dengue_test)
dengue_train$Pred_MARS <- predict(mars_mod, newdata = dengue_train)
dengue_test$Pred_MARS <- predict(mars_mod, newdata = dengue_test)

# Performance Metrics
evaluate_metrics <- function(actual, predicted) {
  list(
    MAPE = mape(actual, predicted),
    RMSE = rmse(actual, predicted),
    MAE = mae(actual, predicted)
  )
}

train_metrics_gam <- evaluate_metrics(dengue_train$Dhaka_DI, dengue_train$Pred_GAM)
test_metrics_gam <- evaluate_metrics(dengue_test$Dhaka_DI, dengue_test$Pred_GAM)
train_metrics_mars <- evaluate_metrics(dengue_train$Dhaka_DI, dengue_train$Pred_MARS)
test_metrics_mars <- evaluate_metrics(dengue_test$Dhaka_DI, dengue_test$Pred_MARS)

# Print metrics
cat("Training Metrics (GAM):\n")
print(train_metrics_gam)
cat("\nTest Metrics (GAM):\n")
print(test_metrics_gam)
cat("\nTraining Metrics (MARS):\n")
print(train_metrics_mars)
cat("\nTest Metrics (MARS):\n")
print(test_metrics_mars)

# Visualization
ggplot() +
  geom_line(data = dengue_train, aes(x = Date, y = Dhaka_DI, color = "Actual Train"), size = 1) +
  geom_line(data = dengue_train, aes(x = Date, y = Pred_GAM, color = "GAM Train"), size = 1, linetype = "dashed") +
  geom_line(data = dengue_train, aes(x = Date, y = Pred_MARS, color = "MARS Train"), size = 1, linetype = "dotted") +
  geom_line(data = dengue_test, aes(x = Date, y = Dhaka_DI, color = "Actual Test"), size = 1) +
  geom_line(data = dengue_test, aes(x = Date, y = Pred_GAM, color = "GAM Test"), size = 1, linetype = "dashed") +
  geom_line(data = dengue_test, aes(x = Date, y = Pred_MARS, color = "MARS Test"), size = 1, linetype = "dotted") +
  labs(title = "Dengue Forecast with GAM and MARS Models", x = "Date", y = "Dengue Count") +
  theme_minimal() +
  scale_color_manual(values = c("Actual Train" = "black", "GAM Train" = "blue", "MARS Train" = "green", 
                                "Actual Test" = "grey", "GAM Test" = "red", "MARS Test" = "orange")) +
  guides(color = guide_legend(title = "Legend"))


summary(gam_mod)
summary(mars_mod)

plot(mars_mod)

plot(gam_mod)

gam.check(gam_mod)
$\endgroup$
2
  • 1
    $\begingroup$ Your GAM assumes $Y$ is conditionally distributed Gaussian, so it is not surprising that a response that can only take values in the set of Natural numbers (non negative integers) violates that assumption. But you aren't restricted to modelling conditionally Gaussian distributed responses; you can use them like GLMs, but also the mgcv 📦 implements methods that extend REML smoothness selection to some non-exponential family distributions. So I would start with family = nb() in your model to fit a GAM with $Y$ assumed to be conditionally distributed negative binomial and proceed from there. $\endgroup$ Commented Mar 8 at 8:21
  • $\begingroup$ @GavinSimpson I attempted to use family = nb() and ZiP, but they showed the worst performance! In fact, I have watched your YouTube videos multiple times to resolve my issues, but I'm not sure what went wrong! $\endgroup$
    – Quantam
    Commented Mar 8 at 8:43

1 Answer 1

2
$\begingroup$

Thanks for the reproducible example. I agree with Gavin that a Gaussian observation model is not really appropriate here. But you have other complexities to deal with as well. I'll use an R package, {mvgam}, to illustrate. This package was designed for exactly these kinds of problems, i.e. using GAMs to analyse and forecast time series (installation instructions and examples are found here: https://nicholasjclark.github.io/mvgam/).

First I load the data, add a series and a time indicator and plot features of the series:

url <- "https://gist.githubusercontent.com/JishanAhmed2019/4f5e915a4b490e1dfeb1ba4defcc85cf/raw/4bcf2f76d54f532f9287a8d40b46125b5a684a92/seasons_data.xlsx"
dengue <- read.csv(url, sep="\t", header=T)


dengue <- dengue %>% 
  mutate(
    Date = as.Date(Date),
    YearDay = as.numeric(format(Date, "%j")),
    Year = as.numeric(format(Date, "%Y")),
    # Ensure 'season' is a factor with appropriate levels
    season = factor(season, levels = c("Summer", "Rainy season", "Autumn", "Late autumn", "Winter", "Spring"))
  )

# Load mvgam; add a 'series' indicator and look at the data
library(mvgam)
dengue %>%
  mutate(time = dplyr::row_number(),
         series = as.factor('Dhaka')) -> dengue

# Huge autocorrelation and many zeros to think about
plot_mvgam_series(data = dengue, 
                  y = 'Dhaka_DI')

Time series features produced with mvgam

As Gavin mentioned, your response is non-negative integers and there are many zeros. A Gaussian observation model won't respect these constraints. Next I split the data into training and testing as you have done, and fit a model that only includes the seasonal effects as well as a latent dynamic process (as an AR1 process):

# Splitting the data
split_date <- as.Date("2023-08-22")  # Adjust accordingly
dengue_train <- dengue %>% filter(Date < split_date)
dengue_test <- dengue %>% filter(Date >= split_date)

# Fit a simple model with Poisson observations over a cyclic YearDay 
# smooth, a parametric season effect and a latent AR1 process
gam_mod <- mvgam(Dhaka_DI ~ s(YearDay, k = 12, bs = "cc") + 
                   season, 
                 trend_model = AR(),
               data = dengue_train, 
               family = poisson()) 

By default, mvgam() will estimate model parameters using Stan for full Bayesian inference. This takes a while to fit (120 seconds on my i9 processor), which is a bad sign given the simplicity of the model. The summary doesn't tell us too much, apart from showing that the seasonal spline is a bit difficult to estimate (look at the low effective sample sizes of the spline coefficients, and the warnings about treedepths):

summary(gam_mod)
GAM formula:
Dhaka_DI ~ s(YearDay, k = 12, bs = "cc") + season

Family:
poisson

Link function:
log

Trend model:
AR()

N series:
1 

N timepoints:
598 

Status:
Fitted using Stan 
4 chains, each with iter = 1000; warmup = 500; thin = 1 
Total post-warmup draws = 2000


GAM coefficient (beta) estimates:
                    2.5%    50% 97.5% Rhat n_eff
(Intercept)         2.30  3.100  3.80 1.01   528
seasonRainy season -1.30 -0.530  0.21 1.00  1099
seasonAutumn       -1.10 -0.097  0.80 1.00  1295
seasonLate autumn  -0.75  0.470  1.60 1.00  1255
seasonWinter       -1.80 -0.490  0.72 1.00  1275
seasonSpring       -2.50 -1.300 -0.27 1.00   830
s(YearDay).1       -2.80 -1.800 -0.83 1.02   298
s(YearDay).2       -3.20 -2.200 -1.10 1.02   388
s(YearDay).3       -3.00 -1.900 -0.77 1.02   517
s(YearDay).4       -1.80 -0.760  0.38 1.01   561
s(YearDay).5       -0.12  1.000  2.20 1.01   435
s(YearDay).6        1.10  2.400  3.60 1.02   325
s(YearDay).7        1.70  3.000  4.00 1.02   314
s(YearDay).8        1.50  2.800  4.10 1.00   374
s(YearDay).9        0.95  2.300  3.60 1.00   414
s(YearDay).10      -0.10  1.100  2.30 1.01   433

Approximate significance of GAM observation smooths:
            edf Chi.sq p-value   
s(YearDay) 7.46  16008  0.0078 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Latent trend parameter AR estimates:
         2.5%  50% 97.5% Rhat n_eff
ar1[1]   0.86 0.91  0.95 1.01   500
sigma[1] 0.50 0.55  0.59 1.01   768

Stan MCMC diagnostics:
n_eff / iter looks reasonable for all parameters
Rhat looks reasonable for all parameters
0 of 2000 iterations ended with a divergence (0%)
399 of 2000 iterations saturated the maximum tree depth of 12 (19.95%)
 *Run with max_treedepth set to a larger value to avoid saturation
E-FMI indicated no pathological behavior

Samples were drawn using NUTS(diag_e) at Sun Mar 10 9:56:02 PM 2024.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split MCMC chains
(at convergence, Rhat = 1)

Plotting the cyclic smooth doesn't really indicate any issues:

plot(gam_mod, type = 'smooths')

Cyclic smooth using bs = 'cc' with mvgam

But showing the conditional cyclic smooth against the observed data is more useful (these plots use support from the {marginaleffects} package:

plot_predictions(gam_mod, condition = c('YearDay'),
                 points = 0.6) +
  theme_classic()

conditional smooths from mvgam model, plotted with marginaleffects

Here you can see that the seasonality changes quite a lot between years, and the cyclic smooth is forced to make a compromise between the two years. As a result, it doesn't fit the data very well and the latent AR1 process is forced to compensate by being extremely flexible (almost behaving as a Random Walk):

plot(gam_mod, type = 'trend')

latent AR1 dynamic trend estimated with mvgam

The resulting forecast is therefore not great:

fc <- forecast(gam_mod, newdata = dengue_test)
plot(fc)

probabilistic forecasts from mvgam model

I suggest you look into the data a bit more and try to work out why seasonality changes so much from year to year. In the past I've had success doing this with distributed lags of climate variables, for example to predict time-varying seasonal curves of tick paralysis incidence in Australian dogs: https://www.youtube.com/watch?v=Yu5zD4WSrmA&list=PLzFHNoUxkCvtIGABuakH_T5CLVZPedaXQ. But I can't say whether this would work in your context

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.