18
$\begingroup$

We know that a paired t-test is just a special case of one-way repeated-measures (or within-subject) ANOVA as well as linear mixed-effect model, which can be demonstrated with lme() function the nlme package in R as shown below.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

When I run the following paired t-test:

t.test(x1, x2, paired = TRUE)

I got this result (you will get a different result because of the random generator):

t = -2.3056, df = 9, p-value = 0.04657

With the ANOVA approach we can get the same result:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Now I can obtain the same result in lme with the following model, assuming a positive-definite symmetrical correlation matrix for the two conditions:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

Or another model, assuming a compound symmetry for the correlation matrix of the two conditions:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

With the paired t-test and one-way repeated-measures ANOVA, I can write down the traditional cell mean model as

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

where i indexes condition, j indexes subject, Yij is the response variable, μ is constant for the fixed effect for overall mean, αi is the fixed effect for condition, βj is the random effect for subject following N(0, σp2) (σp2 is population variance), and εij is residual following N(0, σ2) (σ2 is within-subject variance).

I thought that the cell mean model above would not be appropriate for the lme models, but the trouble is that I can't come up with a reasonable model for the two lme() approaches with the correlation structure assumption. The reason is that the lme model seems to have more parameters for the random components than the cell mean model above offers. At least the lme model provides exactly the same F-value, degrees of freedom, and p-value as well, which gls cannot. More specifically gls gives incorrect DFs due to the fact that it does not account for the fact that each subject has two observations, leading to much inflated DFs. The lme model most likely is overparameterized in specifying the random effects, but I don't know what the model is and what the parameters are. So the issue is still unresolved for me.

$\endgroup$
  • 2
    $\begingroup$ Not quite sure what you're asking. The model you wrote down is precisely the model for the random effects model; the correlation structure is induced by the random effect. $\endgroup$ – Aaron Feb 22 '12 at 4:29
  • $\begingroup$ @Aaron: the random-effect βj in the cell mean model is supposed to follow N(0, σp2). My confusion is, how is this term (with only one parameter σp2) associated with the correlation structure specified by either compound symmetry or a simple symmetrical matrix in the lme model? $\endgroup$ – bluepole Feb 22 '12 at 4:41
  • $\begingroup$ When you calculate the correlation between the two observations on the same subject, the correlation is sigma_p^2/(sigma_p^2+sigma^2) because they share the same beta_j. See Pinheiro/Bates p.8. Also, the random effect model as you wrote it is equivalent to compound symmetry; other correlation structures are more complex. $\endgroup$ – Aaron Feb 22 '12 at 15:33
  • $\begingroup$ @Aaron: Thanks! I've already read Pinheiro/Bates book about this, and still couldn't figure out the specifics about the random effects. The more relevant pages seem to be the example at P.160-161. Also, the random-effects output from lme() with compound symmetry assumption does not seem to agree with the correlation of σp2 / (σp2 + σ2) in the cell mean model. Still baffled about the model structure. $\endgroup$ – bluepole Feb 22 '12 at 16:07
  • $\begingroup$ Well, almost equivalent to compound symmetry; in CS the correlation can be negative but not with random effects. Perhaps that's where your difference arises. See stats.stackexchange.com/a/14185/3601 for details. $\endgroup$ – Aaron Feb 22 '12 at 16:11
14
$\begingroup$

The equivalence of the models can be observed by calculating the correlation between two observations from the same individual, as follows:

As in your notation, let $Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}$, where $\beta_j \sim N(0, \sigma_p^2)$ and $\epsilon_{ij} \sim N(0, \sigma^2)$. Then $Cov(y_{ik}, y_{jk}) = Cov(\mu + \alpha_i + \beta_k + \epsilon_{ik}, \mu + \alpha_j + \beta_k + \epsilon_{jk}) = Cov(\beta_k, \beta_k) = \sigma_p^2$, because all other terms are independent or fixed, and $Var(y_{ik}) = Var(y_{jk}) = \sigma_p^2 + \sigma^2$, so the correlation is $\sigma_p^2/(\sigma_p^2 + \sigma^2)$.

Note that the models however are not quite equivalent as the random effect model forces the correlation to be positive. The CS model and the t-test/anova model do not.

EDIT: There are two other differences as well. First, the CS and random effect models assume normality for the random effect, but the t-test/anova model does not. Secondly, the CS and random effect models are fit using maximum likelihood, while the anova is fit using mean squares; when everything is balanced they will agree, but not necessarily in more complex situations. Finally, I'd be wary of using F/df/p values from the various fits as measures of how much the models agree; see Doug Bates's famous screed on df's for more details. (END EDIT)

The problem with your R code is that you're not specifying the correlation structure properly. You need to use gls with the corCompSymm correlation structure.

Generate data so that there is a subject effect:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Then here's how you'd fit the random effects and the compound symmetry models.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

The standard errors from the random effects model are:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

And the correlation and residual variance from the CS model is:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

And they're equal to what is expected:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Other correlation structures are usually not fit with random effects but simply by specifying the desired structure; one common exception is the AR(1) + random effect model, which has a random effect and AR(1) correlation between observations on the same random effect.

EDIT2: When I fit the three options, I get exactly the same results except that gls doesn't try to guess the df for the term of interest.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(The intercept is different here because with the default coding, it's not the mean of all subjects but instead the mean of the first subject.)

It's also of interest to note that the newer lme4 package gives the same results but doesn't even try to compute a p-value.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283
$\endgroup$
  • $\begingroup$ Thanks again for the help! I know this part from the perspective of cell mean model. However, with the following result from lme() with compound symmetry: Random effects: Formula: ~x - 1 | subj Structure: Compound Symmetry StdDev xx1 1.1913363 xx2 1.1913363 Corr: -0.036 Residual 0.4466733. I still can't reconcile these numbers with the cell mean model. Maybe you can further help me sort these numbers out? $\endgroup$ – bluepole Feb 22 '12 at 16:50
  • $\begingroup$ Also, any thoughts about the model formulation with other correlation structures such as simple symmetrical matrix? $\endgroup$ – bluepole Feb 22 '12 at 16:58
  • $\begingroup$ I see! I should have read your response in the other thread more carefully. I thought about using gls() before, but failed to figure out the correlation specification. It's interesting that lme() with compound symmetry structure for the random effect still renders the same t-value, but it seems the variances for the random effects are not directly interpretable. I really appreciate your help! $\endgroup$ – bluepole Feb 22 '12 at 18:28
  • $\begingroup$ After some 2nd thought, I feel that my original confusion is still unresolved. Yes, gls can be used to demonstrate the correlation structure and mean squared rums, but the model underneath it is not exactly the same as the paired-t test (or one-way repeated-measures ANOVA in general), and such an assessment is further supported by the incorrect DFs and p-value from gls. In the contrast, my lme command with compound symmetry provides the same F, DFs, and p-value. The only thing I'm puzzled about is how the lme model is parameterized as stated in my original post. Any help out there? $\endgroup$ – bluepole Feb 25 '12 at 22:27
  • $\begingroup$ Not sure how to help you. Could you write out what you think the two different models are? Something is wrong in how you're thinking about one of them. $\endgroup$ – Aaron Feb 26 '12 at 13:04
3
$\begingroup$

You might also consider using function mixed in package afex to return p values with Kenward-Roger df approximation, which returns identical p values as a paired t test:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

Or

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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