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)?
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).