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

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?

CodePudding user response：

I think you can do this much more easily with nlme::gnls:

nlme::gnls(y~1/(b x^p),
           params = list(p~group-1, b~1), 
           data=d, 
           start = list(b=1, p = rep(3,3)))

Results:

Generalized nonlinear least squares fit
  Model: y ~ 1/(b   x^p) 
  Data: d 
  Log-likelihood: 62.05887

Coefficients:
p.groupA p.groupB p.groupC        b 
2.262383 2.895903 3.475324 1.407561 

Degrees of freedom: 17 total; 13 residual
Residual standard error: 0.007188101

The params argument allows you to specify fixed-effect submodels for each nonlinear parameter. Using p ~ b-1 parameterizes the model with a separate estimate for each group, rather than fitting a baseline (intercept) value for the first group and the differences between successive groups. (I still don't understand the differences between this fit and the nls-based one, I will investigate further.)

CodePudding user response：

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