1
$\begingroup$

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)

enter image description here

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`[[1]] )
##           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`[1]^0.5
## [1] 3.059879
tTable[2, "t-value"]
## [1] 3.05988

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

lsm <- lsmeans::lsmeans(AOV, pairwise ~ Timepoint)
lsm$contrasts
##  contrast  estimate       SE df t.ratio p.value
##  t1 - t2  -2.559512 0.836475  2   -3.06  0.0923

Side note: Testing $\delta$ with t.test on the group means

The design being balanced, we find the same results by applying the paired $t$-test to the group means:

DTmeans <- DTsmall[, .(means=mean(y)), by=.(Timepoint,Group)]
tt <- t.test(means ~ Timepoint, data=DTmeans, paired=TRUE)
tt$statistic
##         t 
## -3.059879
tt$p.value
## [1] 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
fit$sigma^2
## [1] 0.723177

The question is: how to get the estimates of $\sigma^2_b$ and $\rho$ ?

Here are the ones given by lme:

VarCorr(fit)
## Group = pdCompSymm(0 + Timepoint) 
##             Variance StdDev    Corr 
## Timepointt1 7.421826 2.7243029      
## Timepointt2 7.421826 2.7243029 0.883
## Residual    0.723177 0.8503982

Can we get these estimates from the tables below (or from the AOV object) ?

sAOV
## 
## Error: Group
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  2  113.2   56.62               
## 
## Error: Group:Timepoint
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Timepoint  1  39.31   39.31   9.363 0.0923 .
## Residuals  2   8.40    4.20                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 18  13.02  0.7232
$\endgroup$
  • 1
    $\begingroup$ I'd like an explanation of this downvote ^^ $\endgroup$ – Stéphane Laurent Apr 15 '17 at 14:03
  • $\begingroup$ Can you show the varCorr results from fit (the small dataset) so we can see what those estimates are? $\endgroup$ – rvl Apr 16 '17 at 17:12
  • $\begingroup$ @rvl They are given at the end, just before sAOV. $\endgroup$ – Stéphane Laurent Apr 16 '17 at 17:13
  • $\begingroup$ Duh - I'm blind I guess. It's striking that the estimated rho is so high given that the true value is 0. $\endgroup$ – rvl Apr 16 '17 at 17:38
2
$\begingroup$

The model fitted by aov is of the form $$ y_{ijk} = \mu + \tau_k + G_i + \tau G_{ik} + E_{ijk} $$ which is a different parameterization of the model shown fitted using lme. The last three terms in the model are independent normal random variables, with variances $\sigma^2_G$, $\sigma^2_{\tau G}$, and $\sigma^2_E$ respectively. Consulting a standard design text, one finds that $$ E(MS(G)) = 8\sigma^2_G + 4\sigma^2_{\tau G} + \sigma_E^2 $$ $$ E(MS(\tau G)) = 4\sigma^2_{\tau G} + \sigma_E^2 $$ $$ E(MS(E)) = \sigma^2_E $$ Equating these expressions with the respective observed residual MSs for the three strata and solving, one obtains the method-of-moments estimates of these variances: $\hat\sigma^2_G = (56.62 - 4.20)/8 = 6.5525$, $\hat\sigma^2_{\tau G} = (4.20-.7232)/4 = 0.869$, $\hat\sigma^2_E=.7232$.

Under the model fitted by lme, the residual errors are the same $.7232$; but $\sigma^2_b = \sigma^2_G + \sigma^2_{\tau G}$ so that $\hat\sigma_b^2 = 6.553 + .869 = 7.42$. Meanwhile, $\rho$ is the intraclass correlation, $\hat\rho = 6.553/7.42=.882$.

$\endgroup$
  • $\begingroup$ Thank you infinitely! What is this book please? $\endgroup$ – Stéphane Laurent Apr 16 '17 at 18:29
  • 1
    $\begingroup$ Remark: lmer(y ~ Timepoint + (1|Group/Timepoint)) gives these estimates, and the Kenward-Roger number of degrees of freedom is correct. On the other hand, aov does not really fit this model, because it possibly yields a negative estimate of $\rho$, as opposed to lmer. $\endgroup$ – Stéphane Laurent Apr 16 '17 at 18:46
  • 1
    $\begingroup$ One such book is Oehlert, A First Course in Design of Experiments, avail via Creative Commons, users.stat.umn.edu/~gary/book/fcdae.pdf. See Ch.12, especially 12.6. Yes, the $\rho$ estimate can go negative, in the same cases that it is estimated as zero by lmer; then you discard offending terms and get back on the same page. $\endgroup$ – rvl Apr 16 '17 at 20:56
  • $\begingroup$ Thank you for the reference. Another question: are the numbers $8$ and $4$ universal or are they derived from $I$ and $J$ ? $\endgroup$ – Stéphane Laurent Apr 16 '17 at 22:36
  • $\begingroup$ I see, these are $2J$ and $J$. $\endgroup$ – Stéphane Laurent Apr 17 '17 at 11:28

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.