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

user_16491398 0
Viewed 213

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?

nls grouping r parameters
2 Answers

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

fit2 <- 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. (In R's formula language, -1 or +0 signify "fit a model without intercept/set the intercept to 0", which in this case corresponds to fitting all three groups separately.)

I'm quite surprised that gnls and nls don't give identical results (although both give reasonable results); would like to dig in further ...

points and fitted curves for nls and gnls

Parameter estimates (code below):

  term    nls  gnls
1 b      1.41  1.40
2 pA     2.28  2.28
3 pB     3.19  3.14
4 pC     3.60  3.51

par(las = 1, bty = "l")
plot(y~x, data = d, col = d$group, pch = 16)
xvec <- seq(0, 1, length = 21)
f <- function(x) factor(x, levels = c("A","B","C"))
## fit1 is nls() fit
ll <- function(g, c = 1) {
  lines(xvec, predict(fit1, newdata = data.frame(group=f(g), x = xvec)), col = c)
}
Map(ll, LETTERS[1:3], 1:3)
d2 <- expand.grid(x = xvec, group = f(c("A","B","C")))
pp <- predict(fit2, newdata = d2)
ll2 <- function(g, c = 1) {
  lines(xvec, pp[d2$group == g], lty = 2, col = c)
}
Map(ll2, LETTERS[1:3], 1:3)
legend("bottomleft", lty = 1:2, col = 1, legend = c("nls", "gnls"))

library(tidyverse)
library(broom)
library(broom.mixed)
(purrr::map_dfr(list(nls=fit1, gnls=fit2), tidy, .id = "pkg")
  %>% select(pkg, term, estimate)
  %>% group_by(pkg)
  ## force common parameter names
  %>% mutate(across(term, ~ c("b", paste0("p", LETTERS[1:3]))))
  %>% pivot_wider(names_from = pkg, values_from = estimate)
)
user_1902770

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
user_164913980
Related