add residual diagnostics and summary functions

DESCRIPTION

@@ -10,9 +10,9 @@ Description: This package is for fitting Bayesian Generalised Additive Models to
              The primary purpose of the package is to build a version of dynamic factor models that incorporates 
              the flexibility of GAMs in the linear predictor and latent dynamic trend components that are useful 
              for time series analysis and forecasting. Estimation is performed using Markov Chain Monte Carlo (MCMC) by
-             the Gibbs sampling software JAGS. The package also includes utilities for online updating of 
-             forecast distributions with a recursive particle filter that uses sequential importance sampling to 
-             assimilate new observations as they become available.
+             the Gibbs sampling software JAGS (requires version 4.3.0 or higher). The package also includes 
+             utilities for online updating of  forecast distributions with a recursive particle filter that 
+             uses sequential importance sampling to assimilate new observations as they become available.
 Maintainer: Nicholas J Clark <[email protected]>
 License: MIT + file LICENSE
 Depends: R (>= 3.6.0), mgcv, rjags, parallel

NAMESPACE

@@ -11,10 +11,12 @@ export(pfilter_mvgam_init)
 export(pfilter_mvgam_online)
 export(pfilter_mvgam_smooth)
 export(plot_mvgam_fc)
+export(plot_mvgam_resids)
 export(plot_mvgam_smooth)
 export(plot_mvgam_trend)
 export(plot_mvgam_uncertainty)
 export(roll_eval_mvgam)
 export(series_to_mvgam)
 export(sim_mvgam)
+export(summary_mvgam)
 importFrom(magrittr,"%>%")

NEON_manuscript/amb_hyp_tests.R

@@ -47,9 +47,9 @@ hyp3 = y ~
   s(season, by = series, m = 1, k = 8) - 1
 
 # Fit each hypothesis
-n.burnin = 100000
-n.iter = 5000
-thin = 5
+n.burnin = 1000
+n.iter = 1000
+thin = 1
 
 # Use default priors for latent drift terms and for smooth penalties
 fit_null <- fit_mvgam(data_train = all_data$data_train,
@@ -315,22 +315,18 @@ for(i in seq_along(levels(all_data$data_train$series))){
       mai = c(.38,.4,.45,.05),
       mgp = c(2, 1, 0))
   plot_mvgam_fc(object = fit_null$out_gam_mod, data_test = all_data$data_test,
-                data_train = all_data$data_train,
                 series = i)
   mtext(side = 3, 'Null', cex = 1,
         line = 1.55)
   plot_mvgam_fc(object = fit_hyp1$out_gam_mod, data_test = all_data$data_test,
-                data_train = all_data$data_train,
                 series = i)
   mtext(side = 3, 'Hyp1', cex = 1,
         line = 1.35)
   plot_mvgam_fc(object = fit_hyp2$out_gam_mod, data_test = all_data$data_test,
-                data_train = all_data$data_train,
                 series = i)
   mtext(side = 3, 'Hyp2', cex = 1,
         line = 1.55)
   plot_mvgam_fc(object = fit_hyp3$out_gam_mod, data_test = all_data$data_test,
-                data_train = all_data$data_train,
                 series = i)
   mtext(side = 3, 'Hyp3', cex = 1,
         line = 1.55)

NEON_manuscript/ixodes_hyp_tests.R

@@ -186,7 +186,7 @@ ggplot(plot_dat_ranks,
         panel.grid.minor = element_blank()) -> plot2
 
 dir.create('NEON_manuscript/Figures', recursive = T, showWarnings = F)
-pdf('NEON_manuscript/Figures/Fig3_Ixodes_performances.pdf', width = 6.25, height = 5)
+pdf('NEON_manuscript/Figures/Fig4_Ixodes_performances.pdf', width = 6.25, height = 5)
 cowplot::plot_grid(plot1, plot2, ncol = 1)
 dev.off()
 

NEON_manuscript/neon_analysis.R

@@ -10,7 +10,7 @@ library(mvgam)
 all_data <- prep_neon_data(species = 'Ixodes_scapularis', split_prop = 0.8)
 
 # Make plots for manuscript
-pdf('NEON_manuscript/Figures/Fig4_ixodes_example.pdf', width = 6.25, height = 5.85)
+pdf('NEON_manuscript/Figures/Fig5_ixodes_example.pdf', width = 6.25, height = 5.85)
 par(mfrow = c(2, 2),
     mgp = c(2.5, 1, 0),
     mai = c(0.6, 0.6, 0.2, 0.2))
@@ -39,7 +39,7 @@ dev.off()
 # Plot some of the varying seasonal shapes from the Amblyomma model
 load('NEON_manuscript/Results/amb_models.rda')
 all_data <- prep_neon_data(species = 'Ambloyomma_americanum', split_prop = 0.8)
-pdf('NEON_manuscript/Figures/Fig5_amb_seasonalities.pdf', width = 6.25, height = 5.85)
+pdf('NEON_manuscript/Figures/Fig6_amb_seasonalities.pdf', width = 6.25, height = 5.85)
 par(mfrow = c(2, 2),
     mgp = c(2.5, 1, 0),
     mai = c(0.6, 0.6, 0.2, 0.2))
@@ -58,14 +58,13 @@ for(i in c(1, 6, 9, 15)){
 dev.off()
 
 # Plot some of the different uncertainty decompositions for Amblyomma model
-pdf('NEON_manuscript/Figures/Fig6_amb_uncertainties.pdf', width = 6.25, height = 5.85)
+pdf('NEON_manuscript/Figures/Fig7_amb_uncertainties.pdf', width = 6.25, height = 5.85)
 par(mfrow = c(2, 2),
     mgp = c(2.5, 1, 0),
     mai = c(0.6, 0.6, 0.2, 0.2))
 for(i in c(4, 6, 9, 15)){
   plot_mvgam_uncertainty(fit_hyp3$out_gam_mod, series = i,
                           data_test = all_data$data_test,
-                          data_train = all_data$data_train,
                          legend_position = 'bottomleft')
 
 }

NEON_manuscript/neon_utility_functions.R

@@ -135,9 +135,8 @@ plot_mvgam_season = function(out_gam_mod, series, data_test, data_train,
                                                                                  levels(data_train$series)[series])])))
   Xp <- predict(out_gam_mod$mgcv_model, newdata = pred_dat, type = 'lpmatrix')
   betas <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'b')
-  gam_comps <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'gam_comp')[,series]
-  plot(gam_comps[1] * (Xp %*% betas[1,]) ~ pred_dat$season, ylim = range(gam_comps[1] * (Xp %*% betas[1,]) + 2.5,
-                                                                         gam_comps[1] * (Xp %*% betas[1,]) - 2.5),
+  plot((Xp %*% betas[1,]) ~ pred_dat$season, ylim = range((Xp %*% betas[1,]) + 2.5,
+                                                                         (Xp %*% betas[1,]) - 2.5),
 
        col = rgb(150, 0, 0, max = 255, alpha = 10), type = 'l',
        ylab = paste0('F(season) for ', levels(data_train$series)[series]),
@@ -146,9 +145,9 @@ plot_mvgam_season = function(out_gam_mod, series, data_test, data_train,
   rect(xleft = 4, xright = 18, ybottom = -100, ytop = 100, col = 'gray90',
        border = NA)
   abline(h = 0, lwd=1)
-  text(x=11,y=max(gam_comps[1] * (Xp %*% betas[1,]) + 2.3), labels = 'Peak tick season')
+  text(x=11,y=max((Xp %*% betas[1,]) + 2.3), labels = 'Peak tick season')
   for(i in 2:1000){
-    lines(gam_comps[i] * (Xp %*% betas[i,]) ~ pred_dat$season,
+    lines((Xp %*% betas[i,]) ~ pred_dat$season,
 
           col = rgb(150, 0, 0, max = 255, alpha = 10))
   }
@@ -173,10 +172,9 @@ plot_mvgam_gdd = function(out_gam_mod, series, data_test, data_train,
                                                                                  levels(data_train$series)[series])])))
   Xp <- predict(out_gam_mod$mgcv_model, newdata = pred_dat, type = 'lpmatrix')
   betas <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'b')
-  gam_comps <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'gam_comp')[,series]
   preds <- matrix(NA, nrow = 1000, ncol = length(pred_dat$cum_gdd))
   for(i in 1:1000){
-    preds[i,] <- rnbinom(length(pred_dat$cum_gdd), mu = exp(gam_comps[i] * (Xp %*% betas[i,])),
+    preds[i,] <- rnbinom(length(pred_dat$cum_gdd), mu = exp((Xp %*% betas[i,])),
                          size = MCMCvis::MCMCsummary(out_gam_mod$jags_output, 'r')$mean)
   }
   int <- apply(preds,

NEON_manuscript/portal_example.R

@@ -20,23 +20,24 @@ ds_resids = function(truth, fitted, size){
 
 # Convert PP catches to timeseries and then to mvgam format
 catches <- portal_dat$PP
-plot(catches, type = 'l')
 
 # Keep only data from the year 2004 onward for this example to make the series more manageable
 series <- xts::as.xts(subset(ts(catches, start = c(1979, 9), frequency = 12),
                              start = c(293)))
 head(series)
 colnames(series) <- c('PP')
+plot(series)
 PP_data <- series_to_mvgam(series, freq = 12, train_prop = 0.9)
 rm(series, portal_dat)
 
-# Fit standard GAM to the training data using an AR3 model (the maximum order that mvgam allows) for the residuals.
-# For a GAMM model, we cannot use negative binomial so will stick with Poisson (even though)
-# a negative binomial is likely more appropriate here. We add a 'time' index for the
-# form of the residual autocorrelation
+# Fit standard GAM to the training data using an AR3 model for the residuals.
+# For a GAMM model, we cannot use negative binomial so will stick with Poisson (even though
+# a negative binomial is likely more appropriate here). We add a 'time' index for the
+# form of the residual autocorrelation and use a 1st derivative penalty for the year effect
+# to prevent linear extrapolation of the yearly trend when forecasting
 PP_data$data_train$time <- seq(1:NROW(PP_data$data_train))
 gam_mod <- gamm(y ~ s(season, bs = 'cc', k = 8) +
-                 s(year, bs = 'gp', k = 4) + ti(season, year),
+                 s(year, bs = 'gp', k = 4, m = 1) + ti(season, year),
                data = PP_data$data_train,
                family = 'poisson', correlation = corARMA(form = ~ time, p = 3))
 summary(gam_mod$gam)
@@ -51,11 +52,14 @@ gam_resids <- ds_resids(truth = gam_mod$gam$y,
                         size = 1000)
 acf(gam_resids, na.action = na.pass)
 
+# Standardized residuals also show autocorrelation
+acf(residuals(gam_mod$lme,type="normalized"),main="standardized residual ACF")
+
 # Which kind of model might be suitable for the trend?
-plot(gam_resids, type = 'l')
+plot(residuals(gam_mod$lme,type="normalized"), type = 'l')
 forecast::auto.arima(gam_resids)
 
-# DS residuals show strong autocorrelation. Do Pearson residuals?
+# DS and standardized residuals both still show strong autocorrelation. Do Pearson residuals?
 gam_pearson_resids <- (gam_mod$gam$y - fv)/sqrt(fv)
 acf(gam_pearson_resids)
 forecast::auto.arima(gam_pearson_resids)
@@ -65,25 +69,29 @@ forecast::auto.arima(gam_pearson_resids)
 # which also highlights existing autocorrelation
 gam.check(gam_mod$gam)
 
-# Equivalent model using mvgam's dynamic trend component.
+# Now a model using mvgam
 # These models generally require at least 5 - 15K iterations
 # of adaptation to tune the samplers efficiently, so it can be a little slow. Note that
 # the model will predict for the data_test observations by considering the outcomes as
-# missing for these observations
+# missing for these observations. We also incorporate prior knowledge about upper bounds on
+# this series, which exist due to limits on the number of traps used in each trapping session
 PP_data$data_train$time <- NULL
 df_gam_mod <- mvjagam(data_train = PP_data$data_train,
                       data_test = PP_data$data_test,
                formula = y ~ s(season, bs = c('cc'), k = 8) +
-                 s(year, bs = 'gp', k = 4) + ti(season, year),
+                 s(year, bs = 'gp', k = 4, m = 1) + ti(season, year),
                trend_model = 'AR3',
-               family = 'poisson',
+               family = 'nb',
                n.burnin = 10000,
-               n.iter = 2000,
-               thin = 2,
-               auto_update = F)
+               n.iter = 1000,
+               thin = 1,
+               auto_update = F,
+               upper_bounds = 100)
+summary_mvgam(df_gam_mod)
+summary(df_gam_mod$mgcv_model)
 
 # ACF of DS residuals shows no remaining autocorrelation
-acf(df_gam_mod$resids$PP, na.action = na.pass)
+plot_mvgam_resids(df_gam_mod, 1)
 
 # Smooth function plots from each model
 par(mfrow=c(1,2))
@@ -123,8 +131,10 @@ polygon(c(seq(1:(NCOL(cred_ints))), rev(seq(1:NCOL(cred_ints)))),
 lines(cred_ints[2,], col = rgb(150, 0, 0, max = 255), lwd = 2, lty = 'dashed')
 points(PP_data$data_test$y[1:12], pch = 16)
 
-# The mgcv model overpredicts for the out of sample test set, though I don't know
-# how to incorporate the autocorrelation component into this forecast
+# The mgcv model overpredicts for the out of sample test set and has unreasonably narrow
+# prediction intervals, though admittedly I don't know how to easily
+# incorporate the autocorrelation component into this forecast.
+
 # Forecast for the dynamic gam model on the same y-axis scale
 fits <- MCMCvis::MCMCchains(df_gam_mod$jags_output, 'ypred')
 fits <- fits[,(NROW(df_gam_mod$obs_data)+1):(NROW(df_gam_mod$obs_data)+13)]
@@ -154,21 +164,23 @@ par(mfrow = c(1,1))
 
 ## Online forecasting from the dynamic gam
 # Initiate particle filter by assimilating the next observation in data_test
-pfilter_mvgam_init(object = df_gam_mod, n_particles = 20000, n_cores = 3,
+pfilter_mvgam_init(object = df_gam_mod, n_particles = 40000, n_cores = 4,
                    data_assim = PP_data$data_test)
 
-# Assimilate some of the next observations
-pfilter_mvgam_online(data_assim = PP_data$data_test[1:3,], n_cores = 3,
+# Assimilate some of the next observations in the series using the
+# kernel smoothing sequential monte carlo algorithm
+pfilter_mvgam_online(data_assim = PP_data$data_test[1:5,], n_cores = 4,
                      kernel_lambda = 1)
 
 # Forecast from particles using the covariate information in remaining data_test observations
-fc <- pfilter_mvgam_fc(file_path = 'pfilter', n_cores = 3,
+fc <- pfilter_mvgam_fc(file_path = 'pfilter', n_cores = 4,
                        data_test = PP_data$data_test, ylim = c(0, 100),
                        plot_legend = F)
 
-# Compare to original forecast
+# Compare the updated forecast to the original forecast
 par(mfrow=c(1,2))
 plot_mvgam_fc(df_gam_mod, 1, data_test = PP_data$data_test, ylim = c(0, 100))
 fc$PP()
 points(c(PP_data$data_train$y,
              PP_data$data_test$y), pch = 16)
+

NEON_manuscript/simulation_analysis.R

@@ -1,4 +1,105 @@
 #### Simulation analysis ####
+# Simulate a series with a nonlinear trend to visualise how a spline extrapolates
+library(mvgam)
+library(xts)
+library(forecast)
+data("AirPassengers")
+set.seed(1e8)
+dat <- floor(AirPassengers + cumsum(rnorm(length(AirPassengers),
+                                          sd = 35)))
+dat <- dat + abs(min(dat))
+series <- ts(dat, start = c(1949, 1), frequency = 12)
+fake_data <- series_to_mvgam(series, freq = 365, train_prop = 0.74)
+
+# Fit a GAM to the data using mgcv; use a wiggly Gaussian Process smooth for the
+# trend and a cyclic smooth for the seasonality
+gam_mod <- gam(y ~ s(year, k = 9, bs = 'tp') +
+                 s(season, bs = 'cc', k = 10) +
+                 ti(season, year),
+               data = fake_data$data_train,
+               family = nb(),
+               method = 'REML')
+
+# Calculate predictions for the mgcv model to inspect extrapolation behaviour
+Xp <- predict(gam_mod, newdata = rbind(fake_data$data_train,
+                                       fake_data$data_test), type = 'lpmatrix')
+vc <- vcov(gam_mod)
+sim <- MASS::mvrnorm(1000,
+                     mu = coef(gam_mod), Sigma = vc)
+dims_needed <- dim(exp(Xp %*% t(sim)))
+fits <- rnbinom(n = prod(dims_needed), mu = as.vector(exp(Xp %*% t(sim))),
+                size = gam_mod$family$getTheta(TRUE))
+fits <- t(matrix(fits, nrow = dims_needed[1], ncol = dims_needed[2]))
+
+# Plot the smooth function and forecast 95% and 68% HPD intervals
+pdf('NEON_manuscript/Figures/Fig1_extrapolation_example.pdf', width = 6.25, height = 5.85)
+par(mfrow = c(2, 2),
+    mgp = c(2.5, 1, 0),
+    mai = c(0.6, 0.6, 0.2, 0.2))
+plot(gam_mod, select = 1, shade = T, main = '2nd derivative penalty')
+cred_ints <- apply(fits, 2, function(x) hpd(x, 0.95))
+ylims <- c(0,
+           max(cred_ints) + 2)
+plot(cred_ints[3,] ~ seq(1:NCOL(cred_ints)), type = 'l',
+     col = rgb(1,0,0, alpha = 0),
+     ylim = ylims,
+     ylab = 'Predicted counts',
+     xlab = 'Time')
+polygon(c(seq(1:(NCOL(cred_ints))), rev(seq(1:NCOL(cred_ints)))),
+        c(cred_ints[1,],rev(cred_ints[3,])),
+        col = rgb(150, 0, 0, max = 255, alpha = 100), border = NA)
+cred_ints <- apply(fits, 2, function(x) hpd(x, 0.68))
+polygon(c(seq(1:(NCOL(cred_ints))), rev(seq(1:NCOL(cred_ints)))),
+        c(cred_ints[1,],rev(cred_ints[3,])),
+        col = rgb(150, 0, 0, max = 255, alpha = 180), border = NA)
+lines(cred_ints[2,], col = rgb(150, 0, 0, max = 255), lwd = 2, lty = 'dashed')
+lines(as.vector(series))
+abline(v = NROW(fake_data$data_train),
+       lty = 'dashed')
+text(x=NROW(fake_data$data_train) + 5,
+     y=3000, labels = 'Forecast horizon', srt = -90)
+
+# Repeat by shifting the training window slightly so see how the extrapolation of the
+# mgcv model is simply reacting to autocorrelation
+fake_data <- series_to_mvgam(series, freq = 365, train_prop = 0.74)
+gam_mod <- gam(y ~ s(year, k = 9, bs = 'tp', m = 1) +
+                 s(season, bs = 'cc', k = 10) +
+                 ti(season, year),
+               data = fake_data$data_train,
+               family = nb(),
+               method = 'REML')
+plot(gam_mod, select = 1, shade = T, main = '1st derivative penalty')
+Xp <- predict(gam_mod, newdata = rbind(fake_data$data_train,
+                                       fake_data$data_test), type = 'lpmatrix')
+vc <- vcov(gam_mod)
+sim <- MASS::mvrnorm(1000,
+                     mu = coef(gam_mod), Sigma = vc)
+dims_needed <- dim(exp(Xp %*% t(sim)))
+fits <- rnbinom(n = prod(dims_needed), mu = as.vector(exp(Xp %*% t(sim))),
+                size = gam_mod$family$getTheta(TRUE))
+fits <- t(matrix(fits, nrow = dims_needed[1], ncol = dims_needed[2]))
+cred_ints <- apply(fits, 2, function(x) hpd(x, 0.95))
+ylims <- c(0,
+           max(cred_ints) + 2)
+plot(cred_ints[3,] ~ seq(1:NCOL(cred_ints)), type = 'l',
+     col = rgb(1,0,0, alpha = 0),
+     ylim = ylims,
+     ylab = 'Predicted counts',
+     xlab = 'Time')
+polygon(c(seq(1:(NCOL(cred_ints))), rev(seq(1:NCOL(cred_ints)))),
+        c(cred_ints[1,],rev(cred_ints[3,])),
+        col = rgb(150, 0, 0, max = 255, alpha = 100), border = NA)
+cred_ints <- apply(fits, 2, function(x) hpd(x, 0.68))
+polygon(c(seq(1:(NCOL(cred_ints))), rev(seq(1:NCOL(cred_ints)))),
+        c(cred_ints[1,],rev(cred_ints[3,])),
+        col = rgb(150, 0, 0, max = 255, alpha = 180), border = NA)
+lines(cred_ints[2,], col = rgb(150, 0, 0, max = 255), lwd = 2, lty = 'dashed')
+lines(as.vector(series))
+abline(v = NROW(fake_data$data_train),
+       lty = 'dashed')
+dev.off()
+
+#### Analysis plots from the simulation component ####
 load('NEON_manuscript/Results/sim_results.rda')
 run_parameters <- expand.grid(n_series = c(4, 12),
                               T = c(72),
@@ -81,7 +182,7 @@ ggplot(drps_plot_dat %>%
   scale_fill_viridis(discrete = T, begin = 0.2, end = 1, guide = FALSE) +
   theme_bw() + coord_flip() + labs(x = '', y = 'Normalised DRPS calibration (lower is better)',
                                    title = 'Strong trend') -> plot2
-pdf('NEON_manuscript/Figures/Fig1_simulation_drps_missing_plot.pdf')
+pdf('NEON_manuscript/Figures/Fig2_simulation_drps_missing_plot.pdf')
 cowplot::plot_grid(plot1, plot2, ncol = 1)
 dev.off()
 
@@ -154,6 +255,6 @@ ggplot(coverage_plot_dat %>%
   ylim(0, 1) +
   theme_bw() + coord_flip() + labs(x = '', y = '90% interval coverage',
                                    title = 'Strong trend') -> plot2
-pdf('NEON_manuscript/Figures/Fig2_simulation_coverage_nseries_plot.pdf')
+pdf('NEON_manuscript/Figures/Fig3_simulation_coverage_nseries_plot.pdf')
 cowplot::plot_grid(plot1, plot2, ncol = 1)
 dev.off()