I am struggling to calculate / define confidence bands to use with a custom function that I must fit to my data. I have seen plenty of examples for standard models (linear regression, log functions, etc.). However, calculating the intervals to add bands to my custom function is proving tricky. Below I have the data and my function (working in R):

data <- matrix(c(0.08, 0.1, 0.12, 0.13, 0.49, 0.11, 0.12, 0.15, 
                 0.22, 0.47, 7, 8, 9, 21, 30, 3, 8, 13, 15, 17 ), ncol=2)
mycurve <- function(x){ a + (b*log(x)) }
a    <- 31
b    <-  9
data <- as.data.frame(data) # So that ggplot can use it
ggplot(data, aes(V1,V2)) + geom_point() + stat_function(fun=mycurve, color="red")

And the result is:

enter image description here

The red curve is the custom function that I must use for this dataset. This is one example, but I have a few more datasets and must fit a different function to to each dataset. So I'm looking for an approach for calculating confidence bands that I can use for any function.

I have spent quite a while on this, but haven't yet found a "generic" solution.

  • 1
    $\begingroup$ 1. I don't see any fitting going on there; this makes it a bit hard to guess what the situation actually is -- how are the coefficients obtained? 2. Are you seeking a pointwise interval or is the interval a simultaneous one for the whole curve in the range of the data? $\endgroup$ – Glen_b -Reinstate Monica Jul 22 '15 at 17:25
  • $\begingroup$ @Glen_b Please see my explanation in response to Dalton's comment below. I'm looking to draw the interval for the whole curve i.e. confidence bands rather than just intervals around each point. $\endgroup$ – Dirk Snyman Jul 22 '15 at 17:42
  • $\begingroup$ You misunderstand -- you will get bands either way. The distinction refers to how the proportion falling in the bands is computed. "Pointwise" means that at each point there's a CI calculation (this is the more common one); the CIs taken together give a smooth curve. . $\endgroup$ – Glen_b -Reinstate Monica Jul 22 '15 at 17:56
  • $\begingroup$ @Glen_b are "pointwise" more common? When I think confidence band, I think a continuous curve within the range of the observed $x$. And by each point, do you mean an estimate of the mean $y$ at each observed $x$ value? I assume you do, but it might not be clear to other readers. $\endgroup$ – Dalton Hance Jul 22 '15 at 18:11
  • $\begingroup$ @DaltonHance See here, which has this diagram of the distinction in regression. As you see, both are "continuous bands", but the pointwise interval is narrower. $\endgroup$ – Glen_b -Reinstate Monica Jul 22 '15 at 22:54

As Glen_b pointed out, you aren't actually fitting a relationship here. You're simply plotting a curve. There is no obvious relationship between this curve and the data.

However, if what you want to do is fit a curve that follows the form you've given: $y = a+b*\mathrm{log}(x)$, where $x$ and $y$ are data and you want to find $a$ and $b$ then I can answer your question. A generic solution to the problem of fitting a non-linear relationship between some predictor, $x$, and some predicted value $y$ is to use non-linear least squares. This carries some of the usual assumptions of linear regression, that is that the $x$ values are measured without error and that the variance of the residuals of the relationship are i.i.d. The computation of the relationship is no where near as easy as for simple linear regression, but luckily we have computers for that.

This is implemented in R by the nls function. You will need to provide starting values for the coeffiecients. Typically this is done by linearizing the relationship (that is finding a transformation for $y$ such that the relationship on the RHS is linear, then fitting a linear regression and using those estimates as starting values.

To get confidence bands on the relationship, you will need to use an error propagating function. An implementation of that for nls in R is available in the propagate package in the function predictNLS.

EDIT: Working example and a note that the fit curve you gave does not match the point estimate of the best fit.


DirkSnyman <- data.frame(x = c(0.08, 0.1, 0.12, 0.13, 0.49, 
                 0.11, 0.12, 0.15, 0.22, 0.47), 
            y = c(7, 8, 9, 21, 30, 
                  3, 8, 13, 15, 17))

DS.fit <- nls(y ~ a+b*log(x), data = DirkSnyman, start = list(a = 31, b = 9))

preds <- data.frame(x = seq(0.05, 0.5, 0.01))
prop <- predictNLS(DS.fit, newdata = preds)

preds$mean <- prop$summary[,2]
preds$lcl <- prop$summary[,5]
preds$ucl <- prop$summary[,6]

mycurve <- function(x){ a + (b*log(x)) }
a    <- 31
b    <-  9

ggplot(DirkSnyman, aes(x,y)) + geom_point() + 
  stat_function(fun=mycurve, color="red") + 
  geom_line(data = preds, aes(x = x, y = mean), color= "blue") +
  geom_ribbon(data = preds, 
              aes(x = x, y = mean, ymin = lcl, ymax = ucl), 
              color= "blue", alpha = .5)

enter image description here

  • $\begingroup$ I fitted the original curve in TableCurve software. So a and b are the estimated parameters from that process. TableCurve gave the equation as y=a+blnx with a=31.02 (s.e.= 5.4) and b=9.9 (s.e.=2.8). So I used the equation and parameter values for drawing my graphs in R as posted above. $\endgroup$ – Dirk Snyman Jul 22 '15 at 17:41
  • $\begingroup$ Why not just do the whole process in R? nls gives the same answer, and you can do everything in one environment and achieve your confidence bands. $\endgroup$ – Dalton Hance Jul 22 '15 at 18:06
  • $\begingroup$ your code was great for log and exp functions, but I struggled with functions like y ~ a*sinpi(2*x/c+b) -- I added an edit to your answer with the new data, function and associated hassles. $\endgroup$ – Dirk Snyman Jul 24 '15 at 6:41

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.