Why do we trust the p-value when fitting a regression on a single sample?

Question

I have code below that builds a linear model for a set of data:

x = rnorm(100,5,1)
b = 0.5
e = rnorm(100,0,3)
beta_0= 2.5
beta_1= 0.5
y = beta_0 + beta_1*x + e
plot(x,y)
m1 = lm(y~x)
abline(m1)
summary(m1)

When I run this block of code multiple times, the p-value can vary from 0.05 to around ~-0.7. So my question is why do we trust that a coefficient is statistically significant based only on one sample when it can vary when fitting on a different sample?

The phrase "trust the p-value" seems strange to me (what are we trusting that it should do?). You're aware that (i) the p-value is a random variable? (i) that under the null it has a uniform distribution? (iii) that under the alternative there's not some "population p-value" that you're estimating? (i.e. as sample sizes go up, it doesn't converge on some particular value, but just tends to get typically smaller, while still having some - albeit decreasing - chance of large values) — Glen_b, Sep 15 '16 at 05:16
Practically, when analyzing a data set, we report the p-value of a coefficient and indicate if it is statistically significant at some alpha level. But we are analyzing only one data set (some sample from the population). If you sample 100 subjects again, fit the regression, again and report a p-value higher than the alpha level, then the coefficient is not statistically significant. I am just confused how this works in practice (i.e. medical studies when you are analyzing only one sample typically) and you interpret this difference. — zorny, Sep 15 '16 at 05:50

score 3 · Accepted Answer · edited Apr 13 '17 at 12:44

I assume that you talk about the p-value on the estimated coefficient $\hat{\beta}_1$. (but the reasoning would be similar for $\hat{\beta}_0$).

The theory on linear regression tells us that, if the necessary conditions are fulfilled, then we know the distribution of that estimator namely, it is normal, it has mean equal to the ''true'' (but onknown) $\beta_1$ and we can estimate the variance $\sigma_{\hat{\beta}_1}$. I.e. $\hat{\beta}_1 \sim N(\beta_1, \sigma_{\hat{\beta}_1})$

If you want to ''demonstrate'' (see What follows if we fail to reject the null hypothesis? for more detail) that the true $\beta_1$ is non-zero, then you assume the opposite is true, i.e. $H_0: \beta_1=0$.

Then by the above, you know that, if $H_0$ is true that $\hat{\beta}_1 \sim N(\beta_1=0, \sigma_{\hat{\beta}_1})$.

In your regression result you observe a value for $\hat{\beta_1}$ and you can compute its p-value. If that p-value is smaller than the significance level that you decide (e.g. 5%) then you reject $H_0$ en consider $H_1$ as ''proven''.

In your case the ''true'' $\beta_1$ is $\beta_1=0.5$, so obviously $H_0$ is false, so you expect p-values to be below 0.05.

However, if you look at the theory on hyptothesis testing, then they define ''type-II'' errors, i.e. accepting $H_0$ when it is false. So in some cases you may accept $H_0$ even though it is false, so you may have p-values above 0.05 even though $H_0$ is false.

Therefore, even if in your true model $\beta_1=0.5$ it can be that you accept the $H_0: \beta_1=0$, or that you make a type-II error.

Of course you want to minimize the probability of making such type-II errors where you accept that $H_0: \beta_1=0$ holds while in reality it holds that $\beta=0.5$.

The size of the type-II error is linked to the power of your test. Minimizing the type-II error means maximising the power of the test.

You can simulate the type-II error as in the R-code below:

Note that:

if you take $\beta_1$ further from the value under $H_0$ (zero) then the type II error decreases (execute the R-code with e.g. beta_1=2) which means that the power increases.
If you put beta_1 equal to the value under $H_0$ then you find $1-\alpha$.

R-code:

x = rnorm(100,5,1)
b = 0.5
beta_0= 2.5
beta_1= 0.5

nIter<-10000
alpha<-0.05

accept.h0<-0

for ( i in 1:nIter) {
  e = rnorm(100,0,3)

  y = beta_0 + beta_1*x + e

  m1 = lm(y~x)
  p.value<-summary(m1)$coefficients["x",4]

  if ( p.value > alpha) accept.h0<- accept.h0+1 
}

cat(paste("type II error probability: ", accept.h0/nIter))

score 3 · Answer 2 · answered Sep 15 '16 at 08:16

"Trusting" the p-value may very well mean misunderstanding it. You make up a model with considerable error and sometimes the regression will detect the linear relation, some times not. The risk is determined by choosing the p-value-threshold alpha.

In the case you have proposed. Each p-value under 0.05 is "right", and each above 0.05 lacks observations. Try larger samples then n=100 and with increasing numbers you will find decreasing occurence of p-values above 0.05. So your question is essentially about the power of the test.

To find a significant correlation between x and y with a power of 90% there has to be a correlation of at least r=0.31

> library(pwr)
> pwr.r.test(n=100, sig.level = 0.05, power=0.9)

     approximate correlation power calculation (arctangh transformation) 

              n = 100
              r = 0.3164205
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

The correlation of your data is somewhere around 0.16. So the problem is not the trust in p-values but that your "study" is massively underpowered.

Find a sample of n=500 to see "wrong" p-values about one in twenty:

> pwr.r.test(r=0.16, power=.95)

     approximate correlation power calculation (arctangh transformation) 

              n = 501.0081
              r = 0.16
      sig.level = 0.05
           power = 0.95
     alternative = two.sided

Lesson learned: Never trust a not-significant p-value without a sound power analysis.

Why do we trust the p-value when fitting a regression on a single sample?

2 Answers2