Linear model with multiple measures

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


• @Jonas Berge: runif(1000) gives 1000 random draws from a uniform distribution, rep(runif(10), each=100) only 10 draws that are repeated 100 times, so for rep(runif(10), each=100) I have only ten different values. The chance that I have only 10 different values runif(1000) are rather small I think ?
– user83346
Commented Sep 9, 2015 at 7:50
• @Elvis: well I have used it long time ago, I think it is the other way around, mixed effect models are a special kind of var-covar models namely (if my memories are well) the mixed effect model estimates a compound symmetry var-covar. Anyhow, if you have repeated measures (for the same individual) then you introduce dependencies among the observations and your regression assumptions are violated. see also stats.stackexchange.com/questions/166434/…
– user83346
Commented Sep 9, 2015 at 8:52
• I tried to add subject id as a random effect using a linear mixed model to account for the repeated measures. But the p-values for u were still lower than expected when repeating the procedure as you did. I think that the problem is that there are so few subjects but I hope that somebody else (f coppens?) can explain this.
– JonB
Commented Sep 9, 2015 at 9:38
• It is an interesting example that may have practical implications in similar situations. It should be the same situation if you have 10 schools with 100 students in each school, and you include school-level predictors as well as school id as a random effect. This might increase the chance of type 1 errors for the school level predictors?
– JonB
Commented Sep 9, 2015 at 9:38
• I deleted my answer and added an entirely new one. Try the code in my answer now, I think the mystery is solved. :)
– JonB
Commented Sep 9, 2015 at 15:06

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. • (+1) you should add a simulation where the standard deviation of the rnorm that is added to het$y\$ is e.g. sd=5
– user83346
Commented Sep 9, 2015 at 10:59
• Good idea, I'll try it out right away.
– JonB
Commented Sep 9, 2015 at 11:03
• Yeah, I too thought that including the random individual effect would account for x. I think that the answer has to do with what you said earlier, something about there always being a correlation. If we do a<-rnorm(10); b<-rnorm(10); cor.test(a,b) there will be an insignificant correlation that is != 0. If we then: a <- rep(a, 100); b <- rep(b, 100); cor.test(a,b) the correlation will be significant. I think we have a similar situation here, with a correlation made highly significant by it being repeated in each repeated measure, and the mixed model is unable to ..
– JonB
Commented Sep 9, 2015 at 12:18
• .. unable to separate the individual random effect from this artificial correlation. I hope that somebody can provide a more technical explanation.
– JonB
Commented Sep 9, 2015 at 12:19
• @Elvis: I edited my answer, see at the bottom, I have to thni about it but maybe that can help you further ?
– user83346
Commented Sep 9, 2015 at 13:38