1
$\begingroup$

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)

qq plot

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

$\endgroup$
1
  • 1
    $\begingroup$ Are you sure that $u$ and $y$ are independent ? $\endgroup$
    – user83346
    Sep 9, 2015 at 7:24

2 Answers 2

1
$\begingroup$

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)

enter image description here

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)

enter image description here

$\endgroup$
30
  • $\begingroup$ @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 ? $\endgroup$
    – user83346
    Sep 9, 2015 at 7:50
  • 1
    $\begingroup$ @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/… $\endgroup$
    – user83346
    Sep 9, 2015 at 8:52
  • 1
    $\begingroup$ 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. $\endgroup$
    – JonB
    Sep 9, 2015 at 9:38
  • 1
    $\begingroup$ 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? $\endgroup$
    – JonB
    Sep 9, 2015 at 9:38
  • 1
    $\begingroup$ I deleted my answer and added an entirely new one. Try the code in my answer now, I think the mystery is solved. :) $\endgroup$
    – JonB
    Sep 9, 2015 at 15:06
1
$\begingroup$

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 

enter image description here

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

$\endgroup$
12
  • $\begingroup$ (+1) you should add a simulation where the standard deviation of the rnorm that is added to het $y$ is e.g. sd=5 $\endgroup$
    – user83346
    Sep 9, 2015 at 10:59
  • $\begingroup$ Good idea, I'll try it out right away. $\endgroup$
    – JonB
    Sep 9, 2015 at 11:03
  • 1
    $\begingroup$ 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 .. $\endgroup$
    – JonB
    Sep 9, 2015 at 12:18
  • 1
    $\begingroup$ .. unable to separate the individual random effect from this artificial correlation. I hope that somebody can provide a more technical explanation. $\endgroup$
    – JonB
    Sep 9, 2015 at 12:19
  • 1
    $\begingroup$ @Elvis: I edited my answer, see at the bottom, I have to thni about it but maybe that can help you further ? $\endgroup$
    – user83346
    Sep 9, 2015 at 13:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.