forecast() of time series regression model does not work in R shiny

I am making an application that has a user input (selectInput) for multiple data frames . The user should obtain a dygraph of the data and its prediction interval (30%, 50%, 70%) for the time series regression model when a data frame input is selected. The forecast function and its dygraph for this model work fine in R but they aren't working in R shiny as I keep getting the Warning: Error in <-: object is not a matrix message. I suspect it might have something to do with the forecast() and the reactive() functions as I could not print out the output of the forecast function either (it also returns the same error message). The almost same code for the regression model is also used for the ARIMA model and the dygraph plus the forecast function for this model work! The differences in code between the two models are:

For Arima Model,

  1. fit <- auto.arima(train, stepwise = FALSE, approximation = FALSE)
  2. ARIMA.mean <- forecast(fit(), level = c(30,50,70), h = 36, biasadj = TRUE)

For Regression Model,

  1. fit2 <- tslm(train ~ trend + relevel(season, ref = which.min(tapply(train, cycle(train, FUN = sum)))))
  2. fit.30 <- forecast(fit.train, h = 36, level = 30, biasadj = TRUE)

fit.50 <- forecast(fit.train, h = 36, level = 50, biasadj = TRUE)

fit.70 <- forecast(fit.train, h = 36, level = 70, biasadj = TRUE)

Note: I have to split the forecast function into three different functions for each prediction level because for some reason levels = c(30,50,70) does not work for the regeression model.

The full code for the app is as follow:

ui <- fluidPage(
  headerPanel(h1("Time Series Analysis")),
  headerPanel(""),
  sidebarPanel(
    div(style = "font-size:85%;",uiOutput(outputId = "DATA")),
    width = 2),
  mainPanel(
    tabsetPanel(
      tabPanel("Tab 2", dygraphOutput(outputId = "ARIMA", width = "100%"))
      tabPanel("Tab 3",
                   dygraphOutput(outputId = "TSLM", width = "100%"))
      )))
)

server <- function(input, output){
 output$DATA <- renderUI({
    selectInput(inputId = "dataset",
                label = "please choose a data set", width = "200px",
                selected = NULL, multiple = FALSE, 
                choices = list(...))
  })
########################### TAB 2 #########################################
fit <- reactive({
    if(is.null(input$dataset)){return()}
    df <- get(input$dataset)
    tsdata <- ts(df$FixedCounts, frequency = 12,
                 start = c(min(df$year), min(df[df$year == min(df$year), "month"])),
                 end = c(max(df$year), max(df[df$year == max(df$year), "month"])))
    train <- window(tsdata,
                    start = c(start(time(tsdata))[1], match(month.abb[cycle(tsdata)][1], month.abb)),
                    end = c(floor(time(tsdata)[floor(length(tsdata)*0.8)]),
                            match(month.abb[cycle(tsdata)][floor(length(tsdata)*0.8)], month.abb)))
    fit <- auto.arima(train, stepwise = FALSE, approximation = FALSE)
  })
  output$ARIMA <- renderDygraph({
    ARIMA.mean <- forecast(fit(), level = c(30,50,70), h = 36, biasadj = TRUE)
    graph <- cbind(actuals = ARIMA.mean$x, pointfc_mean = ARIMA.mean$mean,
                   lower_70 = ARIMA.mean$lower[,"70%"], upper_70 = ARIMA.mean$upper[,"70%"],
                   lower_50 = ARIMA.mean$lower[,"50%"], upper_50 = ARIMA.mean$upper[,"50%"],
                   lower_30 = ARIMA.mean$lower[,"30%"], upper_30 = ARIMA.mean$upper[,"30%"])
    dygraph(graph, main = ARIMA.mean$method, ylab = "Monthly Visitors") %>%
      dySeries(name = "actuals") %>%
      dySeries(name = "pointfc_mean", label = "forecast") %>%
      dySeries(name = c("lower_30", "pointfc_mean", "upper_30"), label = "30% PI") %>%
      dySeries(name = c("lower_50", "pointfc_mean", "upper_50"), label = "50% PI") %>%
      dySeries(name = c("lower_70", "pointfc_mean", "upper_70"), label = "70% PI")
  })
########################### TAB 3 #########################################
  graph <- reactive({
# get the selected data
    if(is.null(input$dataset)){return()}
    df <- get(input$dataset)
# convert the data into a time series object
    tsdata <- ts(df$FixedCounts, frequency = 12,
                 start = c(min(df$year), min(df[df$year == min(df$year), "month"])),
                 end = c(max(df$year), max(df[df$year == max(df$year), "month"])))
# split the data into training (80%) and test (20%) sets
    train <- window(tsdata,
                    start = c(start(time(tsdata))[1], match(month.abb[cycle(tsdata)][1], month.abb)),
                    end = c(floor(time(tsdata)[floor(length(tsdata)*0.8)]),
                            match(month.abb[cycle(tsdata)][floor(length(tsdata)*0.8)], month.abb)))
# fit the model using regression (with seasonal dummy) method
    fitted <- tslm(train ~ trend + relevel(season, ref = which.min(tapply(train, cycle(train, FUN = sum)))))
    fit.30 <- forecast(fitted, h = 36, level = 30, biasadj = TRUE)
    fit.50 <- forecast(fitted, h = 36, level = 50, biasadj = TRUE)
    fit.70 <- forecast(fitted, h = 36, level = 70, biasadj = TRUE)
    graph2 <- cbind(actuals = fit.30$x, pointfc_mean = fit.30$mean,
                   lower_70 = fit.70$lower, upper_70 = fit.70$upper,
                   lower_50 = fit.50$lower, upper_50 = fit.50$upper,
                   lower_30 = fit.30$lower, upper_30 = fit.30$upper)
  })
  output$TSLM <- renderDygraph({
    graph <- graph()
      dygraph(graph) %>%
        dySeries(name = "actuals") %>%
        dySeries(name = "pointfc_mean", label = "forecast") %>%
        dySeries(name = c("lower_30", "pointfc_mean", "upper_30"), label = "30% PI") %>%
        dySeries(name = c("lower_50", "pointfc_mean", "upper_50"), label = "50% PI") %>%
        dySeries(name = c("lower_70", "pointfc_mean", "upper_70"), label = "70% PI")
  })
}

shinyApp(ui = ui, server = server)

Any help would be greatly appreciated.
Many thanks.

If I was going to debug this (after fixing your syntax errors and restyling), I would start by pulling out the big fit reactive into its own function:

fit <- reactive({
  if (is.null(input$dataset)) {
    return()
  }
  df <- get(input$dataset)
  tsdata <- ts(df$FixedCounts,
    frequency = 12,
    start = c(min(df$year), min(df[df$year == min(df$year), "month"])),
    end = c(max(df$year), max(df[df$year == max(df$year), "month"]))
  )
  train <- window(tsdata,
    start = c(start(time(tsdata))[1], match(month.abb[cycle(tsdata)][1], month.abb)),
    end = c(
      floor(time(tsdata)[floor(length(tsdata) * 0.8)]),
      match(month.abb[cycle(tsdata)][floor(length(tsdata) * 0.8)], month.abb)
    )
  )
  fit <- auto.arima(train, stepwise = FALSE, approximation = FALSE)
})

This appears to take a single data frame as an input:

build_arima <- function(df) {
  tsdata <- ts(df$FixedCounts,
    frequency = 12,
    start = c(min(df$year), min(df[df$year == min(df$year), "month"])),
    end = c(max(df$year), max(df[df$year == max(df$year), "month"]))
  )
  train <- window(tsdata,
    start = c(start(time(tsdata))[1], match(month.abb[cycle(tsdata)][1], month.abb)),
    end = c(
      floor(time(tsdata)[floor(length(tsdata) * 0.8)]),
      match(month.abb[cycle(tsdata)][floor(length(tsdata) * 0.8)], month.abb)
    )
  )
  auto.arima(train, stepwise = FALSE, approximation = FALSE)
})

Now you can interactively experiment with this function to figure out what's going wrong. If it turns out that build_arima() is fine, repeat the same operation for the graph reactive.

1 Like

This topic was automatically closed 54 days after the last reply. New replies are no longer allowed.

Hi Hadley. Thank you for your reply!

My supervisor told me today the problem lies in the
fit2 <- tslm(train ~ trend + relevel(season, ref = which.min(tapply(train, cycle(train, FUN = sum))))) function. If I change that tofit2 <- tslm(train ~ trend + season) then it works perfectly fine. But that also means there is a bug in there somewhere because it should work whether or not I change the baseline level for the regression model. I also realized I didn't need to have all of my code inside the reactive function since input$dataset is already reactive.

I will definitely have a look into your suggestion as I can already see it is way more efficient and tidied than my original script. Hopefully I will be able to find out where the problem is, so when I change my baseline level, it might then work.

Thanks again!