Context

De standaard gepaarde t-test wordt typisch gebruikt om te testen of er een significant verschil is tussen de meetwaarden van beide behandelingen. We kunnen dit probleem echter ook formuleren als een mixed model. Het voordeel hiervan is dat we een meer informatieve modeloutput krijgen, waardoor we onder andere gemakkelijker kunnen nagaan of aan de modelassumpties voldaan is. Zo veronderstellen we bij de gepaarde t-test dat het verschil tussen behandeling a en behandeling b normaal verdeeld is. Indien niet aan de assumpties voldaan is, laat het mixed model framework gemakkelijk toe om het model aan te passen, bijvoorbeeld een andere distributie kiezen voor de respons (tellingen: Poisson of negatief binomiale verdeling, binaire data of proporties: binomiale verdeling, …).

We tonen hier aan hoe de analogie tussen de gepaarde t-test en de equivalente mixed model formulering in elkaar zit.

Gebruikte packages

library(lme4)
library(pander)
library(tidyr)
library(broom)

Data

set.seed(124)

gegevens <- data.frame(
  plot = rep(1:10, 2),
  respons = c(rnorm(10), rnorm(10, 3, 1.5)),
  behandeling = rep(c("a","b"),each = 10)
)

gegevens$behandeling <- as.factor(gegevens$behandeling)
gegevens$plot <- as.factor(gegevens$plot)

# in breed formaat
gegevensbreed <- spread(gegevens, key = behandeling, value = respons)


pandoc.table(gegevens)
plot respons behandeling
1 -1.385 a
2 0.03832 a
3 -0.763 a
4 0.2123 a
5 1.426 a
6 0.7445 a
7 0.7002 a
8 -0.2294 a
9 0.1971 a
10 1.207 a
1 3.478 b
2 0.8643 b
3 2.392 b
4 4.493 b
5 4.438 b
6 4.377 b
7 2.774 b
8 1.165 b
9 1.697 b
10 1.436 b 
pandoc.table(gegevensbreed)
plot a b
1 -1.385 3.478
2 0.03832 0.8643
3 -0.763 2.392
4 0.2123 4.493
5 1.426 4.438
6 0.7445 4.377
7 0.7002 2.774
8 -0.2294 1.165
9 0.1971 1.697
10 1.207 1.436

De gepaarde t-test

ttest <- with(gegevensbreed, t.test(y = a, x = b, paired = TRUE))

pandoc.table(tidy(ttest))
estimate statistic p.value parameter conf.low conf.high
2.497 5.157 0.0005972 9 1.402 3.592

Alternatieve, maar equivalente formulering via lineair mixed model

Plot identificeert de gepaarde waarnemingen. Een random effect van plot laat toe om deze afhankelijkheid in rekening te brengen.

mm <- lmer(respons ~ behandeling + (1|plot), data = gegevens)

De parameterschatting van behandeling b geeft het verschil ten opzichte van behandeling a (= intercept), rekening houdend met de gepaardheid. Dit verschil is hetzelfde als de estimate van de gepaarde t-test.

pandoc.table(tidy(mm))
term estimate std.error statistic group
(Intercept) 0.2148 0.3711 0.5787 fixed
behandelingb 2.497 0.4841 5.157 fixed
sd_(Intercept).plot 0.4534 NA NA plot
sd_Observation.Residual 1.082 NA NA Residual

De anovatabel geeft een test voor behandeling in termen van een F-test. De t-test is gebaseerd op de t-statistiek. Beide zijn aan elkaar gerelateerd: \(F = t^2\).

pandoc.table(anova(mm, test = "F"))
  Df Sum Sq Mean Sq F value
behandeling 1 31.17 31.17 26.6 
anova(mm, test = "F")[["F value"]]
## [1] 26.59879
unname(ttest[["statistic"]]) ^ 2
## [1] 26.59879

We kunnen het betrouwbaarheidsinterval van de t-test berekenen door ons te baseren op de t-distributie.

verschil <- data.frame(
  verschil = summary(mm)$coefficients[2,1],
  se = summary(mm)$coefficients[2,2])

verschil$lwr <- verschil$verschil - qt(p = 1 - 0.05/2, df = 9) * verschil$se # overeenkomstig het BI van de t-test
verschil$upr <- verschil$verschil + qt(p = 1 - 0.05/2, df = 9) * verschil$se # overeenkomstig het BI van de t-test


pandoc.table(verschil)
verschil se lwr upr
2.497 0.4841 1.402 3.592

De geijkte procedure om een betrouwbaarheidsinterval te berekenen voor parameters van mixed models is echter met de confint() functie. Dit kan benaderend (Wald statistiek) of via profile likelihood (meer rekentijd nodig). Deze intervallen wijken lichtjes af van de output van de t-test.

Via profile likelihood:

pandoc.table(confint(mm, parm = "behandelingb", method = "profile")) 
## Computing profile confidence intervals ...
  2.5 % 97.5 %
behandelingb 1.505 3.488

Benaderend:

pandoc.table(confint(mm, parm = "behandelingb", method = "Wald"))
  2.5 % 97.5 %
behandelingb 1.548 3.446