Non-linear regression, nls, in R: singular gradient

Question

I want to fit my data to a specific function that has already been optimized using Matlab.

I get the following error: 'Warning message: Computation failed in stat_smooth(): singular gradient '

Please help! Here is my R code:

tibble
       x     y     SEM
 1     1 0.0342 0.00532
 2     3 0.0502 0.00639
 3     5 0.0700 0.0118 
 4    10 0.123  0.0269 
 5    20 0.154  0.0125 
 6    30 0.203  0.0190 
 7    40 0.257  0.0255 
 8    50 0.287  0.0266 
 9    60 0.345  0.0347 
10    90 0.442  0.0398 
11   120 0.569  0.0570 
12   180 0.726  0.0406 
13   240 0.824  0.0150 
14   360 0.868  0.00821
15  1440 0.890  0.0246 

tibble %>% 
  ggplot(aes(x, y)) +
  geom_point()+
  geom_errorbar(aes(ymin=y-SEM, ymax=y+SEM), width=25)+
  geom_ribbon(aes(ymin = y-2.575*SEM, ymax = y+2.575*SEM), alpha = 0.1)+
  geom_smooth(method="nls", 
              formula= y ~ (1-((k2/(k2-k1))*exp(-k1*x))+((k1/(k2-k1))*exp(-k2*x))),
              se=F,
              method.args = list(start=list(k1=0.006999, k2=849.6)))

You might want to try a simple exponential model instead. Fit looks pretty good on the plot: `fm <- nls(y ~ a * (1 - exp(-b * x)), DF, start = list(a = 1, b = 1)); plot(DF[1:2]); lines(fitted(fm) ~ x, DF)` — G. Grothendieck, Nov 13 '21 at 14:41
@G.Grothendieck the formula actually simplifies to `1 - exp(-k1 * x)` as `k2` tends to infinity. The sum of squares also falls to an asymptote as `k2` tends to infinity, so you will always get a better fit by just using `1 - exp(-k1 * x)`. Of course, you could add in the `a` parameter as per your suggestion to get an even better fit if this makes sense for the model. I've updated my answer to reflect this. — Allan Cameron, Nov 13 '21 at 15:27
`fm2 <- nls(y ~ (1 - exp(-b * x)), DF, start = list(b = coef(fm)[2])); anova(fm2, fm)` indicates that adding `a` to the model has a p value of 0.001106. — G. Grothendieck, Nov 13 '21 at 15:38
@G.Grothendieck yes, I know it's a better fit, and I demonstrate that graphically in my answer. I meant that `1 - exp(-k1 * x)` is necessarily better than the model in the OP,s question. The model `a * (1 - exp(-b * x))` will always be a better fit than `1 - exp(-k1 * x)` (or at worst, equal). My point is that we don't know whether a model of the underlying process would be free to set a value `a`, or whether we need to assume the asymptote is 1 on theoretical grounds. At least both options are now open to the OP. — Allan Cameron, Nov 13 '21 at 16:59
It's not necessarily true that it would necessarily fit *significantly* better though but it does. — G. Grothendieck, Nov 13 '21 at 17:36

Allan Cameron · Accepted Answer · 2021-11-13T15:28:36.613

As it stands, nls can't decide on an optimal value of k2 because the sum of squares decreases asymptotically towards a value of about 0.0225 as k2 tends to infinity. There is therefore no finite value of k2 that minimizes the sum of squares. As k2 tends to infinity, it effectively cancels out, to leave the formula equivalent to:

y ~ 1 - exp(-k1*x)

which means that this formula will always fit the data better than the original formula for any finite value of k2.

In short, k2 is a redundant parameter that can only worsen your fit.

Your plot could therefore be made like this:

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_errorbar(aes(ymin = y - SEM, ymax = y + SEM), width = 25) +
  geom_ribbon(aes(ymin = y - 2.575 * SEM, ymax = y + 2.575 * SEM), alpha = 0.1) +
  geom_smooth(method = "nls", 
              formula = y ~ 1 - exp(-k1 * x),
              se = FALSE,
              method.args = list(start = list(k1 = 0.006999)))

Or, as G. Grothendieck suggests, use an extra parameter to optimize the fit like this (assuming it makes physical sense in your use case)

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_errorbar(aes(ymin = y - SEM, ymax = y + SEM), width = 25) +
  geom_ribbon(aes(ymin = y - 2.575 * SEM, ymax = y + 2.575 * SEM), alpha = 0.1) +
  geom_smooth(method = "nls", 
              formula = y ~ k2 * (1 - exp(-k1 * x)),
              se = FALSE,
              method.args = list(start = list(k1 = 0.006999, k2 = 1)))

Non-linear regression, nls, in R: singular gradient

1 Answers1