Since Volume is heavily right-skewed, we will use the cubic root of the volume as the response variable. The resulting variable approximately resembles the diameter.
fit linear mixed model
library(lmerTest)D$group <-factor(D$group)# if we only analyze the first repetition:# fit <- lmer(diameter ~ group + (1|cage), data=D)# otherwise:fit <-lmer(diameter ~ group + repetition + (1|cage), data=D)
anova and two-sided t-test
anova(fit)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 17.386 17.386 1 42.403 6.5197 0.014352 *
repetition 36.722 36.722 1 14.215 13.7702 0.002272 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit) # look at the p-value of groupWT
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: diameter ~ group + repetition + (1 | cage)
Data: D
REML criterion at convergence: 185
Scaled residuals:
Min 1Q Median 3Q Max
-2.05849 -0.80721 0.05774 0.76071 1.86338
Random effects:
Groups Name Variance Std.Dev.
cage (Intercept) 0.2457 0.4957
Residual 2.6668 1.6330
Number of obs: 48, groups: cage, 17
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.4218 0.7860 18.6525 11.987 3.32e-10 ***
groupWT 1.2757 0.4996 42.4025 2.553 0.01435 *
repetition -1.2116 0.3265 14.2151 -3.711 0.00227 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) gropWT
groupWT -0.403
repetition -0.854 0.010
So groupWT is 1.276 bigger than groupKO. The p-value is 0.014.
one_sided_ttest <-function(fit, which_term="group", greater_than_reference=TRUE){ COEFS <-coef(summary(fit))stopifnot(any(grepl(which_term, names(fixef(fit))))) which_term <-grep(which_term, rownames(COEFS)) t_val <- COEFS[which_term, "t value"] df <- COEFS[which_term, "df"] p_val <-pt(t_val, df, lower.tail=!greater_than_reference)cat("p-value for one-sided t-test for the alternative hypothethis:\n", levels(D$group)[1], if(greater_than_reference) "<"else">", levels(D$group)[2], ":", p_val, "\n")invisible(p_val)}# PVALUE: one-sided t-test for groupWT > groupKOone_sided_ttest(fit, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
KO < WT : 0.007176008
# PVALUE: one-sided t-test for groupWT < groupKOone_sided_ttest(fit, greater_than_reference=FALSE)
p-value for one-sided t-test for the alternative hypothethis:
KO > WT : 0.992824
The one-sided p-value corresponds to half the two-sided p-value, if and only if the reference group is the one supposed to be smaller. The reference group is a term in R, and is (unless specified) the groupname that comes first in the alphabet. In this case, the reference group is groupKO. So the one-sided p-value for groupWT > groupKO is 0.007 which also equals 0.014/2.
Group Sequential Design
The consultancy case that has motivated this post, is doing a group sequential design. The idea is to stop the experiment early if the effect is very strong. The interim analysis is done after each repetition. The stopping rule is to stop if the one-sided p-value is smaller than 0,0221 (using the correction method described in Pocock (1977)).
# model with only the first repetitionfit1 <-lmer(diameter ~ group + (1|cage), data=D, subset= repetition==1)one_sided_ttest(fit1, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
KO < WT : 0.004332269
# model with the first two repetitionsfit2 <-lmer(diameter ~ group + repetition + (1|cage), data=D, subset= repetition %in%1:2)one_sided_ttest(fit2, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
KO < WT : 0.01616867
# model with all three repetitionsfit3 <-lmer(diameter ~ group + repetition + (1|cage), data=D, subset= repetition %in%1:3)one_sided_ttest(fit3, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
KO < WT : 0.007176008
The first p-value is smaller than 0.0221, so we could have stopped the experiment after the first repetition.