I’ve been asked this a couple of times, so perhaps it is worth a blog post.

For a simple training/test split, the Mean Absolute Scaled Error (MASE) (Hyndman & Koehler, 2006) is defined as and the Root Mean Squared Scaled Error (RMSSE) is defined as In both cases, for non-seasonal data, where is the seasonal period for seasonal data, and the sum in the numerator is over the test set, while the sum in the denominator is over the training set. The notation means the forecast of given data .

These measures are discussed in my forecasting textbook with George Athanasopoulos. The denominator is a scaling factor, introduced so that you can compare MASE or RMSSE values across series of different units. For example, are the forecasts of widget sales more accurate than the forecasts of electricity demand? If all your forecasts are in the same units, then you don’t need to remove the scale, and it is simpler to just use MAE or RMSE (i.e., only the numerators of the above equations).

Now, the question is, when we are doing time series cross-validation, and computing forecast accuracy over a series of training/test sets with a rolling origin, does it make sense to compute a different scaling factor each time? For example, a cross-validated MASE for -step forecasts is where the first observations form the smallest training set, and subsequent training sets increase one observation at a time.

An alternative approach would be to compute the scaling factor across all available data, rather than calculate it separately for each training set. Then MASE would become

Let’s think about the advantages of this alternative:

1. It is (slightly) faster. But the denominators are very fast to compute, so this really doesn’t make much difference.
2. It removes a source of variation from the calculation. The scaled errors will be more variable when the denominator changes with each test set, and that makes it harder to see the difference between forecasting methods.
3. It uses more data in computing the scaling factor, which reduces the uncertainty in the estimate of the accuracy measure. This is potentially important if , so that some training sets are relatively small compared to the available data.

As for disadvantages:

1. The measures are no longer true measures of forecast accuracy because the calculation potentially involves future observations. On the face of it, this seems important, but in reality it isn’t. The future observations aren’t affecting the forecasts, only the scaling factor, so there is no leakage involved.
2. One of my correspondents suggested that it changed the interpretation. I think the interpretability of these scaled measures is over-rated. They can only be interpreted as a ratio of out-of-sample accuracy to in-sample accuracy, and the two are not necessarily even over the same horizons. How the scaling factor is calculated doesn’t really affect the interpretation, because there is not much value in interpretation either way.

So I suggest that computing the scaling factor across all available training data when calculating cross-validated MASE and RMSSE values is a good idea.

That raises the question as to why shouldn’t we use all available data in the denominator when doing a simple training/test split? In that context, the only advantage that is relevant is #3 above, and unless the test set is particularly large, it shouldn’t make much difference. Nevertheless, I can’t see any real problem in using all available data in computing the scaling factor.

However, version 9.0 of the forecast package has now been released on CRAN. Thanks to Maximilian Muecke for helping with this release, and for modernising some of the old code.

Here are the main new features and changes.

USAccDeaths |> ets() |> forecast()
USAccDeaths |> auto.arima() |> forecast()

USAccDeaths |> snaive()

USAccDeaths |> rw_model(lag = 12) |> forecast()

library(purrr)
model_fns <- list(ets, auto.arima, rw_model, rw_model, mean_model)
args <- vector("list", length(model_fns))
args[[3]] <- list(lag = 12)
args[[4]] <- list(lag = 12, drift = TRUE)
models <- map2(model_fns, args, ~ do.call(.x, c(list(y=USAccDeaths), .y)))
names(models) <- map(models, function(u) u$method)
forecasts <- map(models, forecast)

The WAPE was introduced by Kolassa & Schütz (2007) who called it the MAD/Mean ratio. It is defined as where

• is the actual value at time ,
• is the forecast value at time

You can think of it as a weighted percentage error by writing it as where the weights are given by

It can also be considered a relative MAE where the comparison method has all forecasts equal to zero. In fact, would give the skill score relative to the forecast of zero for all .

Or you could think of it as like a MASE (Hyndman & Koehler, 2006) but with scaling based on the sum of absolute values on the test set rather than the sum of absolute differences on the training set.

This has some obvious advantages over the MAPE:

1. The MAPE is undefined when any actual value in the test set is zero. The WAPE is defined even when some actuals are zero. It is only undefined when all the actuals used in the denominator are zero.
2. Optimising the MAPE does not lead to a sensible point forecasts (Gneiting, 2011), but optimising the WAPE will lead to the median forecast.

However, I think there are a couple of problems that do not seem to have been widely recognized.

1. The resulting estimate is only consistent when the time series is stationary. So it should not be used with data that has trends, or seasonality, or heteroscedasticity.
2. It is quite possible to have all actuals in the test set equal to zero, especially with intermittent demand time series and small test sets. Then the denominator is zero, and the WAPE is undefined.

For these reasons, I think the Mean Absolute Scaled Error (MASE) (Hyndman & Koehler, 2006) is a better choice than the WAPE. The MASE is defined as where for non-seasonal data, and is the seasonal period for seasonal data, and the sum in the numerator is over the test set (of size ), and the sum in the denominator is over the training set (of size ). Because the denominator is defined on the training data, not the test data, it avoids the above problems with the WAPE:

• It is a consistent estimator provided the series is difference stationary, which is a much weaker condition than stationarity.
• The training data is also usually much longer than the test data, so it is much less likely to contain only zeros. In fact, if the training data did contain only zeros, then the obvious forecasts would all be zeros too.

I don’t want to suggest there are no problems with the MASE. Kolassa & Schütz (2007) point out one potential drawback of MASE — when there are structural breaks or outliers in the training data. Also, like the WAPE, optimising the MASE will lead to the median forecast, which is probably not what you want, especially when you have intermittent demand time series.

These days, if I want a scale-free accuracy measure, I prefer the Root Mean Squared Scaled Error (RMSSE): where again the sum in the numerator is over the test set, and the sum in the denominator is over the training set. This has all the advantages (and most disadvantages) of the MASE, but optimising it leads to the mean forecast rather than the median. It also aligns better with how models are estimated. Almost all models are estimated by minimising the sum of squared errors, so it makes sense to evaluate them using squared errors as well.

Point forecast reconciliation also works on means rather than medians, and is optimised using least squares, so it is more natural to evaluate using a squared error measure.

The fpp3 package is a companion to the book Forecasting: Principles and Practice (3rd edition, Hyndman & Athanasopoulos, OTexts). When you load the package, it loads the data and functions needed for the examples in the book.

library(fpp3)

── Attaching packages ──────────────────────────────────────────── fpp3 1.0.0 ──

✔ tibble      3.2.1     ✔ tsibble     1.1.4
✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
✔ tidyr       1.3.1     ✔ feasts      0.3.2
✔ lubridate   1.9.3     ✔ fable       0.3.4
✔ ggplot2     3.5.1     ✔ fabletools  0.4.2

── 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()

This provides a shorthand way of attaching the packages needed for most time series forecasting tasks, using packages from the tidyverts, just as the tidyverse package does for data manipulation and visualization.

In the most recent update, we have added many new data sets that can be used as examples by instructors, or for users to practice their forecasting skills.

Each of these new data sets were previously used in an exam for the forecasting subjects taught by George Athanasopoulos and me at Monash University. Three of our exams are available online for those who are interested to see how we assess forecasting students.

• 2021
• 2022
• 2023

Each year, we take a new data set, and write the exam around how to analyse the time series. We strongly believe in assessing students on their forecasting skills with real data, rather than using artificial data, or asking technical questions that are not relevant to real-world forecasting. As a result, our exams are quite different from anything we have seen elsewhere.

We also think it is best to use “fresh and local”¹ data sets. So the data are almost always from Australia, and included data up to about a month before the exam was written.

In the latest update to the fpp3 package, we have expanded these data sets to cover more series than we did in the exams, and in most cases updated them to include more recent observations.

The new data sets are:

• aus_births: Australian births data
• aus_fertility: Australian fertility rates
• aus_inbound: Monthly short term (<1 year) visitor arrivals to Australia
• aus_migration: Australian migration data
• aus_mortality: Australian mortality data
• aus_outbound: Monthly short term (<1 year) resident departures in Australia
• aus_tobacco: Australian cigarette and tobacco expenditure
• aus_vehicle_sales: Australian vehicle sales
• melb_walkers: Average daily total pedestrian count in Melbourne
• nsw_offences: Monthly offences in NSW
• ny_childcare: New York childcare data
• otexts_views: OTexts page views

Thanks to Nuwani Palihawadana and Shanika Wickramasuriya for most of the work on this latest update, completed at the NUMBAT hackathon last month.

I get this question a lot, so I thought it might help to explain some issues with AIC calculation.

First, the equation for the AIC is given by where is the likelihood of the model and is the number of parameters that are estimated (including the error variance). For a linear regression model with iid errors, fitted to observations, the log-likelihood can be written as where is the residual for the th observation. The AIC is then Since we don’t know , we estimate it using the mean squared error (the maximum likelihood estimator), giving where is a constant that depends only on the sample size and not on the model. This constant is often ignored. Thus, different software implementations can lead to different AIC values for the same model, since they may include or exclude the constant .

Now, let’s look at what R returns in a simple case using the lm() function.

set.seed(2023)
library(fpp3)
df <- tibble(
  time = seq(100),
  x = rnorm(100),
  y = x + rnorm(100)
)
fit <- lm(y ~ x, data = df)
AIC(fit)

[1] 275.6267

We can check how this is calculated by computing it ourselves.

mse <- mean(residuals(fit)^2)
n <- length(residuals(fit))
k <- length(fit$coefficients) + 1
# With constant
2*k + n*log(mse) + n*log(2*pi) + n

[1] 275.6267

# Without constant
2*k + n*log(mse)

[1] -8.161047

Clearly, AIC() applied to the output from lm() is using the version with the constant.

Now compare that with what we obtain using the TSLM() function from the fable package.

df |>
  as_tsibble(index = time) |>
  model(TSLM(y ~ x)) |>
  glance() |>
  pull(AIC)

[1] -8.161047

This is the AIC without the constant.

The situation is even more confusing with ARIMA models, and some other model classes, because some functions use approximations to the likelihood, rather than the exact likelihood.

Thus, AIC values can be compared across models fitted using the same functions, but not necessarily when models have been fitted using different functions.

I received this email today:

My team recently used some techniques found in your writings to perform forecasts … Our work has been well received by reviewers, but one commenter asked two questions that I was hoping you may be able to provide insight on.

First, they wanted to know if we could provide P-values for our prediction intervals. In our work, we said, “Observed rates were deemed significantly different from expected rates when they did not fall within the 95% PI.” This same language has been used by others published in the same journal. I am curious to hear your thoughts on giving P-values for these PIs and what the appropriate method for doing so would be (if any).

Second, they asked about making a correction for multiple comparisons. … I believe we could apply a Bonferroni correction to the PIs, but that feels too liberal. Moreover, I am curious if this is even called for given our statement of what is deemed significant and the fact that our prediction interval construction relies on a non-parametric method.

Here is my reply:

1. I don’t think this makes any sense. A p-value is the probability of obtaining observations at least as extreme as those observed given a null hypothesis. What’s the hypothesis here? In forecasting, we don’t usually have a hypothesis. Instead, we fit a model to the data, and make predictions based on the model. I guess you could make the null hypothesis “The future observations come from the forecast distributions”, and then the p-value for each future time period would be the probability of the tails beyond the observations. But it is well-known that the estimated prediction intervals are almost always too narrow due to them not taking into account all sources of variance. So the size of this test would not be well-calibrated. I think you’re better off pushing back rather than trying to meet the request.

2. A Bonferroni correction assumes independence between the intervals, and that is not true for PIs from a forecasting model. The future forecast errors are all correlated (with the strength of the correlation depending on the model and the DGP). Usually we just say that these are pointwise PI, and so we expect 5% of observations to fall outside the 95% prediction intervals. It is possible to generate uniform PI, which contain 95% of all future sample paths, but this is a little tricky due to the correlations between horizons. It could be done via simulation – simulate a 1000 future sample paths and compute the envelope that contains 950 of them.

It sounds like the reviewers are only familiar with inferential statistics, and not with predictive modelling. You could point them to Shmueli’s excellent 2010 paper “To explain or predict”, highlighting the differences between the two paradigms.

There is a widespread myth that NASA did not discover the hole in the ozone layer above the Antarctic because they had been throwing away anomalous data that would have revealed it. This is not true, but the real story is also instructive (Pukelsheim 1990; Christie 2001, 2004).

NASA had been collecting satellite data on Antarctic ozone levels using a Total Ozone Mapping Spectrometer (TOMS) since 1979, while British scientists had collected ozone data using ground sensors at the Halley Research Station, on the edge of the Brunt Ice Shelf in Antarctica, since 1957. Figure 1 shows average daily values from the NASA measurements in blue, and from the British observations in orange. There is a clear downward trend in the British data, especially from the late 1970s, which is confirmed with the NASA data. So why wasn’t the “ozone hole” discovered until 1985?

The British scientists had noticed the low ozone values as early as 1981, but it took a few years for the scientists to be convinced that the low values were real and not due to instrument problems, and then there were the usual publication delays. Eventually, the results were published in Farman et al. (1985).

Meanwhile, NASA was flagging observations as anomalous when they were below 180 DU (shown as a horizontal line in Figure 1). As is clear from the figure, this is much lower than any of the plotted points before the early 1980s. However, the 180 threshold was used for the daily measurements, which are much more variable than the monthly averages that are plotted. Occasionally daily observations did fall below 180, and so it was a reasonable threshold for the purpose of identifying instrument problems.

In fact, NASA had checked the unusually low TOMS values obtained before 1985 by comparing them against other available data. But the other data available to them showed ozone values of about 300 DU, so it was assumed that the satellite sensor was malfunctioning. The British Halley data were not available to them, and only after the publication of Farman et al. (1985) did the NASA scientists realise that the TOMS results were accurate.

In 1986, NASA scientists were able to confirm the British finding, also demonstrating that the ozone hole was widespread across the Antarctic (Stolarski et al. 1986).

This example reveals some lessons about anomaly detection:

• The NASA threshold of 180 was based on daily data, and was designed to identify instrument problems, not genuine systematic changes in ozone levels. The implicit assumption was that ozone levels varied seasonally, but that otherwise the distribution of observations was stable. All anomaly detection involves some implicit assumptions like this, and it is well to be aware of them.
• Sometimes what we think are anomalies are not really anomalies, but the result of incorrect assumptions.
• Often smoothing or averaging data will help to reveal issues that are not so obvious from the original data. This reduces the variation in the data, and allows more systematic variation to be uncovered.
• Always plot the data. In this case, a graph such as Figure 1 would have revealed the problem in the late 1970s, but it seems no-one was producing plots like this.

References

Christie, M. 2001. The Ozone Layer: A Philosophy of Science Perspective. Cambridge University Press.

Christie, M. 2004. “Data Collection and the Ozone Hole: Too Much of a Good Thing?” History of Meteorology 1: 99–105.

Farman, J C, B G Gardiner, and J D Shanklin. 1985. “Large Losses of Total Ozone in Antarctica Reveal Seasonal ClO/NO Interaction.” Nature 315 (6016): 207–10.

Pukelsheim, F. 1990. “Robustness of Statistical Gossip and the Antarctic Ozone Hole.” The IMS Bulletin 19 (4): 540–45.

Stolarski, R S, A J Krueger, M R Schoeberl, R D McPeters, P A Newman, and J C Alpert. 1986. “Nimbus 7 Satellite Measurements of the Springtime Antarctic Ozone Decrease.” Nature 322 (6082): 808–11.

This question arises most time I teach a forecasting workshop, and it was raised again in the following email I received today:

I have a time series that I have split into training and test datasets with an 80%-20% ratio. I fit a series of different models (ETS, BATS, ARIMA, NN etc) to the training data and generate my forecasts from each model. When evaluating the forecasts against the test set I find the model that gives the best outcome is an ARIMA(1,1,1) that was selected using the auto.arima function. My question is this, should I proceed to fit an ARIMA(1,1,1) to the whole data set, or should I use the auto.arima function again which may give me a slightly different (p,d,q) order as there is now an extra 20% of unseen data available to the forecast model? Any guidance would be greatly received.

If you only have one class of model to consider (e.g., only ETS or only ARIMA), then it is easy enough to select the model on all available data using the AIC, and use the selected model to forecast the future. But if you are selecting between model classes, then you need to use either a training/test split, or (preferably) a time-series cross-validation procedure.

If you use time-series cross-validation, then there would usually be different models selected for each training set, and the cross-validated error is a measure of how well the model class works for your data. In that case, there is no single model for the training data, and you are selecting the model class rather than a specific model. This makes it clear that you should then apply the selected model class to all the data, when forecasting beyond the end of the available data. In other words, if you choose ARIMA over ETS, then you would then fit an ARIMA model to all the data, and use that model to forecast the future.

You can think of a simple training/test split as a special case of time-series cross-validation, where there is a single fold. So the same argument applies. That is, you are selecting the model class that works best for your data, and so you should apply that model class to all the data, when forecasting beyond the end of the available data.

This example also illustrates why it is important to use a time-series cross-validation procedure, rather than a simple training/test split. In this case, the ARIMA model was selected because it happened to work best for the particular training/test split that was used. But if a different split had been used, then a different model might have been selected. So the model selection is not stable. By averaging over multiple folds using a time-series cross-validation procedure, you can get a more stable estimate of the model class that works best for your data.

We have taught from the book many times, but this year we decided to pre-record short videos for each section. Our students often prefer a video explanation than reading the textbook, and we thought other readers might appreciate hearing from us as well.

These videos are embedded in most sections of the book. So far, we’ve covered the sections that we include in our own courses, but we hope to eventually have videos for all sections. Most of these were done in a single take, so they are sometimes a little rough, but hopefully still useful.

You can view the entire playlist on YouTube.

Places are limited, so please book in early.

• New York, 11-12 July 2023. Followed by the New York R conference. (Register for the workshop as part of the NYR conference registration.)
• Chicago, 17-18 September 2023. Followed by the Posit conference. (Register for the workshop as part of the Posit conference registration.)

The workshop introduces the tidyverts set of packages. Further details about the workshop are here.

Australia has a problem with government data. Actually it has three problems with government data: 1. It is often kept secret. 2. If it is available, it is often out-of-date. 3. If it is available and timely, it is often in a form that makes any analysis difficult.

I think it would be better for the country if government data was available freely, immediately, and in a form that is useful for analysis. Of course, we should make an exception if there are privacy issues, or some other harm that would be caused by releasing it. Let me explain using some examples.

Take mortality data. During the pandemic, it has been important to know how many people died of any cause, so we could know the effect of the pandemic overall. Obviously some people were dying of COVID-19, but others might have been dying because they were unable to get treated when medical staff were overwhelmed by COVID patients. On the other hand, lockdowns may have reduced deaths due to road crashes, but perhaps they also affected deaths due to suicide. If we could compare the total deaths each week during the pandemic, with the corresponding totals in previous years, we could determine the overall effect of the pandemic on Australian mortality.

You would think, that during a global pandemic, having good mortality data would be important. But in June 2020, nearly six months after the start of COVID-19, the most recent available mortality data in Australia was from 2018. Eighteen months out of date! Think about that. For the first six months of the biggest public health event in 100 years, we had no official data on the effect of COVID-19 on Australian mortality. Eventually the Australian Bureau of Statistics got their act together and started producing provisional mortality data more frequently, but only after several of us complained loudly and publicly. Even now, the provisional mortality data available from the ABS is more than 3 months out of date. Contrast that to other countries. I could find mortality data on 38 countries, and Australia was the 5th worst for producing timely mortality data.

Another example concerns COVID-19 case numbers. There is still no reliable Australian government repository of daily COVID-19 cases by state. Some states are now producing historical data, but for most of 2020, when we really needed reliable information, the public information was incomplete. For much of the first two years of the pandemic, the state health departments were putting out their little dashboard images containing the numbers, but these were preliminary numbers, and did not include cases that were registered late, and some other data revisions. To do any serious analysis, you needed daily case numbers from the beginning of the pandemic, but these were not available on government websites until relatively recently. Some media organizations, and some individuals, were collating the case numbers from the dashboard images and putting them online in the form of spreadsheets, and people were using them to do analysis, but these data were usually inaccurate and subject to revisions. The state health departments generally didn’t update the initial numbers that were released, even though they had more reliable information. So the public data was inaccurate, and most people wanting to do any data analysis were relying on media outlets, or a few 14 year old boys running covidlive.com.au, to get even that.

For nearly three years, I have been part of the forecasting team appointed to provide advice to all of the Chief Health Officers of the states and territories of Australia. Every week, we produce forecasts of COVID daily case numbers for all states and territories. For that purpose, we were able to put together a relatively good data set of case numbers for all states, but we were explicitly forbidden to make the data publicly available, even though our data was more accurate than what was appearing in the media.

Similarly, our forecasts were kept secret even though they were being used to make policy decisions. Premiers would justify their policies by vaguely referring to “the modelling”, or occasionally “the Doherty modelling” (even though most of us are not at the Doherty institute), but we would have preferred to have our forecasts available. So the good data and the forecasts are kept secret, and what is available is of poorer quality, or out-of-date.

Why? There are no privacy issues here. No harm would be done by working more transparently. On the contrary, if everyone had access to the best available data, then the independent modelling that was being done would have been of a higher quality.

We use a forecasting ensemble, where we have several forecasting models, and we combine them to produce the final forecasts that are submitted to the various state governments each week. Because we can’t share the data, the only forecasts that are included are those from members of our team. Generally in forecasting, it is better to use a wide range of models, not rely on a select few. But we can’t do that in Australia because of government obsession with secrecy.

Compare that to the United States where there was an official repository of data set up early in the pandemic, and anyone could download it and produce forecasts, and submit those forecasts to the Centre for Disease Control for inclusion in their analysis. Therefore, the US forecasting ensemble that was being used for policy decisions was based on a much larger range of models, and anyone could contribute to it. The resulting forecasts are then published publicly, so anyone can see what is being forecast, and what information a government has available when making policy decisions.

I’ve focused on COVID, but similar problems arise in many other areas in Australia. We have a culture of secrecy around data that is damaging to our public discourse, it leads to worse analysis, it means less transparency in government, and it feeds distrust of government because it is not clear why decisions are being made. Making more data publicly available leads to a better society.

We assume that the residuals from the method are uncorrelated and homoscedastic, with mean 0 and variance . Let denote the time series observations, and let be the estimated forecast mean (or point forecast). Then we can write where is a white noise process. Let be the estimated -step forecast variance.

Recently I spent a few weeks visiting Professor Tommaso Di Fonzo at the University of Padova (Italy), and one of the things we discussed was finding a notation we were both happy with so we could be more consistent in our future papers.

This is what we came up with. Hopefully others will agree and use it too!

For readers new to forecast reconciliation, Chapter 11 of FPP3 provides an introduction.

We observe time series at time , written as . The base forecasts of given data are denoted by .

The focus this year is on communicating with data. As with all WOMBAT events, the purpose is to bring together analysts from academia, industry and government to learn and discuss new open source tools for business analytics and data science.

I’ve been using Disqus for more than 13 years, largely because it was the only available solution at the time I added comments. To make Disqus interface a little cleaner, I disabled all the advertising and as much of the other noise as possible, but it still looked like something from mySpace (for those of you who remember the 20th century).

But now there are several alternatives, and I’ve opted for giscus which is very lightweight, is built on Github Discussions, and is open source with no tracking or advertising. The other system I considered was utterances which is also hosted on Github, but uses issues rather than discussions. Consequently, comments on utterances can’t be threaded (with replies to previous comments). Also, giscus appears to have a much more active development team behind it.

The first step was to set up giscus on my blog. With quarto, this simply requires adding a few lines to the _metadata.yml file in the relevant folder. Here is what it looks like for me:

comments:
  giscus:
    repo:  robjhyndman/robjhyndman.com
    repo-id: "R_kgDOH5G3Uw"
    category: "Announcements"
    category-id: "DIC_kwDOH5G3U84CRUp9"
    mapping: "pathname"
    reactions-enabled: true
    loading: lazy
    input-position: "bottom"
    theme: "light"

Then I needed to set up giscus on the Github repo that hosts the website (robjhyndman/robjhyndman.com). The instructions on the giscus website make it very simple.

The last step was the hardest – how to migrate 4000 comments from Disqus to giscus. Here I followed the nice blog post of Maëlle Salmon to download the Disqus comments as an xml file, and wrangle them into a tibble. Then I needed to use the GraphQL API for Github Discussions to generate all the comments on the Github repo. Fortunately, Mitch O’Hara-Wild came to my rescue (as usual), and helped with some of this code. The resulting code is here if anyone wants to try to do the same. You will need to change some specific details in lines 9-13. Everything else should work as it is.