1
$\begingroup$

I've used scipy's linregress function to fit lines through the points in the four panels below. The function outputs a p-value, shown on each panel. The p-values are large for panels a and b (I can't reject the null hypothesis that the slope is zero) and small for panels c and d (I can reject the null hypothesis.

I am not happy with this. I want panels a and c to pass, and panels b and d to fail the test (i.e. the spread of the points relative to the line should determine whether the test passes or fails, irrespective of what the slope is). I am looking for a function that will output a p-value, which I can then use to determine whether the test passes or fails (depending on whether it's below or greater than say 0.05 or 0.01).

Find the code below the figure...

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress

plt.figure(figsize=(7, 7))
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0))
ax4 = plt.subplot2grid((2, 2), (1, 1))

# axes 1                                                                                                                                                      
x = np.random.randn(40)
y = np.random.randn(40)*0.1
ax1.scatter(x, y, c='k')

slope, intercept, r_value, p_value, std_err = linregress(x, y)
print p_value
ax1.plot(x, intercept + slope*x, 'r')
ax1.text(0.65, 0.1, 'p = ' +  str(np.round(p_value, 3)), transform=ax1.transAxes)
ax1.set_title('(a)')

# axes 2                                                                                                                                                      
x = np.random.randn(40)
y = np.random.randn(40)*0.8
ax2.scatter(x, y, c='k')
slope, intercept, r_value, p_value, std_err = linregress(x, y)
print p_value
ax2.plot(x, intercept + slope*x, 'r')
ax2.text(0.65, 0.1, 'p = ' +  str(np.round(p_value, 3)), transform=ax2.transAxes)
ax2.set_title('(b)')

# axes 3                                                                                                                                                      
x = np.random.randn(40)
y = x + np.random.randn(40)*0.1
ax3.scatter(x, y, c='k')

slope, intercept, r_value, p_value, std_err = linregress(x, y)
print p_value
ax3.plot(x, intercept + slope*x, 'r')
ax3.text(0.65, 0.1, 'p = ' +  str(np.round(p_value, 3)), transform=ax3.transAxes)
ax3.set_title('(c)')

# axes 4                                                                                                                                                      
x = np.random.randn(40)
y = x + np.random.randn(40)*0.8
ax4.scatter(x, y, c='k')

slope, intercept, r_value, p_value, std_err = linregress(x, y)
print p_value
ax4.plot(x, intercept + slope*x, 'r')
ax4.text(0.65, 0.1, 'p = ' +  str(np.round(p_value, 3)), transform=ax4.transAxes)
ax4.set_title('(d)')

ax1.axis([-2, 2, -2, 2])
ax2.axis([-2, 2, -2, 2])
ax3.axis([-2, 2, -2, 2])
ax4.axis([-2, 2, -2, 2])
plt.tight_layout()
plt.ion(); plt.show() 
$\endgroup$
  • $\begingroup$ The p-values presented in the graphs come from the test of the null hypothesis that slop = 0. What you wanted is to test the variance of error terms. $Y=X\beta + \epsilon$ Your null hypothesis may be like $\mathrm{Var}(\epsilon) = a$ and alternative hypothesis is $\mathrm{Var}(\epsilon) > a$. But you need to specify $a$ first. $\endgroup$ – user158565 Jun 3 '17 at 2:45
2
$\begingroup$

A test on the slope is meant to check if the slope is significantly different from zero (i.e. is there a slope or not). If what you want to test is whether deviations from the line of best fit are significant or not, your random variable of interest would be the regression residuals

$$\epsilon = y - X\beta$$

Since $\epsilon$ isn't observed, you can use $\hat{\epsilon}$. Let's suppose the OLS assumption $\mathbb{E}(\epsilon) = 0$ is met and furthermore, we have homoscedasticity of errors

$$\mathbf{Var}(\epsilon\mid X) = \mathbf{Var}(\epsilon) = \sigma^2$$

Then asymptotically

$$\hat{\sigma}^2 := \frac{1}{n-p}\sum_i \hat{\epsilon}_i^2 \sim \chi^2_{n-p}$$

where $p$ is the number of parameters fit in the linear regression (eg. $p = 2$ if you fit a slope and intercept term to get $\hat{\epsilon}$).

It sounds from your question that you'd want to test the hypothesis that $0 \leq \sigma^2 < \delta$ where $\delta$ is some small number. In the case where you choose $\delta > 0$ for your test, you can follow the standard routine for obtaining a p-value.

However, you can see from this image a la Wikipedia that the $\chi^2$ distribution has no density at exactly zero (actually this is true for all distributions with continuous support):

Chi-Squared PDF

So, in the case you want to test $\delta = 0$, check out @whuber's comment directly following this post.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks. I guess I need to use the residuals from Orthogonal Least Squares not Ordinary Least Squares? Otherwise the error will be more forgiving when there is little or no slope. $\endgroup$ – Oliver Angelil Jun 3 '17 at 20:57
  • $\begingroup$ Whether you use Orthogonal Least Squares (which has its own set of assumptions) or OLS, the procedure is about the same: Residuals will be left behind and a $\chi^2$ test on variance can be used. The difference between the two methods is how residuals are measured and, certainly, it can be shown that the residuals are non-zero in the OLS case if an only if they are non-zero in the Orthogonal Least Squares case. $\endgroup$ – Mustafa S Eisa Jun 5 '17 at 3:11

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.