Fitting the paired t-test with replicates and obtaining the explicit variance estimates: lme() vs. aov() in R

Summary of the question

I will show how to fit a certain model in R with lme and how to test the nullity of the parameter of interest $\delta$ with aov. The results are equivalent, up to the wrong number of degrees of freedom given by lme (and assuming lme does not return an almost null estimate of the between-group variance $\sigma^2_b$).

The lme fitting provides estimates of the variance components of the model. My question is how to find, if possible, these estimates from the results of aov.

Data and model

The dataset is simulated by the function SimData below

library(data.table)
SimData <- function(I, J, delta, sigmab, sigmaw, rho){
Mu <- setNames(c(0, delta), c("t1", "t2"))
Sigma <- rbind(
c(sigmab^2, rho*sigmab^2),
c(rho*sigmab^2, sigmab^2)
)
### simulation between-groups ###
mu <- mvtnorm::rmvnorm(I, Mu, Sigma)
### simulation within-groups ###
y1 <- c(vapply(mu[,"t1"], function(x) rnorm(J, x, sigmaw),
FUN.VALUE = numeric(J)))
y2 <- c(vapply(mu[,"t2"], function(x) rnorm(J, x, sigmaw),
FUN.VALUE = numeric(J)))
### constructs the dataset ####
Timepoint <- gl(2L, I*J, labels=c("t1", "t2"))
Group <- as.factor(paste0("grp", rep(gl(I,J), times=2L)))
dat <- data.table(
Timepoint = Timepoint,
Group     = Group,
y         = c(y1, y2),
key = c("Timepoint", "Group")
)
return(dat)
}

Below is a plot of a small simulated dataset:

library(ggplot2)
set.seed(105L)
DTsmall <- SimData(I=3L, J=4L, delta=0, sigmab=3, sigmaw=1, rho=0)
ggplot(DTsmall, aes(x=Timepoint, y=y, color=Group)) + geom_point(size=2) Say that the observations are measurements taken on $I$ groups at two timepoints, and $J$ measurements are recorded for each group at each timepoint. We make the assumptions:

• the within-group variance is common to the $I$ groups at the two timepoints;

• the between-group variance is common to the two timepoints.

We use the indexes $i$, $j$ and $k$ to respectively denote the timepoint, the group and the observation.

We assume a correlation between the records of a same group taken at the two timepoints.

Mathematically, the theoretical pairs of means $(\mu_{1j}, \mu_{2j})$ of the groups are random effects following a bivariate normal distribution: $$\begin{pmatrix} \mu_{1j} \\ \mu_{2j} \end{pmatrix} \sim_{\text{iid}} {\cal N}\left(\begin{pmatrix} \mu_{1} \\ \mu_{2} \end{pmatrix}, \begin{pmatrix} \sigma^2_{b} & \rho\sigma_{b}^2 \\ \rho\sigma_{b}^2 & \sigma^2_{b} \end{pmatrix} \right),$$ centered around the theoretical pair of means $(\mu_1, \mu_2)$ at the two timepoints. Then one assumes that for each timepoint $i$, the observations follow a normal distribution within each group $j$, with, as said before, a common within-variance $\sigma^2_{w}$: $$(y_{ijk} \mid \mu_{ij}) \sim_{\text{iid}} {\cal N}(\mu_{ij}, \sigma^2_{w}).$$

We set $\mu_1=0$ in the function SimData. We are interested in $\delta:=\mu_2-\mu_1$.

Fitting the model with lme

Firstly, let's simulate a dataset with high values of $I$ and $J$ in order to check that the lme fitting correctly estimates the true parameters.

set.seed(31415L)
DTbig <- SimData(I=300, J=200, delta=1, sigmab=3, sigmaw=1, rho=0.5)
library(nlme)
fit_DTbig <- lme(y~Timepoint, data=DTbig,
random= list(Group= pdCompSymm(form= ~0+Timepoint)))

We find good estimates of $\sigma_b$, $\sigma_w$ and $\rho$:

VarCorr(fit_DTbig)
## Group = pdCompSymm(0 + Timepoint)
##             Variance StdDev   Corr
## Timepointt1 9.647973 3.106119
## Timepointt2 9.647973 3.106119 0.522
## Residual    1.005556 1.002774

As well as good estimates of $\mu_1$ and $\delta$:

summary(fit_DTbig)$tTable ## Value Std.Error DF t-value p-value ## (Intercept) 0.04328747 0.1793786 119699 0.2413191 8.093082e-01 ## Timepointt2 0.99384798 0.1754768 119699 5.6637011 1.484807e-08 We will use the small dataset from now on: fit <- lme(y~Timepoint, data=DTsmall, random= list(Group= pdCompSymm(form= ~0+Timepoint))) sfit <- summary(fit) ( tTable <- sfit$tTable )
##                 Value Std.Error DF    t-value    p-value
## (Intercept) -1.350037 1.5919192 20 -0.8480561 0.40644218
## Timepointt2  2.559512 0.8364745 20  3.0598803 0.00618172

Recall that we are interested in the second row: the inference on $\delta$.

Note: The number of degrees of freedom in the above table is not correct; it should be $I-1$. Consequently the $p$-value is not correct.

Testing $\delta$ with aov

AOV <- aov(y~ Timepoint + Error(Group/Timepoint), data=DTsmall)
sAOV <- summary(AOV)
( fTable <- sAOV$Error: Group:Timepoint[] ) ## Df Sum Sq Mean Sq F value Pr(>F) ## Timepoint 1 39.307 39.307 9.3629 0.09226 . ## Residuals 2 8.396 4.198 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 This$F$-value is the square of the$t$-value given by lme: fTable$F value^0.5
##  3.059879
tTable[2, "t-value"]
##  3.05988

We can get the $t$-table with lsmeans:

lsm <- lsmeans::lsmeans(AOV, pairwise ~ Timepoint)
##         t
## -3.059879
tt$p.value ##  0.09226219 Question We can get the estimate of$\sigma^2_w$from the aov results: sAOV$Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 18  13.02  0.7232