Paired t-test vs general mixed model with R

Dataset

Consider the following dataset:

> x1 <- c(3,4,3,5,6,6,6,7,7,6,8,9,11)
> x2 <- c(3,4,5,5,8,6,8,7,9,8,8,9,14)
> cbind(x1,x2)
      x1 x2
 [1,]  3  3
 [2,]  4  4
 [3,]  3  5
 [4,]  5  5
 [5,]  6  8
 [6,]  6  6
 [7,]  6  8
 [8,]  7  7
 [9,]  7  9
[10,]  6  8
[11,]  8  8
[12,]  9  9
[13,] 11 14

where x1 and x2 are paired sample.

Paired t-test

The process of a typical paired t-test with R is following:

> t.test(x1, x2, paired=T)

 Paired t-test

data:  x1 and x2
t = -3.1225, df = 12, p-value = 0.008814
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.6977786 -0.3022214
sample estimates:
mean of the differences 
                     -1

The results show that 1) average of $x_1 - x_2$ equals to -1; 2) $t_{12} = -3.1225$, $p = 0.008814$.

Mixed model for paired sample

We can build a mixed model for this paired-sample example. Consider the following dataset:

> y <- c(x1,x2)
> trt <- gl(2,13,26)
> subject <- gl(13,1,26)
> cbind(y, trt, subject)
       y trt subject
 [1,]  3   1       1
 [2,]  4   1       2
 [3,]  3   1       3
 [4,]  5   1       4
 [5,]  6   1       5
 [6,]  6   1       6
 [7,]  6   1       7
 [8,]  7   1       8
 [9,]  7   1       9
[10,]  6   1      10
[11,]  8   1      11
[12,]  9   1      12
[13,] 11   1      13
[14,]  3   2       1
[15,]  4   2       2
[16,]  5   2       3
[17,]  5   2       4
[18,]  8   2       5
[19,]  6   2       6
[20,]  8   2       7
[21,]  7   2       8
[22,]  9   2       9
[23,]  8   2      10
[24,]  8   2      11
[25,]  9   2      12
[26,] 14   2      13

where y is dependent variable, trt is a two-level fixed factor, and subject is experimental units considered as blocks (random factor). A mixed model can fit this dataset:

\[ y_{ij} = \tau_i + b_{j} + \varepsilon_{ij} \]

where $i=\{1,2\}$ showing two-level in fixed factor $\tau$, $j=\{1,2,\ldots , 13\}$ showing 13 subjects as a random factor $b$.

We can find the solution of this mixed model by using lme4:lmer():

> require(lme4)
> mod <- lmer(y ~ trt + (1|subject) )
> summary(mod)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ trt + (1 | subject) 

REML criterion at convergence: 98.5708 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 5.8590   2.4205  
 Residual             0.6667   0.8165  
Number of obs: 26, groups: subject, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept)   6.2308     0.7085   8.794
trt2          1.0000     0.3203   3.123

Correlation of Fixed Effects:
     (Intr)
trt2 -0.226

As you can see, the estimate of trt2 equals 1 showing the average of $x_2 - x_1$; the statistic $t=3.123$ is equivalent to the t-value I get in paired t-test.

We can find the p-value by using car:Anova():

> require(car)
> Anova(mod, type=3, test.statistic="F")
Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: y
                F Df Df.res    Pr(>F)    
(Intercept) 77.34  1 13.288 6.636e-07 ***
trt          9.75  1 12.000  0.008814 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

and then we can find that $F_{1,12} = 9.75 = 3.123^2$ and $p=0.008814$. These outcomes are all equivalent to what I got in paired t-test.