0

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?

Laurent Franckx
  • 161
  • 1
  • 6

2 Answers2

1

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
Onyambu
  • 67,392
  • 3
  • 24
  • 53
  • I think I have been a bit too quick in thinking the problem was solved. I first thought that it was due to something I was mistaken about, but the same problem actually arises with the example code. To see this, run the regression with nlsLM rather than with nls. If you do this, you will get the error message: "Error in `rownames<-`(`*tmp*`, value = "MyArguments") : length of 'dimnames' [1] not equal to array extent". If you dig inside the code of nls.lm, you will see that only the first element of the list varies from iteration to iteration. – Laurent Franckx Jun 12 '18 at 13:07
  • I think I may have the underlying problem. If you look at https://www.rdocumentation.org/packages/minpack.lm/versions/1.2-1/topics/nls.lm , you see that " If par is a list, then each element must be of length 1". However, if you use nls.lm in debug mode, you see that par[[1]] yields [1] 200.0 0.1. So the problem is indeed that only the first element is read. I am not sure this can be solved without modifying nls.lm. Are there any other alternatives to nls? – Laurent Franckx Jun 13 '18 at 08:46
0

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).

Laurent Franckx
  • 161
  • 1
  • 6