using R function nls with long list of variables

Question

I am trying to estimate a non-linear model with the R function nls.

I will illustrate my problem with an example from the "help" facility for nls.

Treated <- Puromycin[Puromycin$state == "treated", ]
weighted.MM <- function(resp, conc, Vm, K)
{
    pred <- (Vm * conc)/(K + conc)
    (resp - pred) / sqrt(pred)
}

Pur.wt <- nls( ~ weighted.MM(rate, conc, Vm, K), data = Treated,
          start = list(Vm = 200, K = 0.1))

In this example, weighted.MM is a function with a very limited number of arguments, and I have no problem implementing analogous approaches for the type of model I am working with.

However, I am trying now to move to a more realistic problem, where I have literally dozens of arguments to pass on to the function. I could of course simply enumerate them, but I find this a bit unwieldy.

I have considered putting them in a separate list first. For instance, using the example from above, I would first define:

MyArguments <- list(Vm, K)  

and then passing MyArguments on as argument to the function (and accessing the individual arguments from within the function). However, that doesn't work, because I get the error message

Error: object 'Vm' not found

Alternatively,

    MyArguments <- list("Vm", "K") 

weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
{
  conc <- c(conc1, conc.1)
  pred <- (thearguments[[1]] * conc)/(thearguments[[2]] + conc)
  (resp - pred) / sqrt(pred)
}

Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
                data = lisTreat, start = list(Vm = 200, K = 0.1))

yields:

Error in thearguments[[1]] * conc : 
  non-numeric argument to binary operator
Called from: weighted.MM1(rate, conc1, conc.1, MyArguments)

Is there any work-around for this?

Answer 1

I have tried to implement a new solution based on the suggestions kindly made by @Onyambu.

This has created new problems, though.

First, I have tried to implement a solution with nls. Here's the actual code that I have used:

DiffModel <- nls(COPERTFreq ~ CalculateProbaVehChoiceDiffusion(MyArguments, years = RegistYear),
                 data = DataSetForModel , start = list(MyArguments = c(ASC_Mat, ASC_No_Size)))

where CalculateProbaVehChoiceDiffusion() is a non-linear function that was defined elsewhere, RegistYear is a constant, and MyArguments is the list of coefficients to estimate, with c(ASC_Mat, ASC_No_Size) as initial values.

This leads to the following error message:

Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

Now, I had read elsewhere that this problem can be solved by using nlsLM instead. This leads to a new error message:

Error in `rownames<-`(`*tmp*`, value = "MyArguments") : 
  length of 'dimnames' [1] not equal to array extent

OK, so I ran the model again with nls.lm in debug mode. This shows that the error message originates from the following line of code:

names(out$par) <- rownames(out$hessian) <- colnames(out$hessian) <- names(out$diag) <- names(par)

However, it is from inspecting the "out" object that becomes clear where the problem lies. First, out$hessian is simply a scalar, while I would expect its number of rows and columns to be equal to the number of parameters. Second, out$par$MyArguments shows that, except for the first element, the value of MyArguments does not change from one iteration to the other.

Is this a known bug or do I have to modify the way I pass on MyArguments to the function call?

Note that, as far as I can see, this problems also occurs when I apply nlsLm to the example provided by @Onyambu:

> undebug(nls.lm)
> Treated <- Puromycin[Puromycin$state == "treated", ]
> 
> lisTreat <- with(Treated,
+                  list(conc1 = conc[1], conc.1 = conc[-1], rate = rate))
> 
> weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
+ {
+   conc <- c(conc1, conc.1)
+   pred <- (thearguments[1] * conc)/(thearguments[2] + conc)
+   (resp - pred) / sqrt(pred)
+ }
> nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
+      data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Nonlinear regression model
  model: 0 ~ weighted.MM1(rate, conc1, conc.1, MyArguments)
   data: lisTreat
MyArguments1 MyArguments2 
   206.83468      0.05461 
 residual sum-of-squares: 14.6

Number of iterations to convergence: 5 
Achieved convergence tolerance: 3.858e-06
> 
> nlsLM( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
+        data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Error in `rownames<-`(`*tmp*`, value = "MyArguments") : 
  length of 'dimnames' [1] not equal to array extent

Thus, while nls works with this example, nlsLM does not, and yields the same error message as with my code. I then run nlsLM again with nls.lm in debug mode. After the following line,

out <- .Call("nls_lm", par, lower, upper, fn1, jac1, ctrl, 
    new.env(), PACKAGE = "minpack.lm")

I inspect the out object and see:

$par
$par$MyArguments
[1] 244.5117   0.1000


$hessian
[1] 0.02859739

$fvec
 [1] -5.5215487 -0.9787468 -0.5543382 -1.5986605  0.4486608 -0.9651245  0.7020058  1.2419040  1.1430780  0.4488084  1.1445818  1.6121474

$info
[1] 1

$message
[1] "Relative error in the sum of squares is at most `ftol'."

$diag
$diag[[1]]
[1] 0.2077949


$niter
[1] 4

$rsstrace
[1] 112.59784  43.41211  42.89350  42.89349  42.89349

Thus, the value of the second argument has not changed after 4 iterations. This could of course be the correct solution. But I find it an intriguing coincidence that the same happens with my model.

Final edit: I have finally decided to fall back on the complete enumeration of all the arguments. As I wrote in the statement of the problem, it's not very elegant, but, at least it works with nls.lm (though still not with nls).

Answer 2

Since you are looking for the values of your arguments, you do not need to define them: Look at the code below:

weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
{
  conc <- c(conc1, conc.1)
  pred <- (thearguments[1] * conc)/(thearguments[2] + conc)
  (resp - pred) / sqrt(pred)
}
nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
                data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Nonlinear regression model
  model: 0 ~ weighted.MM1(rate, conc1, conc.1, MyArguments)
   data: lisTreat
MyArguments1 MyArguments2 
   206.83468      0.05461 
 residual sum-of-squares: 14.6

Number of iterations to convergence: 5 
Achieved convergence tolerance: 3.858e-06