Mixed model formulation of a paired t-test
Learn some insights about mixed models based on a simple example
What you will learn
In this tutorial we explain the analogy between the paired t-test and the corresponding mixed model formulation.
Used packages
library(knitr)
library(lme4)
library(tidyr)
library(broom)
library(DHARMa)
Data
- plot: identifies paired measurements
- response: measurement values
- treatment: identifies two treatments (a and b)
set.seed(124)
paired_data <- data.frame(
plot = rep(1:10, 2),
response = c(rnorm(10), rnorm(10, 3, 1.5)),
treatment = rep(c("a", "b"), each = 10)
)
paired_data$treatment <- as.factor(paired_data$treatment)
paired_data$plot <- as.factor(paired_data$plot)
# in wide format
paired_data_wide <- pivot_wider(
paired_data,
id_cols = plot,
names_from = treatment,
values_from = response
)
kable(paired_data)
| plot | response | treatment |
|---|---|---|
| 1 | -1.3850706 | a |
| 2 | 0.0383232 | a |
| 3 | -0.7630302 | a |
| 4 | 0.2123061 | a |
| 5 | 1.4255380 | a |
| 6 | 0.7444798 | a |
| 7 | 0.7002294 | a |
| 8 | -0.2293546 | a |
| 9 | 0.1970939 | a |
| 10 | 1.2071538 | a |
| 1 | 3.4775051 | b |
| 2 | 0.8643017 | b |
| 3 | 2.3923637 | b |
| 4 | 4.4930798 | b |
| 5 | 4.4382267 | b |
| 6 | 4.3771318 | b |
| 7 | 2.7735456 | b |
| 8 | 1.1653968 | b |
| 9 | 1.6967636 | b |
| 10 | 1.4362720 | b |
kable(paired_data_wide)
| plot | a | b |
|---|---|---|
| 1 | -1.3850706 | 3.4775051 |
| 2 | 0.0383232 | 0.8643017 |
| 3 | -0.7630302 | 2.3923637 |
| 4 | 0.2123061 | 4.4930798 |
| 5 | 1.4255380 | 4.4382267 |
| 6 | 0.7444798 | 4.3771318 |
| 7 | 0.7002294 | 2.7735456 |
| 8 | -0.2293546 | 1.1653968 |
| 9 | 0.1970939 | 1.6967636 |
| 10 | 1.2071538 | 1.4362720 |
The paired t-test
ttest <- with(
paired_data_wide,
t.test(y = a, x = b, paired = TRUE)
)
kable(tidy(ttest))
| estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|
| 2.496692 | 5.157401 | 0.0005972 | 9 | 1.401584 | 3.591799 | Paired t-test | two.sided |
Alternative, but equivalent formulation via a lineair mixed model
Plot identifies the paired measurements. A random effect for plot allows us to take this dependence into account.
mm <- lmer(response ~ treatment + (1 | plot),
data = paired_data
)
The parameter estimates for treatment b gives the difference compared to treatment a (= intercept), accounting for the paired nature of the data. This difference is the same as the estimate for the paired t-test.
kable(tidy(mm))
| term | estimate | std.error | statistic | group |
|---|---|---|---|---|
| (Intercept) | 0.2147669 | 0.3711260 | 0.578690 | fixed |
| treatmentb | 2.4966918 | 0.4840988 | 5.157401 | fixed |
| sd_(Intercept).plot | 0.4534165 | NA | NA | plot |
| sd_Observation.Residual | 1.0824778 | NA | NA | Residual |
The anova output gives us a test for treatment in terms of an F-test. The t-test is based on the t-statistic. Both test statistics are related: (F = t^2).
kable(anova(mm))
| npar | Sum Sq | Mean Sq | F value | |
|---|---|---|---|---|
| treatment | 1 | 31.16735 | 31.16735 | 26.59879 |
anova(mm)[["F value"]]
## [1] 26.59879
unname(ttest[["statistic"]])^2
## [1] 26.59879
We can calculate the confidence interval given as part of the t-test output, based on the t-distributie.
difference <- data.frame(
diff = summary(mm)$coefficients[2, 1],
se = summary(mm)$coefficients[2, 2]
)
difference$lwr <- difference$diff - qt(p = 1 - 0.05 / 2, df = 9) * difference$se
difference$upr <- difference$diff + qt(p = 1 - 0.05 / 2, df = 9) * difference$se
kable(difference)
| diff | se | lwr | upr |
|---|---|---|---|
| 2.496692 | 0.4840988 | 1.401584 | 3.591799 |
The recommended procedure to calculate a confidence interval for
parameters of mixed models is, however, to use the `confint function.
Either an approximation (Wald statistic) or a profile likelihood
confidence interval can be calculated. These intervals are slightly
different from the t-distribution based confidence interval.
Via profile likelihood:
kable(confint(mm, parm = "treatmentb", method = "profile"))
## Computing profile confidence intervals ...
| 2.5 % | 97.5 % | |
|---|---|---|
| treatmentb | 1.505743 | 3.487635 |
Wald-type confidence interval:
kable(confint(mm, parm = "treatmentb", method = "Wald"))
| 2.5 % | 97.5 % | |
|---|---|---|
| treatmentb | 1.547876 | 3.445508 |
Were model assuptions met? Yes.
DHARMa::plotQQunif(mm)
Take home message
The standard paired t-test is typically used to test for a significant differences between two paired treatments. We can formulate the test in terms of a mixed model. The benefit is that we get more informative model output, which allows us among other things to check if model assumptions were met. For the paired t-test, one assumption is that the paired differences between treatments follow a normal distribution. When these assumptions are not met, the flexibility of the mixed model framework allows to improve the model to better fit the requirements for the data at hand. For instance, one can choose from a number of parametric statistical distributions that are likely to fit the data (for counts, the Poisson or negative binomial distribution can be chosen, and for binary or proportional data, a binomial distribution is an obvious choice).