why a "subscript out of bounds" error in Shiny, but not R?

Question

I recently posted a similar inquiry in the shiny google group, but did not find a solution. We are developing a Shiny app and as the subject indicates we get an "error: subscript out of bounds" message upon running the app. However, when we isolate the offending code and run it on its own in RStudio, there is no error.

This makes me wonder if there is a bug in Shiny itself, or if we are missing something.

Please see the instructions below along with a small example that produces the error. We are using Shiny version 0.8.0 and RStudio 0.98.501.

Thanks for your help!

---

To run the app, place ui.R and server.R (see below) in a folder and run

library(shiny)
runApp("<folder path>")

It should produce a user interface with a button on the left, but on the right you will see "error: subscript out of bounds".

However, if just run the following three lines of code (approximately lines 57-59 in server.R)

show=data.frame(ps=c(4,-1,0,1),ns=c(0,1,0,0),ts=c(45842,15653,28535,21656))
best.fit1=regsubsets(ts~ps+ns,data=show,nvmax=1)
pred1=predict.regsubsets(best.fit1,show,id=1)  # line that offends Shiny

in RStudio (need to include the function "predict.regsubsets" - given at the beginning of server.R), then there are no errors.

#####################
## server.R
#####################

library(rms)
library(leaps)
library(shiny)
library(datasets)
library(stringr)
library(ttutils)
library(plyr)
library(utils)
library(ggplot2)

# object is a regsubsets object
# newdata is of the form of a row or collection of rows in the dataset
# id specifies the number of terms in the model, since regsubsets objects 
#  includes models of size 1 up to a specified number
predict.regsubsets=function(object,newdata,id,...){
  form=as.formula(object$call[[2]])

  mat=model.matrix(form,newdata)

  mat.dims=dim(mat)
  coefi=coef(object,id=id)
  xvars=names(coefi)
  # because mat only has those categorical variable categories associated with newdata, 
  # it is possible that xvars (whose variables are defined by the "best" model of size i)
  # has a category that is not in mat
  diffs=setdiff(xvars,colnames(mat))
  ndiffs=length(diffs)
  if(ndiffs>0){
    # add columns of 0's for each variable in xvars that is not in mat
    mat=cbind(mat,matrix(0,mat.dims[1],ndiffs))
    # for the last "ndiffs" columns, make appropriate names
    colnames(mat)[(mat.dims[2]+1):(mat.dims[2]+ndiffs)]=diffs
    mat[,xvars]%*%coefi
  }
  else{
    mat[,xvars]%*%coefi
  }
}

# Define server logic required to summarize and view the selected dataset
shinyServer(function(input, output) {

mainTable1 <- reactive({

  }) 

output$table21 <- renderTable({
    mainTable1()
  })


formulamodel1 <- reactive({
    #ticketsale<-dataset1Input()

  show=data.frame(ps=c(4,-1,0,1),ns=c(0,1,0,0),ts=c(45842,15653,28535,21656))
  best.fit1=regsubsets(ts~ps+ns,data=show,nvmax=1)
  pred1=predict.regsubsets(best.fit1,show,id=1)

  })

output$model1fit <- renderPrint({
    formulamodel1()

  }) 

 })

######################
## end server.R
######################

######################
## ui.R
######################

library(rms)
library(leaps)
library(shiny)
library(datasets)
library(stringr)
library(ttutils)
library(plyr)
library(utils)
library(ggplot2)

shinyUI(pageWithSidebar(

headerPanel("Forecasting ticket sales for xxx"),

sidebarPanel(
        p(strong("Model Fitting")),

    selectInput("order1", "Sort results by:",c("a","b","c")),
    submitButton("Run Model")

    ),

   mainPanel(

    h3(strong("Model fit without using ticket sales") ),
    tableOutput("table21"),
    verbatimTextOutput(outputId = "model1fit")

   )
))

Answer 1

These three lines only seem to work when executed in the global environment. If you take that snippet and run it inside of a local({...}) block you'll see the same error.

The error is coming from the first line of predict.regsubsets, where you look at object$call[[2]]. It's object$call that is very different depending on whether it's executed in the global environment or not; it's created in leaps:::regsubsets.formula by calling sys.call(sys.parent()). Perhaps this needs to be sys.call(sys.parent(0)) (just a guess)?

Answer 2

Thanks to John Harrison for this answer. He attempted to reply via the shiny Google group but the system deleted his answers, as well as my attempt later to post his solution. Here it is.

---

John Harrison says:

The issue is with the regsubsets function:

> test_env <- new.env(parent = globalenv())
> with(test_env, {show=data.frame(ps=c(4,-1,0,1),ns=c(0,1,0,0),ts=c(45842,15653,28535,21656))
+                 best.fit1=regsubsets(ts~ps+ns,data=show,nvmax=1)
+                 #pred1=predict.regsubsets(best.fit1,show,id=1)
+                 #pred1
+                 best.fit1})
Subset selection object
Call: eval(expr, envir, enclos)
2 Variables  (and intercept)
   Forced in Forced out
ps     FALSE      FALSE
ns     FALSE      FALSE
1 subsets of each size up to 1
Selection Algorithm: exhaustive

You can see it gets it Call: output relative to the environment its in:

> getAnywhere(regsubsets.formula)
A single object matching ‘regsubsets.formula’ was found
It was found in the following places
  registered S3 method for regsubsets from namespace leaps
  namespace:leaps
with value

function (x, data, weights = NULL, nbest = 1, nvmax = 8, force.in = NULL, 
    force.out = NULL, intercept = TRUE, method = c("exhaustive", 
        "backward", "forward", "seqrep"), really.big = FALSE, 
    ...) 
{
    formula <- x
    rm(x)
    mm <- match.call()
    mm$formula <- formula
    mm$x <- NULL
    mm$nbest <- mm$nvmax <- mm$force.in <- mm$force.out <- NULL
    mm$intercept <- mm$method <- mm$really.big <- NULL
    mm[[1]] <- as.name("model.frame")
    mm <- eval(mm, sys.frame(sys.parent()))
    x <- model.matrix(terms(formula, data = data), mm)[, -1]
    y <- model.extract(mm, "response")
    wt <- model.extract(mm, "weights")
    if (is.null(wt)) 
        wt <- rep(1, length(y))
    else wt <- weights
    a <- leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, 
        force.in = force.in, force.out = force.out, intercept = intercept)
    rval <- switch(1 + pmatch(method[1], c("exhaustive", "backward", 
        "forward", "seqrep"), nomatch = 0), stop(paste("Ambiguous or unrecognised method name :", 
        method)), leaps.exhaustive(a, really.big), leaps.backward(a), 
        leaps.forward(a), leaps.seqrep(a))
    rval$call <- sys.call(sys.parent())
    rval
}
<environment: namespace:leaps>

rval$call <- sys.call(sys.parent())

is the offending line of code

---

I replied:

I'm in a bit over my head in terms of these R functions, environments, etc. I roughly followed your explanation above but I don't understand it enough to have any real sort of idea of what to do to fix it (or whether it is even fixable). Could you easily point me in the right direction?

---

John replied:

You could define your own regsubsets function:

myregsubsets <- function (x, data, weights = NULL, nbest = 1, nvmax = 8, force.in = NULL, 
                          force.out = NULL, intercept = TRUE, method = c("exhaustive", 
                                                                         "backward", "forward", "seqrep"), really.big = FALSE, 
                          ...){
  formula <- x
  rm(x)
  mm <- match.call()
  mm$formula <- formula
  mm$x <- NULL
  mm$nbest <- mm$nvmax <- mm$force.in <- mm$force.out <- NULL
  mm$intercept <- mm$method <- mm$really.big <- NULL
  mm[[1]] <- as.name("model.frame")
  mm <- eval(mm, sys.frame(sys.parent()))
  x <- model.matrix(terms(formula, data = data), mm)[, -1]
  y <- model.extract(mm, "response")
  wt <- model.extract(mm, "weights")
  if (is.null(wt)) 
    wt <- rep(1, length(y))
  else wt <- weights
  a <- leaps:::leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, 
                           force.in = force.in, force.out = force.out, intercept = intercept)
  rval <- switch(1 + pmatch(method[1], c("exhaustive", "backward", 
                                         "forward", "seqrep"), nomatch = 0), stop(paste("Ambiguous or unrecognised method name :", 
                                                                                        method)), leaps:::leaps.exhaustive(a, really.big), leaps:::leaps.backward(a), 
                 leaps:::leaps.forward(a), leaps:::leaps.seqrep(a))
  rval$call <- sys.call(sys.parent())
  rval$x <- formula
  rval
}

predict.regsubsets=function(object,newdata,id,...){
  form=as.formula(object$x)

  mat=model.matrix(form,newdata)

  mat.dims=dim(mat)
  coefi=coef(object,id=id)
  xvars=names(coefi)
  # because mat only has those categorical variable categories associated with newdata, 
  # it is possible that xvars (whose variables are defined by the "best" model of size i)
  # has a category that is not in mat
  diffs=setdiff(xvars,colnames(mat))
  ndiffs=length(diffs)
  if(ndiffs>0){
    # add columns of 0's for each variable in xvars that is not in mat
    mat=cbind(mat,matrix(0,mat.dims[1],ndiffs))
    # for the last "ndiffs" columns, make appropriate names
    colnames(mat)[(mat.dims[2]+1):(mat.dims[2]+ndiffs)]=diffs
    mat[,xvars]%*%coefi
  }
  else{
    mat[,xvars]%*%coefi
  }
}

---

Later, John added:

The regsubsets function assumed the user was calling it in a certain fashion. The myregsubsets is a replacement for regsubsets.formula. In your predict.regsubsets you access the formula using as.formula(object$call[[2]]). When nested in environments this doesnt give you what is expected. The myregsubsets replacement gets the formula using rval$x <- formula. The changed predict.regsubsets then uses form=as.formula(object$x) rather then as.formula(object$call[[2]]).

---