An oddly skewed distribution of p-values

Question

I stumbled upon an odd result which I have difficulties to explain. In the following code, $x_1$ and $x_2$ are very similar variables. Yet the distribution of p-values for the coefficient in $x_1$ is oddly skewed towards 1 while the equivalent for $x_2$ is uniform as expected.

p1=vector('double',1000)    # p-values for Model 1
p2=vector('double',1000)    # p-values for Model 2
for (i in seq(1,1000)){
  z=rnorm(100)             # A generic covariate affecting both variables
  x=0.2*rnorm(n=100)+z
  y=0.2*rnorm(n=100)+z
  x1=x-z
  x2=residuals(lm(x~1+z))
  m1=lm(y~1+x1)
  m2=lm(y~1+x2)
  p1[i]=coefficients(summary(m1))['x1',4]  
  p2[i]=coefficients(summary(m2))['x2',4]  
}

Any idea of why that is?

Answer 1

Even though x1 and x2 appear similar (their correlation is very high, as remarked by @Henry), the construction of x2 as the residual after regressing x on z imposes structure on x2 that drastically affects parameter estimation. On each iteration of your simulation, the x2 values not only sum to zero, they also are uncorrelated with the z values. As a result, in model m2 the slope estimate has mean zero but considerably smaller variance than in model m1, as you can check by comparing histograms of $\hat\beta_1$ (it's the numerator that gets much tighter). This leads to a t-statistic that has smaller variance than would be expected under the normal-theory null hypothesis, and higher p-values on average.

A simpler situation: Consider testing $\mu=0$ for a Normal($\mu=0$,$\sigma^2$) sample. Of course the null hypothesis is true, so P-values will be uniform over $[0,1]$. But suppose we alter the sample by subtracting off the sample mean from every observation. Since the population mean is zero, the sample mean is pretty small, right? The seemingly innocuous act of subtracting off the sample mean leads to a test statistic that is identically zero, and a P-value that is always 1. Loosely speaking, imposing such constraints yields a population such that normal theory thinks the null hypothesis is really, really true :).

Answer 2

You've constructed a very interesting problem. What is happening here is, essentially, that the calculations of p-values assume a certain amount of randomness in the observed correlation between, in this case, $y$ and $x_2$, under the null hypothesis, but the way you have constructed $x_2$ - as the residual from a regression of $z$ on $x$, which forces the correlation between $x_2$ and $z$ to be exactly equal to $0$ all the time - has reduced that randomness below what the null hypothesis would typically indicate. The end result is that the observed correlations between $y$ and $x_2$ are more tightly clustered around $0$ than would be expected, so the associated p-values are more tightly stacked up near $1$ than the uniform distribution that would be expected.

Mathematically, we can write $y = w + z$, where $w = $0.2*rnorm(n=100) in your code. Then, using the independence of $w$ and $z$, the null hypothesis is "expecting":

$$\rho(y,x) = f(\rho(w,x),\rho(z,x))$$

The exact equation of $f$ isn't important, although it's easy enough to write out. Due to the fact that $x$, as you have constructed it, has zero correlation with $z$, what the calculations are actually seeing is:

$$\rho(y,x) = f(\rho(w,x),0)$$

This will cause the correlation between $x$ and $y$ in the second case to almost always be smaller than the correlation between $x$ and $y$ in the first case. We can see this in the histograms below:

Note the difference in the scale on the x-axes.

Answer 3

Sorry, that I can not comment on that which I would prefer (missing rep), cause I am not so sure, if my answer is sufficient. But when actually plotting y and x1/x2 in a scatterplot you see absolutely no correlation between both. I would only expect the coefficients to be similar (and also the p-values), if there was a somehow obvious relationship between them. Given the missing relationship between both, the differing coefficients are what I would expect: garbage in garbage out.

In general, the p-values I saw when doing one iteration of your model comparisons, were above 0.1 which in many cases is the highest value that seems to be acceptable for refusing H0. So the common thing about the p-values for both x1 and x2 that we generally can take from them is, that there is no significant relationship between xi and y.

You can actually see the relationship that you want to see between xi and y when you change the factor 0.2 you apply to x and y to something greater like 1000. The coefficients and p-values are then nearly the same and you actually see the pattern (that you want to see) emerge in the plot of xi and y.