An oddly skewed distribution of p-values

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?

• x1 is the difference of two std normal distributions, so a std.n.d. again, x2 is the residuals of a linear model with an intercept, are you sure these two variables are similar? – user2974951 Sep 13 '18 at 8:02
• @user2974951 - they are very similar except that x1 has a sample average close to $0$ while x2 has a sample average exactly $0$. The sample correlation between them typically seems to be over $0.99$ – Henry Sep 13 '18 at 8:12
• To see the similarity visually try plot(x1,x2) – Henry Sep 13 '18 at 8:16
• I am seeing plot.ecdf(p1) as close to uniform but plot.ecdf(p2) skewed towards $1$, the opposite of your question – Henry Sep 13 '18 at 8:26
• @Henry Indeed this is what I meant. – user1764944 Sep 13 '18 at 20:18

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

• ... hadn't noticed your answer (+1) appearing while I was working on mine! – jbowman Sep 14 '18 at 16:47

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.

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.

• Well, xi is indeed not correlated with y (which is intentional) and so I would expect the null hypothesis to be true and the distribution of p-values to be uniform. x1 and x2 are highly correlated between each other (the mean correlation between them is above 0.99) and yet less than 0.5% of p-values for x2 are below 0.5 while about 50% of p-values for x1 are below 0.5. – user1764944 Sep 14 '18 at 13:44
• True, i actually have problems to explain why those differences appear in detail, but i guess the strongly differing coefficients you get during each run explain the difference of the p-values. I also do not really get why the coefficients are so different in the first place, given the data is so similar. But somehow this is not the case when the data provides actual patterns to learn. And since this is the case, I am generally interested in, I do not worry to much. Sorry, that I can not give you a final answer to your question, but maybe a direction to investigate. GL! – Elmar Macek Sep 14 '18 at 14:51