How can I incoporate error bars into my P values for linear regression in python?

I'm interested in statistically validating linear regression problems in python. Traditionally, these problems can be solved with scipy's linregress function. For example:

x = np.linspace(0,1,25)
y = 0.5*x + np.random.normal(0,0.15,len(x))
err = np.random.uniform(3.8,0.5,len(x))
plt.scatter(x,y)

then we can use linregress(x,y) to compute our p value. In this case we obtain a pvalue=1.3e-8 so our fit is significant, which seems reasonable given our plot.

However, the picture changes if we also plot the error bars:

Now, given the size of the error, the conclusion that the fit is significant seems suspect. Is there a way to incorporate information about the size of the errors into a pvalue test in python?

As far as i know, the common linear regression just minimizes the squared sum of the errors to the regression line, so it doesn't take into account the individual errors of the data points.

What I think is that you may have a interpretation error of the p-value, even if the error is absolutely huge as it is the case, the correlation and slope looks to be there.

Think it like this, if the error bars are sooo huge, isn't it weird that you have such a well defined ascending line by the points? so that's why the p-value is small.

From the docs:

> pvalue float
>
> The p-value for a hypothesis test whose null hypothesis is that the
> slope is zero, using Wald Test with t-distribution of the test
> statistic.

So for me it looks ok, also you can think of the case where your measurement is precise but not accurate, so you may have extremely big shifts in the y axis (hence errors bars if you want) like in a non calibrated instrument (with a linear response), that would still not a affect to that p-value, and it is kind of similar case to this one.