Predictive Analytics in Tableau Part 9: Automatic Time Series Analysis

Today, we will talk about automatically performing time series analysis using Tableau 8.1's new R functionality.  Up to now, we've looked at ways to determine which type of exponential or ARIMA model fits the data best.  However, the number of models can be overwhelming to an analyst with less statistical exposure.  Today, we'll look at how to let R choose the model for you.  As usual, we will use the Superstore Sales sample data set in Tableau.

Let's imagine that we came up with a candidate list of models and want to see which ones are better for our type of data.

Too Many models

Here, we can see that charts can easily get very messy to see which models are best when there are a lot of models.  So, what better way is there?  Now, we're going to introduce two functions that will help.  The function ets() will create an appropriate Exponential model for your data and auto.arima() will create an appropriate ARIMA model for your data.  Like always, the code is identical until the last couple of lines.  You can find the full code in the appendix at the end of this post.  Here are the last couple of lines:

Exponential Model

fit <- ets(timeser)

c(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])

ARIMA Model

fit <- auto.arima(timeser)
(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])

Now, let's look at our chart.

Automatic Models

As you can see, these calculations quickly gave us two "good" candidate models to look at.  You might notice that neither of these models fits the data extremely well.  However, we could look at these models and see that the ARIMA model doesn't seem to be able to account for any of the variability, which is a very bad thing.  The Exponential model seems much better, but still needs some tweaking.  Therefore, we could throw out all of the ARIMA models and focus on creating a better Exponential model.  In just a couple of minutes, we were able to narrow our modeling down to a single family!  That's a huge time saver.  Perhaps we could even improve this process by adding new functions and new techniques.  But, that's a topic for another post.  Thanks for reading.  We hope you found this informative.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

Appendix

ARIMA (Auto)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- auto.arima(timeser)
    c(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])
",


ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

Exponential (Auto)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- ets(timeser)
    c(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])
",


ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

Predictive Analytics in Tableau Part 8: ARIMA Time Series

Today, we will talk about creating ARIMA time series models using Tableau 8.1's new R functionality.  In layman's terms, an ARIMA model uses three different numeric parameters to make varying types of time series models.  The ARIMA family is one of the most researched and respected families in the field of Time Series Analysis.  It is also used by a number of automated prediction tools like SQL Server Analysis Services.  As usual, we will use the Superstore Sales sample data set from Tableau.

The most important part of creating an ARIMA model is choosing the parameters.  The parameters are as follows:

AR  - Autoregressive
I      - Integrated
MA - Moving Average

We won't delve heavily into what each of these parameters does.  However, we will talk about methods for selecting them.  First, let's take a look at our data.

Sales by Month

Now, let's look at choosing the AR parameter.  In order to do this, we need to look at the Partial Autocorrelation Function (PACF).  Just like with the previous posts, almost all of the code is used to create the time series.  Only the last line is used to give us the PACF values.  The full code can be found in the appendix at the end of the post.

rep(pacf(timeser, plot=FALSE)$acf,5)[1:len.orig]

Now, let's look at these values.

Sales by Month (PACF)

The question we want to ask is "Starting at the beginning, how many consecutive values are above .25 or below -.25?"  We can easily accomplish this with colors.

AR Colors

Sales by Month (PACF with Color)

We can see that there are no values above .25 or below -.25.  So, there doesn't seem to be an AR component to this model.  Next, let's move on to the MA component using a very similar technique.  Moving Average components can be analyzed using the Autocorrelation Function (ACF).  Here's the last line of code.

rep(acf(timeser, plot=FALSE)$acf,5)[1:len.orig]

Now, let's look at the values.

Sales by Month (ACF with Color)

The first value is always guaranteed to be 1 and should be ignored.  We see that the values after one aren't really worth noting.  Therefore, the MA component in the model is 0 and the AR component is 0.  Now, let's look at the Integrated portion of the model.  R has a nice built-in function called ndiffs() that tells us what the I component should be.  Here's the last line of code.

ndiffs(timeser)

Now, let's see the value.

Sales by Month (I Component)

Here, we see that our I component should be zero.  This leaves us with a bit of an interesting conundrum.  An ARIMA(0,0,0) model is basically useless.  It is guaranteed always return the mean (average) value.  So, this means that the ARIMA family is not a good fit for this data.  Just for kicks, let's swap this data out for another data set that does fall in the ARIMA family.  So, let's look at a data set with monthly Unemployment values for the US.

Unemployment by Month

Now, let's look at the ARIMA components.

Unemployment by Month (ARIMA Components)

We can see that there are 3 consecutive AR values above .25 or below -.25, and I of 0, and 2 MA values (the first of which we ignore).  So, our model should be ARIMA(3,0,1).  Now, let's see what our forecasts look like.

Unemployment by Month (Forecast)

We can see that the ARIMA model seems to be pretty good at predicting the up-and-down motion of the data.  However, the recent rise in unemployment has caused some variability that the model doesn't handle very well.  Fortunately, there are plenty of other families of models out there.  Thanks for reading.  We hope you found this informative.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

EDIT: An anonymous user noted an error in the original appendix.  The R variables were defined as 

year.orig <- .arg1
month.orig <- .arg2

while our assignment at the end of the segment was defining

ATTR( MONTH( [Order Date] ) ) = .arg1
ATTR( YEAR( [Order Date] ) ) = .arg2

This has now been fixed an all of the code should be valid.  Many thanks to our readers for pointing that out.

Appendix

AR Component

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg2
    month.orig <- .arg1
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    rep(pacf(timeser, plot=FALSE)$acf,5)[1:len.orig]
",

ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

I Component

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg2
    month.orig <- .arg1
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    ndiffs(timeser)
",

ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

MA Component

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg2
    month.orig <- .arg1
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    rep(acf(timeser, plot=FALSE)$acf,5)[1:len.orig]
",

ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

Forecast (ARIMA(0,0,1))

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg2
    month.orig <- .arg1
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- arima(timeser, order=c(0, 0, 0))
    c(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])
",

ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

Forecast (ARIMA(3,0,1))

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg2
    month.orig <- .arg1
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- arima(timeser, order=c(3, 0, 1))
    c(rep(NA, len.new), forecast(fit)[[4]][1:hold.orig])
",

ATTR( MONTH( [Order Date] ) ), ATTR( YEAR( [Order Date] ) ), SUM( [Sales] ), [Months to Forecast] )

Predictive Analytics in Tableau Part 7: Double and Triple Exponential Time Series

Today, we will talk about creating Double and Triple Exponential Time Series using Tableau 8.1's new R functionality.  If you read our previous post, Single Exponential Time Series, you remember that we were only able to predict a constant value for all future observations.  Trend and Seasonality were not considered using the Single Exponential.  These models should fix all of that.  Once again, we will use the Superstore Sales sample data set from Tableau.

At its core, a single exponential time series can predict values without trend (increase or decrease over longer periods of time) or seasonality (repetition of patterns at regular intervals, i.e. yearly, monthly, etc.).  The double exponential time series adds trend to the model while the triple exponential time series adds trend and seasonality.  We're also adding the model that accounts for seasonality without accounting for trend.  We're not sure what this model is called, so we'll call it the Reverse Double Exponential Model for now.  If you have any idea what it's called or if it's invalid for some reason, please let us know in the comments.  The code is almost identical like before, just with the second to last line changed.  You can find the full code in the appendix at the end of this post.

Double Exponential

fit <- HoltWinters(timeser, gamma=FALSE)

Reverse Double Exponential

fit <- HoltWinters(timeser, beta=FALSE)

Triple Exponential

fit <- HoltWinters(timeser)

Now, let's look at these predictions.

Sales by Month (Predictions)

We can see that the green line (Double Exponential) is straight.  This is exactly what we expected.  However, there is some definite seasonality in this data that we need to account for.  The orange line (Triple Exponential) seems to follow that seasonality, but the trend seems to have pushed it too high.  The red line (Reverse Double Exponential) appears to be the "Goldilocks" model.  It accounts for the seasonality, but doesn't let the trend drag it too high.  For further examination, let's look at the 95% confidence intervals for these models to see which models fit the data more tightly.

Sales by Month (Bounds)

We can easily see that the green bounds (Double Exponential) are so spread out that they become basically worthless.  Also, the orange bounds (Triple Exponential) and the red bounds (Reverse Double Exponential) seem to the same except that the orange bounds are shifted up slightly higher, which is exactly what we saw with the predictions as well.  Therefore, we have good evidence to say that the Reverse Double Exponential model is the best fit for this data.  Now, let's use the same trick as last week to predict new values by adding blank rows to our data set.

Sales by Month (Reverse Double Exponential with Bounds)

Now, we have seemingly accurate predictions using a well-respected time series model.  However, we're not done with this yet.  There are other families of time series models we could use.  There's also artificial neural networks, naive bayes, and much more.  Thanks for reading.  We hope you found this informative.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

Appendix

Double Exponential

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$mean)
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),
 [Months to Forecast] )

Double Exponential (Lower 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Double Exponential (Upper 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Reverse Double Exponential

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$mean)
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),
 [Months to Forecast] )

Reverse Double Exponential (Lower 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Reverse Double Exponential (Upper 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Triple Exponential

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$mean)
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Triple Exponential (Lower 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Triple Exponential (Upper 95%)

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Predictive Analytics in Tableau Part 6: Single Exponential Time Series

Before we start today's topic, we want to direct you towards a great source of information about Tableau 8.1's R functionality.  This page was compiled by members of the Tableau community and can be found here.  Now, on to the demonstration.

Today, we will talk about creating time series analyses using Tableau 8.1's new R functionality.  In our previous demonstrations, we used varying regression techniques to predict values.  These techniques took into account the relationships between the variables, but not across time.  Time Series Analyses approach from the opposite direction.  They attempt to predict new values by looking at previous (or future) values in time, but only considering one variable.  For this demonstration, we will use the Superstore Sales sample data set from Tableau.

As always, the first step is to look at the data we want to forecast.  Let's look at Total Sales per Month.

Sales by Month

We can see that there is quite a bit of variation in this data.  Now, we run into our first conundrum.  I don't know of an easy way to allow Tableau to see dates outside of what's in the data.  Therefore, we can only forecast for dates we already have in the data.  This is actually a good and bad thing for us, as we'll soon see.  If you have a system that includes budget data for the future, then the dates will already be available and you won't have this problem.  Alas, we do not have such luck.  So, let's first hold out the last six months of our data and create the model using only data up to June 2012.  That way, we have something to compare our forecasts to.  In order to do this, we create an integer parameter.

Months to Forecast

Now, we can dynamically set the number of months we want to forecast.  Next, we need to create the code to do this.  But, before we do that, we need to install the "forecast" package on the R Server.

Install Packages

If you don't have the access to do this, then you will probably need to talk to the administrator on your R server.  Fortunately, we host the server on our local machine; so, there is no issue.  Now, back to the forecasting.

This code is long and somewhat complex.  So, we'll only touch on the important parts.  The code for Sales (Forecast) is as follows:

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1],
                     dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$mean)
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),
 [Months to Forecast] )

First, we need to install the forecast library using library(forecast).  Then, we need to sort the data by date, and hold out the last few months.  How many months get held out depends on the value of our [Months to Forecast] parameter.  Finally, we fit our time series as a single exponential model and output the results.  Let's see what it looks like.

Sales by Month (with Forecast)

As you can see, the forecast is pretty close to the true values.  This type of model cannot pick up on the trend or seasonality though.  We'll deal with that in a later post.  For now, let's see how well this model fits by also plotting the 80% and 95% confidence intervals.  You should remember confidence intervals from our posts about regression.  The code for the intervals are identical to the code for the forecast, with the exception of the last line of code.  I'll post the final lines here.  If you want to see the full code, please refer to the appendix at the end of this post.

Sales (Forecast) (Lower 80%)

c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,1])

Sales (Forecast) (Lower 95%)

c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,2])

Sales (Forecast) (Upper 80%)

c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,1])

Sales (Forecast) (Upper 95%)

c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,2])

Now, let's see what the intervals look like

Sales by Month (with Forecast and Bounds)

We can see that all of the actual values fall well within the 95% bound.  This means that our model fits our data reasonably well.  Now, let's talk about how to use this model to predict future values.  The simplest way we could find to do this is to add some filler lines to your data set.  In our case, we're using an Excel spreadsheet.  So, we add the following lines to our data.

Adding Empty Dates

You should note that if you are doing forecasts for a specific type of product, customer, etc., you will need to duplicate these dates for every combination of dimensions you want to forecast over.  Now, when we refresh our worksheet, we have our forecasts.

Sales by Month (with Future Forecast and Bounds)

One very important thing to notice about these forecasts is that the bounds are extremely far from the graph.  This would seem to imply that the model is not a good fit for the data.  Therefore, we would want to use a different type of time series model for it.  Unfortunately, that's going to have to wait until the next post.  Thanks for reading.  We hope you found this informative.

P.S.

If anybody can think of a better way to get future forecasts for dates that don't exist in the data, please let us know in the comments.  Also, this method is not robust against missing data.  If you're data has some missing months/week/etc., you may have to alter the method slightly.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

APPENDIX:

The code for Sales (Forecast) (Lower 80%) is as follows:

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,1])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),
 [Months to Forecast] )

The code for Sales (Forecast) (Lower 95%) is as follows:

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$lower[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

The code for Sales (Forecast) (Upper 80%) is as follows:

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,1])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

The code for Sales (Forecast) (Upper 95%) is as follows:

SCRIPT_REAL("
    library(forecast)

    ## Creating vectors

    hold.orig <- .arg4
    len.orig <- length( hold.orig )
    len.new <- len.orig - hold.orig[1]

    year.orig <- .arg1
    month.orig <- .arg2
    sales.orig <- .arg3

    ## Sorting the Data

    date.orig <- year.orig + month.orig / 12
    dat.orig <- cbind(year.orig, month.orig, sales.orig)[sort(date.orig, index.return = TRUE)$ix,]
    dat.new <- dat.orig[1:len.new,]

    ## Fitting the Time Series

    timeser <- ts(dat.new[,3], start = c(dat.new[1,1], dat.new[1,2]), end = c(dat.new[len.new,1], dat.new[len.new,2]), frequency = 12)
    fit <- HoltWinters(timeser, beta=FALSE, gamma=FALSE)
    c(rep(NA,len.new), forecast(fit, hold.orig[1])$upper[,2])
"
, ATTR( YEAR( [Order Date] ) ), ATTR( MONTH( [Order Date] ) ), SUM( [Sales] ),

 [Months to Forecast] )

Predictive Analytics in Tableau Part 5: Polynomial, Exponential and Piecewise Regression

Today, we will talk about some more types of regression using Tableau 8.1's new R functionality.  We previously talked about prediction using Linear Regression.  But, what if the relationship isn't linear?  In the real world, linearity is rarely the truth.  For this demonstration, we will use the same data set we used in Parts 1, 2, and 3.

First, let's look at our data.

DJIA vs. Foreign by Year

As you can see, this data is not related linearly.  There is very little change in DJIA when Foreign is below 140.  However, when Foreign surpasses 140, DJIA skyrockets.  If we tried to use a line to predict these values, we would get something like this:

DJIA vs. Foreign by Year (with Linear Trend)

This line is pretty far off from most of those points.  Fortunately, Tableau offers some more built-in trend lines.  Here are some good ones:

DJIA vs. Foreign by Year (with Quadratic Trend)

DJIA vs. Foreign by Year (with Exponential Trend)

As you can see, both of these models fit the data pretty well.  Now, let's see how to do them using R.

DJIA (Quadratic)

DJIA (Exponential)

We see that all we have to do is make a new variable with whatever function we want and add it to the model.  As far as R knows, these are two completely different variables.  We could add all of the other variables to the mix as well if we wanted to.  We could easily model Foreign as a quadratic while modeling Consumer as a logarithm.  One important thing to note is that if you include a polynomial, you should include all lesser degrees.  Simply put, if you have a x^2, then you must also have x.  If you have x^3, you must also have x^2 and x.  We won't go into detail about why you should do this at this time.

Looking back at the earlier trends, we're not happy with either of these trends.  The quadratic trend has that troubling curvature at the left side while the exponential trend doesn't seem strong enough to capture the upward curvature.  Now, let's make a new model that Tableau doesn't even have!

What if we believe that there are actually two trends here?  Let's imagine that the relationship is linear for small values of Foreign (less than 135) and a different linear relationship for larger values of Foreign?  No problem!  We can combine our models!

SCRIPT_REAL("

    djia <- .arg1

    fore <- .arg2

    th <- 135

    fore2 <- fore^2

    smallfore <- fore[fore<th]

    smalldjia <- djia[fore<th]

    largefore <- fore[fore>=th]

    largedjia <- djia[fore>=th]

    smallfit <- lm( smalldjia ~ smallfore )

    largefit <- lm( largedjia ~ largefore )

    c(smallfit$fitted, largefit$fitted)

",

SUM( [DJIA] ), SUM( [FOREIGN] ) )

Finally, let's see the results.

DJIA vs. Foreign by Year (with Piecewise Prediction)

This model fits our data so much better!  This method is called Piecewise Regression.  The amazing thing about R is that there's a method for predicting anything you want.  You can even create your own prediction method if you need to.  Thanks for reading.  We hope you found this informative.

P.S.

We're curious how Tableau keeps track of which value in the output vector corresponds to each input.  This was a huge obstacle in Part 3 as well.  If you have ideas, let us know in the comments.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

Predictive Analytics in Tableau Part 4: Logistic Regression

Today, we're going to talk about performing Logistic Regression using Tableau 8.1's R functionality.  Logistic Regression is very similar to Linear Regression, which we saw in the previous posts in this series.  However, Logistic Regression is designed to predict binary (Yes/No, 1/0) outcomes.  A very simple example is "Will this customer buy our product if we advertise to them?"  For this exercise, we will use the ubiquitous AdventureWorks data set from Microsoft.  If you've ever seen a Microsoft Data Mining demo, you've seen this data set.  Let's start by looking at our data.

Customer Demographics

As you can see, we have quite a bit of information about our customers, as well as whether or not they purchased a bike from us in the past.  Now, let's look at a logistic model with only one predictor so that we can understand how it works.

Purchased Bike (Predicted by Age)

As you can see, this code is extremely similar to the code for creating a linear regression model.  The only differences are that this model uses a more complex function, glm(), and an extra parameter, family = binomial( logit ).  Now, let's see what the predictions look like.

Predictions by Age

Some of you will immediately ask why this function returns decimals when we asked it to predict a Yes/No response.  Strictly speaking, logistic regression does not predict a Yes/No response, it predicts the probability of a particular response.  In other words, it tells us how likely this person is to buy a bike.  It's up to us to decide how we want to use these probabilities.  Let's look at these predictions in another way.

Predictions by Age (Scatterplot)

As you can see, the probability of buying a bike decreases as the customer gets older.  This is an important, and not very surprising, discovery.  Now, how do we turn these probability into actual predictions?  That's up to us!  An easy way is to say "If the chance is greater than 50%, we say they will buy.  If it's less than 50%, we say they won't."  Let's see what this gets us.

Predictions by Age (Classified)

This procedure doesn't seem to be very accurate.  Perhaps it's because we're only giving it one predictor.  Let's throw the rest of our variables in there and see if it gets better.

Predictions (Classified)

These predictions are much better, but still not as accurate as we'd like.  Unfortunately, we couldn't find an easy way to look at these predictions in aggregate using this method.  So far, it seems that the new R implementation is really good for generating predictions.  However, it falls a little short when it comes to examining the model.  Fortunately, there's even more we could do here.  We could try a different model for this data, such as an artificial neural network or a bayesian model.  Maybe one of our readers can find a really neat way to display this data that sums it up nicely.  The rest is up to your imagination.  We hope that you found this informative.  Thanks for reading.

Brad Llewellyn
Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn

Predictive Analytics in Tableau Part 3: Validating the Accuracy of your Predictions

Today, we will talk about validating the accuracy of your predictions using Tableau 8.1's R functionality.  In the previous posts in this series, we showed you how to create statistical models to predict values for your business.  However, we ignored a very important issue.  A model is supposed to be pretty good at predicting values that you used when you created it.  However, the true test of a model is how well it can predict values it doesn't know the answer to.  For this examination, we will use the same data set and design as we used in the first post in the series, Predictive Analytics in Tableau Part 1: Simple Linear Regression.

When you are creating a model, you typically want to "hold out" a portion of your data to use for testing later.  For instance, if we hold out 25% of our data, then we only use the remaining 75% when we created our model.  The 75% is known as the "Training Set" because it is used to train the model.  The 25% is known as the "Testing Set" because it is used to test the model after it has been trained.  First, let's see how you would create these sets.

Consumer vs. Foreign

This is the data we have for creating our model.  As you can see, we only have 12 observations.  So, we would expect the training set to be 9 observations, and the testing set to be 3 observations.  Let's see how we would create this.

Set

The code is not too complex.  Basically, it creates a 0/1 values for each year, with a 25% chance of 1 and 75% chance of 0.  Then, it changes the labels on these values to be more readable.  Now, let's see what it looks like on our chart.

Consumer vs. Foreign (with Set)

As you can see, it put 8 observations in our training set, and 4 observations in the testing set.  This isn't exactly 25%, but it's only one observation off.  You should be aware that Tableau is going to query the R server every time you refresh this chart.  This means that the observations in the training set WILL change.  This is exactly the way we want it to work.  We want our training set to be as similar as possible to our testing set, with the exception of the size. Now, let's see how we could create our regression model using only our training set.  The code for "Foreign (Predicted)" is below:

SCRIPT_REAL( "

    cons <- .arg1

    fore <- .arg2

    set <- .arg3

    fore.train <- fore[set=='Training']

    cons.train <- cons[set=='Training']

    fit <- lm( fore.train ~ cons.train )

    dat <- data.frame(cbind(fore, cons))

    names(dat) <- c( 'fore.train', 'cons.train' )

    predict(fit, dat, interval = 'prediction')[,1]

"

, SUM( [CONSUMER] ), SUM( [FOREIGN] ), [Set] )

This code simply uses the training data to create the model, then predicts values for all of the data.  Let's see this on our scatterplot.

Consumer vs. Foreign (with Set, Prediction, and Trend)

This is just like the scatterplots we've been looking at, with a couple of additions.  First, we added our Regression model as the grey line.  Then, we added Tableau's Trend Line as the black line.  Our model was not trained using the testing data, but Tableau's trend line was.  So, what's the difference?  Not much according to this graph.  The lines are very close to one another.  This is one piece of evidence saying that this model predicts well.  Let's look at it differently.

Consumer vs. Foreign (with Set, Prediction Interval, and Trend)

Here, we see that not only do all of the training and testing values fall within the bands, but so does the trend line.  This is a very good sign.  For those of you that are less graphical and more numeric, let's see this on a table.

Consumer vs. Foreign (with Set, Prediction, and % Diff)

As you can see, our predictions can be off by as much as 35%.  You might ask "Why did we talk so highly about the model before when it predicts so poorly here?"  The answer is simple.  The model did the best it could with the amount of data it was given.  When you give a model ten values to look at, how accurate do you expect it to be?  The answer is not very.  However, in the business world, it's exceedingly rare to have extremely small data sets.  For instance, if you have monthly sales values for 5 years, that's 60 values.  You could easily hold out 10 or 15 of those and still get a good model.

There are no right answers when it comes to forecasting.  Everything is up to interpretation.  That's why so many companies put so much money into statisticians to develop accurate forecasts.  Even more to that point, perhaps linear regression isn't appropriate for this data set.  Maybe you would want to use a time series method, multiple regression, or a bayesian model.  Your predictive abilities are limited only by your data and your imagination.  We hope you found this informative.  Thanks for reading.

Brad Llewellyn
Associate Data Analytics Consultant
Mariner, LLC
llewellyn.wb@gmail.com
https://www.linkedin.com/in/bradllewellyn