error messages fitting a non-linear exponential model between two variables

Question

I have two variables that I'm trying to model the relationship between and extract the residuals. The relationship between the two variables is clearly a non-linear exponential relationship. I've tried a few different approaches with nls, but I keep getting different error messages.

# dataset
df <- structure(list(y = c(464208.56, 334962.43, 361295.68, 426535.68, 258843.93, 272855.46, 
   166322.72, 244695.28, 227003.03, 190728.4, 156025.45, 72594.24, 56911.4, 175328.95, 161199.76, 
   152520.77, 190610.57, 60734.34, 31620.9, 74518.86, 45524.49, 2950.58, 2986.38, 15961.77, 12484.05, 
   6828.41, 2511.72, 1656.12, 5271.4, 7550.66, 3357.71, 3620.43, 3699.85, 3337.56, 4106.55, 3526.66, 
   2996.79, 1649.89, 4561.64, 1724.25, 3877.2, 4426.69, 8557.61, 6021.61, 6074.17, 4072.77, 4032.95, 
   5280.16, 7127.22), 
   x = c(39.23, 38.89, 38.63, 38.44, 38.32, 38.27, 38.3, 38.4, 38.56, 38.79, 39.06, 39.36, 39.68, 
   40.01, 40.34, 40.68, 41.05, 41.46, 41.93, 42.48, 43.14, 43.92, 44.84, 45.9, 47.1, 48.4, 49.78, 
   51.2, 52.62, 54.01, 55.31, 56.52, 57.6, 58.54, 59.33, 59.98, 60.46, 60.78, 60.94, 60.92, 60.71, 
   60.3, 59.69, 58.87, 57.86, 56.67, 55.33, 53.87, 52.33)), 
   row.names = c(NA, -49L), 
   class = c("tbl_df", "tbl", "data.frame"), 
   na.action = structure(c(`1` = 1L, `51` = 51L), 
   class = "omit"))

# initial model
m <- nls(y ~  a * exp(r * x), 
         start = list(a = 0.5, r = -0.2), 
         data = df)
Error in nls(y ~ a * exp(r * x), start = list(a = 0.5, r = -0.2), data = df,  : singular gradient

# add term for alg
m <- nls(y ~  a * exp(r * x), 
         start = list(a = 0.5, r = -0.2), 
         data = df,
         alg = "plinear")
Error in nls(y ~ a * exp(r * x), start = list(a = 0.5, r = -0.2), data = df,  : 
  step factor 0.000488281 reduced below 'minFactor' of 0.000976562

I would note that for the equation y = a * exp(r * x), the way to think about a is the y value when x = 0. Looking at your plot, when x is around 40, the y values are around 300,000, and those values would be exponentially bigger as x approaches 0, so a = 0.5 is a pretty terrible starting value, off by many orders of magnitude. In fact, since your data doesn't get anywhere close to x = 0, a shifted exponential model might be a much better fit, maybe y ~ a * exp(r * (x - b)) with a possible start values of a = 5e5 and b = 38. (Though I still get a singular gradient there) — Gregor Thomas
– Gregor Thomas, Commented Jun 1, 2023 at 17:55
Is there any reason not to linearize your model by modeling log(y)? lm(log(y) ~ x, data = df) fits just fine, and you can exponentiate the intercept to get a and use the x coefficient as r (which is -0.1915, very close to the r = -0.2 starting value you use). The error structure is a little different, not sure if you care about that... The intercept is 19.1994, so exp(19.1994) = 217868015, which is off from you 0.5 by about 10^9 times. — Gregor Thomas
– Gregor Thomas, Commented Jun 1, 2023 at 18:07
Thanks @GregorThomas. I was trying to avoid a transformation that I'd have to backtransform at a later data. Also, when I plotted a log(y), there's still a strong non-linear pattern in the data. — tnt
– tnt, Commented Jun 1, 2023 at 18:43
Not logging the data doesn't get rid of the nonlinear pattern in the logged data, it just makes it harder to see. If you want more accuracy you could fit log(y) ~ poly(x, 2) .... — Ben Bolker
– Ben Bolker, Commented Jun 1, 2023 at 19:04
@BenBolker, thanks for this suggestion. I ended up using this option. — tnt
– tnt, Commented Jun 2, 2023 at 20:38

Ben Bolker · Accepted Answer · 2023-06-01 18:22:07Z

3

log-Gaussian GLM

As @Gregor Thomas suggests you could linearize your problem (fit a log-linear regression), at the cost of changing the error model. (Basic model diagnostics, i.e. a scale-location plot, suggest that this would be a much better statistical model!) However, you can do this efficiently without changing the error structure by fitting a log-link Gaussian GLM:

m1 <- glm(y ~ x, family = gaussian(link = "log"), data = df)

The model is y ~ Normal(exp(b0 + b1*x), s), so a = exp(b0), r = b1.

I tried using list(a=exp(coef(m1)[1]), r=coef(m1)[2]) as starting values, but even this was too finicky for nls().

There are two ways to get nls to work.

shifted exponential

As @GregorThomas suggests, shifting the x-axis to x=38 also works fine (given a sensible starting value):

m <- nls(y ~  a * exp(r * (x-38)), 
         start = list(a = 3e5, r = -0.35), 
         data = df)

provide nls with a gradient

The deriv function will generate a function with the right structure for nls (returns the objective function, with a ".grad" attribute giving a vector of derivatives) if you ask it nicely. (I'm also using the exponentiated intercept from the log-Gaussian GLM as a starting value ...)

f <- deriv( ~ a*exp(r*x), c("a", "r"), function.arg = c("x", "a", "r"))
m2 <- nls(y ~  f(x, a, r),
         start = list(a = exp(coef(m1)[1]), r = -0.35),
         data = df)

We can plot these to compare the predictions (visually identical):

par(las = 1, bty = "l")
xvec <- seq(38, 60, length = 101)
plot(y ~ x, df)
lines(xvec, predict(m1, newdata = data.frame(x=xvec), type = "response"),
      col = 2)
lines(xvec, predict(m, newdata = data.frame(x=xvec)), col = 4,  lty = 2)
lines(xvec, predict(m2, newdata = data.frame(x=xvec)), col = 5,  lty = 2)

With a little bit of extra work (exponentiating the intercept for the Gaussian GLM, shifting the x-origin back to zero for the nls fit) we can compare the coefficients (only equal up to a tolerance of 2e-4 but that should be good enough, right?)

a1 <- exp(coef(m1)[[1]])
a2 <- coef(m)[[1]]*exp(-38*coef(m)[[2]])
all.equal(c(a = a1, r = coef(m)[[2]]),
          c(a = a2, r = coef(m1)[[2]]), tolerance = 1e-4)
all.equal(c(a = a1, r = coef(m)[[2]]),
          coef(m2), tolerance = 2e-4)

edited Jun 1, 2023 at 18:22

answered Jun 1, 2023 at 18:13

Ben Bolker

230k26 gold badges405 silver badges497 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

tnt Over a year ago

Thanks for answer! It's really helpful to see that all three of these different approaches give a similar answer. Since I'm planning to extract the residual values, would there be a preference between any of the models?

Ben Bolker Over a year ago

At a quick check, they all give identical residuals as well. As I mentioned in my answer, if you care about the statistical properties of this fit (e.g. homoscedasticity), then you should definitely consider fitting lm(log(y) ~ x, df) instead ...

tnt Over a year ago

Thanks Ben. Yes, the heteroskedasticity is quite obvious in all these options.

Collectives™ on Stack Overflow

error messages fitting a non-linear exponential model between two variables

1 Answer 1

log-Gaussian GLM

shifted exponential

provide nls with a gradient

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

log-Gaussian GLM

shifted exponential

provide nls with a gradient

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related