I have the slope and intercept of my model as well as their respective standard error but not the data from which it was estimated (I retrieved it from the literature). I would like to plot the predictions a*x+b with the prediction interval around the curve, to know if my own data falls into it. Do I have enough information to do this, and if yes, how (in R or python)?


migrated from stackoverflow.com Jan 13 '17 at 12:09

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ Is this something you want? tomholderness.wordpress.com/2013/01/10/confidence_intervals $\endgroup$ – periphreal Jan 13 '17 at 11:25
  • 3
    $\begingroup$ No, you don't have enough information. You are at least missing the covariance between the parameter estimates. Also, your data is not expected to fall into the confidence band. Maybe you are confusing that with the prediction band. $\endgroup$ – Roland Jan 13 '17 at 11:32
  • 1
    $\begingroup$ If you want to see if your data lie within an interval, it's not a confidence interval you seek. A confidence interval is an interval for the population mean, not a sample. If the data is not from the same sample as was used to generate the fit you want a prediction interval. $\endgroup$ – Glen_b Jan 13 '17 at 12:14
  • $\begingroup$ Note also that requests for code are generally regarded as off-topic here (not that you're responsible for it ending up on crossvalidated). The remainder of your question is on topic however. $\endgroup$ – Glen_b Jan 13 '17 at 12:17
  • 1
    $\begingroup$ @alpagarou Yes, you could use the bootstrap method, but you'd need the original data for bootstrapping. You could also assume two independent normal distributions for intercept and slope using the estimates and standard errors, but such an assumption of independence is not justified. And you need the covariance to simulate from a multivariate normal distribution. $\endgroup$ – Roland Jan 13 '17 at 14:23

As @Roland points out, you also need the covariance of the parameters. If you have those, for example, you have a used curve_fit in python...

from pylab import *
from scipy.optimize import curve_fit

# fake data
x = arange(0, 30)
y = x + 5*np.random.normal(size=len(x))

# fit
f = lambda x, *p: polyval(p, x)
p, cov = curve_fit(f, x, y, [1, 1])

... you could assume that the fit (p, cov) represents a normal distribution and sample from it, selecting the lower and upper 2.5% quantiles:

# simulated draws from the probability density function of the regression
xi = linspace(np.min(x), np.max(x), 100)
ps = np.random.multivariate_normal(p, cov, 10000)
ysample = np.asarray([f(xi, *pi) for pi in ps])
lower = percentile(ysample, 2.5, axis=0)
upper = percentile(ysample, 97.5, axis=0)

# regression estimate line
y_fit = poly1d(p)(xi)

# plot
plot(x, y, 'bo')
plot(xi, y_fit, 'r-')
plot(xi, lower, 'b--')
plot(xi, upper, 'b--')

Again as @Roland points out, the data doesn't, in general, fall between the confidence intervals of the fit. Also, to see that the covariances are important, you can diagonalize the covariance matrix with cov = np.diag(np.diag(cov)) and re-draw.

There's a more analytical solution specific to the linear fit here (with very slightly different results for some mathy reason).


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.