14
$\begingroup$

I have the following data and would like to fit a negative exponential growth model to it:

Days <- c( 1,5,12,16,22,27,36,43)
Emissions <- c( 936.76, 1458.68, 1787.23, 1840.04, 1928.97, 1963.63, 1965.37, 1985.71)
plot(Days, Emissions)
fit <- nls(Emissions ~ a* (1-exp(-b*Days)), start = list(a = 2000, b = 0.55))
curve((y = 1882 * (1 - exp(-0.5108*x))), from = 0, to =45, add = T, col = "green", lwd = 4)

The code is working and a fitting line is plotted. However, the fit is visually not ideal, and the residual sum of squares seems to be quite huge (147073).

How can we improve our fit? Does the data allow a better fit at all?

We could not find a solution to this challenge on the net. Any direct help or linkage to other websites/posts is greatly appreciated.

$\endgroup$
  • 1
    $\begingroup$ In this case, if you consider a regression model $\text{Emissions}_i=f(\text{Days}_i,a,b)+\epsilon_i$, where $\epsilon_i\sim N(0,\sigma)$, then you obtain similar estimators. By plotting the confidence regions, one can observe how these values are contained in the confindence regions. You cannot expect a perfect fit unless you interpolate the points or use a more flexible nonlinear model. $\endgroup$ – user10525 Jul 23 '12 at 13:12
  • $\begingroup$ I changed the title because "negative exponential model" means something different than described in the question. $\endgroup$ – whuber Jul 23 '12 at 13:15
  • $\begingroup$ Thanks for making the question clearer (@whuber) and thanks for your answer (@Procrastinator). How can I calculate and plot the confidence regions. And, what would be a more flexible nonlinear model? $\endgroup$ – Strohmi Jul 23 '12 at 13:18
  • 4
    $\begingroup$ You need an additional parameter. See what happens with fit <- nls(Emissions ~ a* (1- u*exp(-b*Days)), start = list(a = 2000, b = 0.1, u=.5)); beta <- coefficients(fit); curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T). $\endgroup$ – whuber Jul 23 '12 at 13:20
  • 1
    $\begingroup$ @whuber - maybe you should post that as an answer? $\endgroup$ – jbowman Jul 23 '12 at 13:25
16
$\begingroup$

A (negative) exponential law takes the form $y=-\exp(-x)$. When you allow for changes of units in the $x$ and $y$ values, though, say to $y = \alpha y' + \beta$ and $x = \gamma x' + \delta$, then the law will be expressed as

$$\alpha y' + \beta = y = -\exp(-x) = -\exp(-\gamma x' - \delta),$$

which algebraically is equivalent to

$$y' = \frac{-1}{\alpha} \exp(-\gamma x' - \delta) - \beta = a\left(1 - u\exp(-b x')\right)$$

using three parameters $a = -\beta/\alpha$, $u = 1/(\beta\exp(\delta))$, and $b = \gamma$. We can recognize $a$ as a scale parameter for $y$, $b$ as a scale parameter for $x$, and $u$ as deriving from a location parameter for $x$.

As a rule of thumb, these parameters can be identified at a glance from the plot:

  • The parameter $a$ is the value of the horizontal asymptote, a little less than $2000$.

  • The parameter $u$ is the relative amount the curve rises from the origin to its horizontal asymptote. Here, the the rise therefore is a little less than $2000 - 937$; relatively, that's about $0.55$ of the asymptote.

  • Because $\exp(-3) \approx 0.05$, when $x$ equals three times the value of $1/b$ the curve should have risen to about $1-0.05$ or $95\%$ of its total. $95\%$ of the rise from $937$ to almost $2000$ places us around $1950$; scanning across the plot indicates this took $20$ to $25$ days. Let's call it $24$ for simplicity, whence $b \approx 3/24 = 0.125$. (This $95\%$ method to eyeball an exponential scale is standard in some fields that use exponential plots a lot.)

Let's see what this looks like:

plot(Days, Emissions)
curve((y = 2000 * (1 - 0.56 * exp(-0.125*x))), add = T)

Eyeball fit

Not bad for a start! (Even despite typing 0.56 in place of 0.55, which was a crude approximation anyway.) We can polish it with nls:

fit <- nls(Emissions ~ a * (1- u * exp(-b*Days)), start=list(a=2000, b=1/8, u=0.55))
beta <- coefficients(fit)
plot(Days, Emissions)
curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T, col="Green", lwd=2)

NLS fit

The output of nls contains extensive information about parameter uncertainty. E.g., a simple summary provides standard errors of estimates:

> summary(fit)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 1.969e+03  1.317e+01  149.51 2.54e-10 ***
b 1.603e-01  1.022e-02   15.69 1.91e-05 ***
u 6.091e-01  1.613e-02   37.75 2.46e-07 ***

We can read and work with the entire covariance matrix of the estimates, which is useful for estimating simultaneous confidence intervals (at least for large datasets):

> vcov(fit)
             a             b             u
a 173.38613624 -8.720531e-02 -2.602935e-02
b  -0.08720531  1.044004e-04  9.442374e-05
u  -0.02602935  9.442374e-05  2.603217e-04

nls supports profile plots for the parameters, giving more detailed information about their uncertainty:

> plot(profile(fit))

Here is one of the three output plots showing variation in $a$:

Profile plot

E.g., a t-value of $2$ corresponds roughly to a 95% two-sided confidence interval; this plot places its endpoints around $1945$ and $1995$.

$\endgroup$
  • $\begingroup$ Oh, I almost forgot: res <- residuals(fit); res %*% res tells us that introducing the third parameter $u$ reduces the sum of squares to $2724$ (compared to $147073$ as stated in the question). $\endgroup$ – whuber Jul 23 '12 at 14:23
  • $\begingroup$ All well and good whuber. But maybe the OP had some reason to pick the exponential model (or maybe it is just because it is well known). I think first the residuals should be looked at for the exponential model. Plot them against potential covariates to see if there is structure there and not just large random noise. Before jumping into more sophisticated models try to see if a fancier model could possibly help. $\endgroup$ – Michael Chernick Jul 23 '12 at 14:25
  • 3
    $\begingroup$ Why don't you look at the original plot, Michael? It will make it abundantly obvious why at least one additional parameter is needed. Please note, too, that in a comment to the question the OP has asked, "what would be a more flexible nonlinear model?" One implication of the initial analysis offered in this answer is that it should be considered out of the ordinary to fit an exponential with fewer than three parameters: there must be some inherent constraint operating in such cases (such as an intrinsically determined unit of measurement or an intrinsic location for $x$). $\endgroup$ – whuber Jul 23 '12 at 14:32
  • 2
    $\begingroup$ I was not criticizing your answer! I did not see any residual plots. All I was suggesting is that plots of residuals vs potential covariates should be the first step in finding a better model. If I thought I had an answer to put up there I would have given an answer rather than raised my point as a constant. I thought you gave a great response and i was among the ones who gave you a +1. $\endgroup$ – Michael Chernick Jul 23 '12 at 14:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.