R Tutorial Part 2

Adapted from the state space modeling activity from Michael Dietz’s excellent Ecological Forecasting book

JAGS needs to be installed: https://sourceforge.net/projects/mcmc-jags/files/ rjags R package needs to be installed

Video Tutorial

Text Tutorial

Missing data

  • As we say at the end of the last video, adding NA’s to the end of the time-series produces forecasts from Bayesian state space models

  • This happens because the model attempts to determine the most likely value for NA’s based on the model and the data

  • This means that this type of modeling approach also works naturally with missing data

  • This is different from many types of time-series models that require continuously sampled time-series

  • To illustrate this we’ll add NA’s to a portion of our time series as an example of missing data

  • Go back to the line where we replaced the end of the time-series with NA’s and add some NA’s inside the training portion of the data

  • We’ll use two examples of missing data

  • For the first we’ll just a have data missing from a few random weeks

data$y[c(26, 50, 90, 260, 261, 262)] = NA
  • For the second we’ll mimic an entire missing year of data
data$y[year(time) == 2008] = NA
  • Now rerun the modeling and prediction steps
  • We can see that the model still fits and makes predictions like before
  • Can also see that estimates what the missing values are, along with their uncertainty
  • This is called “imputation”
  • For the single missing data points we can see that the predictions are consistent with the data and the model is somewhat confident in its imputations
  • But when an entire year of data is missing it we see similar behavior to our forecast
  • The model predicts basically no change and uncertainty becomes high
  • It peaks in the middle of the year because that is where there is the last empirical information to constrain the predictions
  • As it gets close to the end of the year the data from the next year starts to constrain it
  • And so we can see that this imputation is capturing the models uncertainty within the training data when no observed values are available
  • It’s actually quite high and this is why our forecasts from the model are highly uncertain, because there are no observations to constrain the model

What went wrong

  • So our first attempt at a forecast didn’t look overly promising
  • What do we do next
  • First let’s explore a little more about where the model when wrong
  • We’ve already plotted our observed and predicted time-series
  • So evaluate the model in another way we’ve learned and look at how the error in the point forecast changes with the forecast horizon
  • First we need to identify the data that is only for the forecast
  • That’s the end of the time series
forecast_data = (length(y)-51):length(y)
  • The we can do the analysis we did before but only with this piece of the time series
     sqrt((predictions[forecast_data] - y[forecast_data])^2))
  • Random walk worked well for a few time steps, but not for longer range forecasts
  • What’s interesting is that normally we’d expect the error to keep going up, but it actually comes back down
  • That’s because the increase in error represents the flu season
  • Once that season is over flu levels return to baseline and since that’s what the model predicted the error declines to near zero again
  • This is kind of pattern can be indicative of a model that is missing cyclic trend
  • The kind of thing we would model with either a seasonal effect or a cyclic predictor variable
  • In the past we’ve modeled this kind of cyclic pattern using seasonal lags
  • There are two other ways of modeling them
  • We can include the date as a predictor
  • Or we can identify the drivers generating the seasonal variation and include them as predictors

Dynamic linear modeling

  • Both of these approaches require adding covariates to our models
  • We describe models that include both a time-series component and predictors as “Dynamic linear models”

Data setup

  • Let’s start by putting together our predictor data
  • First we’ll download some weather data
# install.packages('daymetr')
weather = download_daymet(site = "Orlando",
                               lat = 28.54,
                               lon = -81.34,
                               start = 2003,
                               end = 2016,
                               internal = TRUE)
weather_data = weather$data
  • For date information this data has year and yday
  • yday is the Julian day (starts at 1 on Jan 1 and goes to 365)
  • Make this into a date column so we can combine with the weekly date data in Google Flu
weather_data$date = as.Date(paste(weather_data$year, weather_data$yday, sep = "-"),"%Y-%j")
weather_data = subset(weather_data, weather_data$date %in% time)
  • Let’s add two predictors from weather_data to our data object used by our model
data$Tmin = weather_data$tmin..deg.c.
data$yday = weather_data$yday
  • Finally we’re also going to log transform our response variable to help the JAGS model fit more effectively
data$logy = log(data$y)


  • Now we need to add these to our our random walk model in JAGS
  • But this gets even more complicated once there are predictors
  • So for today we’re going to take a short cut by using a small piece of software to help us do this
  • The ecoforecastR package will generate and run dynamic linear models using JAGS
  • We’ll start by adding minimum temperature to the model
  • We can fit the model using fit_dlm
  • We describe our model using a named list that includes obs for our observation and fixed for our fixed effects in the model
dlm = fit_dlm(model = list(obs="logy", fixed="~ 1 + X + Tmin"), data)
  • We can see that JAGS model that was generated
  • We can look at the parameters to see if things are converging
params = dlm$params
params <- window(dlm$params,start=1000) ## remove burn-in
plot(params, ask = TRUE)
  • We can visualize forecasts in the same way as before
out <- as.matrix(dlm$predict)
ci <- apply(exp(out),2,quantile,c(0.025,0.5,0.975))
plot(time, y)
lines(time, ci[2,])
lines(time, ci[1,], lty = "dashed")
lines(time, ci[3,], lty = "dashed")
  • This is better, but it’s still not great
  • We haven’t really captured the seasonality with temperature
  • So we could try adding in the julian day directly
dlm = fit_dlm(model = list(obs="logy", fixed="~ 1 + X + Tmin + yday"), data)
  • Now rerun the visualization code
  • This looks a lot better!
  • So we’ve got a pretty nice forecast here for the flu
  • But, it’s a little more complicated than this in most cases because Julian day is circular
  • In other words December 31st and January 1st are really close to one another
  • But a linear model on Julian Day thinks of them as far apart
  • So, this happens to work in our case but generally it won’t
  • We can fix this by doing a harmonic transformation of the Julian Day
  • But that’s a subject for another lesson so we’ll just stop here as an example