Error in nls(y~ alpha/(1 + exp(x - xO)/beta), : singular gradient

Question

I have the following error when using nls in R. I want to fit a sigmoid function to my data. Y is from 0.01 to 1 and x is form 0.1 to -3. I am using the code as follows:

dummy_df <- data.frame(
  x = seq(0.1, -3, length.out = 50),
  y = seq(0.01, 1, length.out = 50)
)

  fit <- nls(y ~ alpha / (1 + exp(x - xO) / beta),data = dummy_df,
             start = list(alpha = 1.0009, beta = -10, xO = 0),
    control = nls.control(minFactor = 0.00001, maxiter = 10)) 

This results in the following error: error in nls(y ~ alpha/(1 + exp(x - xO)/beta), : singular gradient

Answer 1

There are several problems:

1) the parameters are not identifiable. It is like asking to solve a+b=1 for a and b. You must fix one of them.

We can write the RHS as

alpha / (1 + exp(x - xO) / beta)
= alpha / (1 + exp(x) * exp(-x0) / beta)

and exp(-x0)/beta is really a single parameter, not two. We can address this by fixing either beta or x0 to an arbitrary value. We have used the nls2 "plinear-random" algorithm which allows us to omit starting values that enter linearly, here alpha, and we have set beta to -1 to make the parameters identifiable. With "plinear" based algorithms we specify the RHS as a matrix whose columns are the nonlinear portions whereas the linear parameters are coefficients that multiply the columns and can optionally be specified as the column names. nls2's "plinear-random" algorithm evaluates the RHS at random values within the rectangle described by start and gives the best of these. We can then use that as a starting value for nls.

2) The starting values are bad. Use nls2 as shown below to address that.

3) The model does not describe the data well. A few plots should convince one of that.

Here is the code; however, as mentioned in (3) we are not going to get good results fitting this model to this data although we can get a result without errors.

set.seed(123)


library(nls2)
set.seed(123)
beta <- -1

fit0 <- nls2(y ~ cbind(alpha = 1 / (1 + exp(x - x0) / beta)),
  data = dummy_df,
  start = data.frame(x0 = c(-20, 20)),
  algorithm = "plinear-random")

fit <- nls(y ~ cbind(alpha = 1 / (1 + exp((x - xO) / beta))),
  data = dummy_df, start = coef(fit0)[1], alg = "plinear")

Different formula

@Ben Bolker in comments suggests that perhaps you meant the formula.
alpha/(1 + exp((x - xO)/beta)). If that is the case it would be easier to fit if re-expressed in the equivalent more linear parameterization shown.

library(nls2)
set.seed(123)

fo2 <- y ~ cbind(alpha = 1/(1 + exp(xO.over.beta - one.over.beta * x)))

fit0 <- nls2(fo2, data = dummy_df,
  start = data.frame(one.over.beta = c(-5, 5), xO.over.beta = c(-5 -5)),
  algorithm = "plinear-random")

fit2 <- nls(fo2, data = dummy_df, 
  start = c(one.over.beta = coef(fit0)[[1]], xO.over.beta = coef(fit0)[[2]]),
  algorithm = "plinear")
fit2

alpha <- coef(fit2)[[3]]
beta <- 1/coef(fit2)[["one.over.beta"]]
xO <- beta * coef(fit2)[["xO.over.beta"]]

In that case we can get a more reasonable fit.

plot(y ~ x, dummy_df)
lines(fitted(fit2) ~ x, dummy_df, col = "red")

(continued after graph)

With this new model @Roland suggested in comments using the self starting model SSlogis which gives the same residual sum of squares with somewhat different coefficients.

nls(y ~ SSlogis(x, alpha, x0, beta), data = dummy_df)

Also note that the drc package has numerous canned growth curves that it will fit without needing starting values. For example, this uses the 5 parameter multi2 model from that package.

library(drc)

fm <- drm(y ~ x, data = dummy_df, fct = multi2())

summary(fm)
sum(resid(fm)^2)  # residual sum of squares

Update

Numerous updates including fixing some errors and following up on comments below the question.