Search code examples
rparametersgroupingnls

Using nls or nlsLM to fit global and group-specific parameters

I would like to use nls to fit a global parameter and group-specific parameters. The closest I have found to a minimum reproducible example is below (found here: https://stat.ethz.ch/pipermail/r-help/2015-September/432020.html)

#Generate some data
d <- transform(data.frame(x=seq(0,1,len=17),
     group=rep(c("A","B","B","C"),len=17)), y =
     round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2))

#Fit to model using nls
nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=rep(3,length(levels(d$group)))))

This gives me an error:

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

I have not been able to figure out if the error is coming from bad guesses for the starting values, or the way this code is dealing with group-specific parameters. It seems the line with p=rep(3,length(levels(d$group))) is for generating c(3,3,3), but switching this part of the code does not remove the problem (same error obtained as above):

#Fit to model using nls
nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3, 3, 3)))

Switching to nlsLM gives a different error which leads be to believe I am having an issue with the group-specific parameters:

#Generate some data
library(minpack.lm)
d <- transform(data.frame(x=seq(0,1,len=17),
                          group=rep(c("A","B","B","C"),len=17)), y =
                 round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2))

#Fit to model using nlsLM
nlsLM(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3,3,3)))

Error:

Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent

Any ideas?

Solution

  • I was able to get this by switching the class of the group from chr to factor. Note the addition of factor() when generating the dataset.

    > d <- transform(data.frame(
+       x=seq(0,1,len=17),
+       group=rep(factor(c("A","B","B","C")),len=17)),
+       y=round(1/(1.4+x^ifelse(group=="A", 2.3, ifelse(group=="B",3.1, 3.5))),2)
+     )
> str(d)
'data.frame':   17 obs. of  3 variables:
 $ x    : num  0 0.0625 0.125 0.1875 0.25 ...
 $ group: Factor w/ 3 levels "A","B","C": 1 2 2 3 1 2 2 3 1 2 ...
 $ y    : num  0.71 0.71 0.71 0.71 0.69 0.7 0.69 0.69 0.62 0.64 ...
> nls(y~1/(b+x^p[group]), data=d, start=list(b=1, p=c(3,3,3)))
Nonlinear regression model
  model: y ~ 1/(b + x^p[group])
   data: d
    b    p1    p2    p3 
1.406 2.276 3.186 3.601 
 residual sum-of-squares: 9.537e-05

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