CODEWITHSUNDEEP
rnls
Tyler
Tyler

Reputation: 10032

Scope of nls selfStart functions

I'm building a package that will:

  1. read in a data set
  2. determine which components are needed to model the data1
  3. build a self-starting NLS model using those components
  4. apply that model to the data

The result is stored in a list that contains the data, model, and nls result. Everything works until step 4. I've stored my constructed models in the list, and R no longer recognizes them as selfStart functions.

Here's a toy example illustrating my problem. The following function (from the SSlogis manual) works fine at top-level:

 Chick.1 <- ChickWeight[ChickWeight$Chick == 1, ]
 fm1 <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = Chick.1)
 fm1

Nonlinear regression model
model: weight ~ SSlogis(Time, Asym, xmid, scal) data: Chick.1
Asym xmid scal
937.02 35.22 11.41
residual sum-of-squares: 76.66

Number of iterations to convergence: 0
Achieved convergence tolerance: 7.778e-07

But when I store the function and data together in a list, R doesn't see the selfstarting function as selfstarting anymore:

myobj <- list()
myobj$model <- SSlogis
myobj$data <- ChickWeight[ChickWeight$Chick == 1, ]

nls(weight ~ myobj$model(Time, Asym, xmid, scal), data = myobj$data)

Error in getInitial.default(func, data, mCall = as.list(match.call(func, :
no 'getInitial' method found for "function" objects

My workflow will eventually include processing dozens of datasets, so I'd like to keep each dataset and its associated model in its own object (and these objects will probably end up in a list of objects). Is there some way to maintain or restore the environment of my selfStart functions even after they're stored inside another list?

UPDATE

In response to Gregor's suggestion, I tried this:

nls(as.formula(sprintf("weight ~ %s(Time, Asym, xmid, scal)",
                       "myobj$model")), data = myobj$data)

Error in nls(as.formula(sprintf("weight ~ %s(Time, Asym, xmid, scal)", :
singular gradient In addition: Warning message:
In nls(as.formula(sprintf("weight ~ %s(Time, Asym, xmid, scal)", :
No starting values specified for some parameters.
Initializing ‘Asym’, ‘xmid’, ‘scal’ to '1.'.
Consider specifying 'start' or using a selfStart model

Update 2

Inspired by @Gregor I came up with a workaround:

nlsDispatch <- function(obj){
  GLOBAL_NLS <<- obj$model
  nls(weight ~ GLOBAL_NLS(Time, Asym, xmid, scal), data = myobj$data) 
}

nlsDispatch(myobj)

Nonlinear regression model
model: weight ~ GLOBAL_NLS(Time, Asym, xmid, scal)
data: myobj$data
Asym xmid scal
937.02 35.22 11.41
residual sum-of-squares: 76.66

Number of iterations to convergence: 0
Achieved convergence tolerance: 6.621e-07

That works, but dropping my function into the global environment is ugly. And it suggests that if I had a better handle on environments I should be able to avoid abusing the global environment to get this done.

1: In my application, this is mostly a matter of counting peaks, and figuring out how many Normal curves are needed to model them.

Upvotes: 1

Views: 2597

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 270075

1) This works:

Model <- myobj$model
nls(weight ~ Model(Time, Asym, xmid, scal), data = myobj$data)

giving

Nonlinear regression model
  model: weight ~ Model(Time, Asym, xmid, scal)
   data: myobj$data
  Asym   xmid   scal 
937.02  35.22  11.41 
 residual sum-of-squares: 76.66

Number of iterations to convergence: 0 
Achieved convergence tolerance: 6.621e-07

2) It does seem scoping is messed up for self starting functions. The problem appears in getInitial.formula which uses this line:

func <- get(as.character(object[[2L]][[1L]]))

Note that there is no second argument to get (the environment) so it does not pay attention to environments.

Thus if you want to put the solution above in a function then it you will need to work around it like this:

f <- function() {
  myobj <- list()
  myobj$model <- SSlogis
  myobj$data <- ChickWeight[ChickWeight$Chick == 1, ]
  Model <<- myobj$model
  nls(weight ~ Model(Time, Asym, xmid, scal), data = myobj$data)
}
f()

3) Another workaound is to attach Model. That would keep it out of the global workspace. Note that we detach it afterwards so the function does not leave any trace.

f2 <- function() {
  on.exit(detach())
  if (any(grepl("list", search()))) stop("list already on search path")
  myobj <- list()
  myobj$model <- SSlogis
  myobj$data <- ChickWeight[ChickWeight$Chick == 1, ]
  attach(list(Model = myobj$model))
  nls(weight ~ Model(Time, Asym, xmid, scal), data = myobj$data)
}
f2()

Note that if the detach does not occur then the next time an attach is done there would be two items attached on the search list (whereas in the prior approach it always overwrites the global variable so that can't happen) so we check that there is no such list on the search path and stop with an error if there is.

Upvotes: 2

Related Questions