Consider the following situation: I have several individuals (say 10), and for each of these individuals there is a covariate $x$ which linearly determines an outcome $y$. I can measure $y$ with some normal error, multiple times per individual (say 100 times). This code should clarify what I mean:
> D <- data.frame( id = rep(1:10, each=100) )
> set.seed(1)
> D$x <- rep(runif(10), each = 100)
> D$y <- 1 + 0.5*D$x + rnorm(1000)
So I have $Y = \alpha + \beta X + \epsilon$, a classical linear model. I can estimate the coefficients, make Wald tests, etc.
> summary( lm(D$y ~ D$x) )
Call:
lm(formula = D$y ~ D$x)
Residuals:
Min 1Q Median 3Q Max
-3.0409 -0.6810 -0.0257 0.6998 3.7819
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.05088 0.06857 15.326 < 2e-16 ***
D$x 0.38846 0.10926 3.555 0.000395 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.035 on 998 degrees of freedom
Multiple R-squared: 0.01251, Adjusted R-squared: 0.01152
F-statistic: 12.64 on 1 and 998 DF, p-value: 0.0003953
So far, so good?
Is that correct to keep these 100 measures per individuals? The reason we take multiple measures is that there is an important measure error; it is natural to summarize them by their mean, leading to a more precise measure. This can be done:
> library(plyr)
> D2 <- ddply(D, c("id","x"), summarize, z = mean(y))
> D2
id x z
1 1 0.26550866 1.215321
2 2 0.37212390 1.178798
3 3 0.57285336 1.336295
4 4 0.90820779 1.490424
5 5 0.20168193 1.042107
6 6 0.89838968 1.429626
7 7 0.94467527 1.252171
8 8 0.66079779 1.304997
9 9 0.62911404 1.334422
10 10 0.06178627 1.067025
Now we go back to our linear model:
> summary( lm(D2$z ~ D2$x) )
Call:
lm(formula = D2$z ~ D2$x)
Residuals:
Min 1Q Median 3Q Max
-0.16567 -0.01444 0.01359 0.05577 0.08674
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.05088 0.05396 19.475 5.02e-08 ***
D2$x 0.38846 0.08598 4.518 0.00196 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.08142 on 8 degrees of freedom
Multiple R-squared: 0.7184, Adjusted R-squared: 0.6832
F-statistic: 20.41 on 1 and 8 DF, p-value: 0.001956
Same parameters estimates (which is logic), but the standard error estimate is different, and so are the degrees of freedom... leading to a different $p$-value.
My first thought is that the first analysis (keep all the measures) is the better solution: by taking the mean we lose some information on the residual variance (?).
But... imagine that I don’t observe $x$, and I have some other variable $u$ which has nothing to do with $y$ — and I want to test it. Something like this:
> D$u <- rep(runif(10), each = 100)
> summary( lm(D$y ~ D$u) )
Call:
lm(formula = D$y ~ D$u)
Residuals:
Min 1Q Median 3Q Max
-3.2210 -0.6985 -0.0245 0.7045 3.7181
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1872 0.0524 22.658 <2e-16 ***
D$u 0.2074 0.1087 1.908 0.0567 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.039 on 998 degrees of freedom
Multiple R-squared: 0.003635, Adjusted R-squared: 0.002637
F-statistic: 3.641 on 1 and 998 DF, p-value: 0.05666
This looks good, but let’s check further: $u$ and $y$ are independent, so we hope that the $p$-value is uniformly distributed:
p <- replicate( 1e3, {D$u <- rep(runif(10), each = 100);
summary(lm(D$y ~ D$u))$coefficients[2,4]})
plot(seq(0,1,length=1000), sort(p), pch=".",
xlab="theoretical quantiles", ylab="empirical quantiles")
abline(0,1,col=3)
...this doesn’t work. There are too many small $p$-values!! If we try with the means of the outcomes (the data frame D2) it is ok.
I understand that both models $y = \alpha + \beta u + \epsilon$ (1) and $y = \alpha + \epsilon$ (2) are false, and that I am implicitly comparing these two models, but I don’t see why it should be so wrong. Note that the normal qq-plot of $y$ looks as linear as you may hope, and you wouldn’t discard model (2) simply by eyeballing.
So, I you read me so far... what are your thoughts? Is that phenomenon well known to more experienced statisticians? Does it have a name? In a practical situation with multiple measures, what do you recommend? Many thanks in advance for any comment.
PS. It can be more natural to check the test behaviour by permutation of the values of $x$, while respecting the group structure. This can be done with D$u <- rep( sample(unique(D$x)), each = 100 ) and leads to a very similar qq-plot.
PPS. After rethinking it, I understand that when $x$ is unobserved, this creates a strong "individual effect" which acts as confounder when testing the effect of $u$. When the multiple measures are condensed on their mean, this individual effect is "absorbed" in the residual variance. When the multiple measures are kept, this won’t work anymore, the residual variance estimation is dominated by the intra-group variance — I’d have to think about it but I guess I’ll finally understand what’s going on here. I am still interested by other explanations, references to this phenomenon in the literature, recommendations, and any kind of comments.
Update after discussion with f coppens (see below his answer) the main problem (the only problem?) is the presence of a structure in the residuals due to the effect of the hidden variable $x$. A possible solution could be to introduce a random individual effect. Any further comments are welcome.
