Linear model with multiple measures

Question

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.

Answer 1

By the way you construct the $y$ and the $u$ I think they are dependent; for each $id$ in $D$ you have the same value for $x$ and the same valuefor $u$, as $y$ is a linear function (+ a random error) of $x$ this 'creates' in my opinion a dependency between $u$ and $y$.

If you make the $u$ independent of $id$ then the p-vlaues are uniform, try this instruction ( a different u for every row of D):

p <- replicate( 1e3, {D$u <- rep(runif(1000), each = 1); 
                      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)

EDIT: After the comments I added this:

I have to think about it, but because of the same values for $x$ I have the feeling thet the variance of the slope coefficient is underestimated. I also have the impression that with a mixed effect model the coefficient becomes insignificant:

library(nlme)

D <- data.frame( id = rep(1:10, each=100) )
set.seed(1)
D$x <- rep(runif(10), each = 100)
D$y <- 1 + 5*D$x + rnorm(1000)

D$u <- rep(runif(10), each = 100)
summary( lm(D$y ~ D$u) )




lme.1<-lme(y ~ u + 1,
    random= ~ 1|id,
    data=D,
    control=lmeControl(opt="optim"),
    weights=varIdent(form=~1|id),
    method="REML")
summary(lme.1)

EDIT 2: It takes some time to execute but I think this is what we expect:

t.I <- replicate( 1000, {D$u <- rep(runif(10), each = 100); 
                         lme.1<-lme(y ~ u + 1,
                                    random= ~ 1|id,
                                    data=D,
                                    control=lmeControl(opt="optim"),
                                    weights=varIdent(form=~1|id),
                                    method="REML");
                         anova(lme.1)[["p-value"]][2]
                         }
                  )


plot(seq(0,1,length=1000), sort(t.I), pch=".", 
     xlab="theoretical quantiles", ylab="empirical quantiles", main="lme corrected for repeated measures")
abline(0,1,col=3)

Answer 2

EDIT 4:

I deleted my previous answer as I had made a mistake that explained the unexpected results, so no need for anyone to waste time working through all of that. Luckily, with the help from f coppens, I realized the mistake and I can now post what I believe will be a final and short answer.

The short answer is that you need to take the repeated measurements into account. This can be done using a mixed model:

library(lme4)
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)
D$u <- rep(runif(10), each = 100)
summary(lmer (y ~ u + (1|id), data=D, REML=F))

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.004164 0.06453 
 Residual             1.073876 1.03628 
Number of obs: 1000, groups:  id, 10

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  1.18725    0.06155 10.00000  19.290 3.05e-09 ***
u            0.20743    0.12769 10.00000   1.625    0.135    

So u is not significant at the 0.05 level anymore.

To address the second question, regarding the theoretical vs empirical p-values, we can now test it with the mixed model using the following (correct) code:

n1 <- 10 # number of individuals
n2 <- 10 # number of repeated measures
rep <- 1000 # number of iterations
D <- data.frame( id = rep(1:n1, each=n2) )
set.seed(1)
p <- foreach (i = 1:(rep), .combine=c) %dopar% {
  D$x <- rep(runif(n1), each = n2)
  D$y <- 1 + 0.5*D$x + rnorm(n1*n2)
  D$u <- rep(runif(n1), each = n2)
  model <- lmer (y ~ u + (1|id), data=D)
  summary(model)$coefficients[10]
}
plot(seq(0,1,length=rep), sort(p), pch=".", 
 xlab="theoretical quantiles", ylab="empirical quantiles")
abline(0,1,col=3)

prop.table(table(ifelse(p<=0.05, 1, 0))) # proportion of type 1 errors at the 0.05 level.
    0     1 
0.966 0.034 

And it is clear that this solved the problem. No need for different variances per subject, as far as I can tell.