# (Bootstrapping) Follow-Up Contrasts for Within-Subject ANOVAs (part 2)

R
statistics
ANOVA
bootstrap
code
mixed-design
permutation
Author

Mattan S. Ben-Shachar

Published

August 14, 2019

A while back I wrote a post demonstrating how to bootstrap follow-up contrasts for repeated-measure ANOVAs for cases where you data violates some / any assumptions. Here is a demo of how to conduct the same bootstrap analysis, more simply (no need to make your data wide!)

## 1. Fit your repeated-measures model with `lmer`

``````library(lme4)

data(obk.long, package = "afex") # data from the afex package

fit_mixed <- lmer(value ~ treatment * gender * phase * hour + (1|id),
data = obk.long)``````

Note that I assume here data is aggregated (one value per cell/subject), as it would be in a rmANOVA, as so it is sufficient to model only a random intercept.

## 2. Define the contrast(s) of interest

For this step we will be using `emmeans` to get the estimates of the pairwise differences between the treatment groups within each phase of the study:

``````library(emmeans)
# get the correct reference  grid with the correct ultivariate levels!
rg <- ref_grid(fit_mixed, mult.levs = rm_levels)
rg``````
``````'emmGrid' object with variables:
treatment = control, A, B
gender = F, M
phase = fup, post, pre
hour = 1, 2, 3, 4, 5``````
``````# get the expected means:
em_ <- emmeans(rg, ~ phase * treatment)``````
``NOTE: Results may be misleading due to involvement in interactions``
``em_``
`````` phase treatment emmean    SE   df lower.CL upper.CL
fup   control     4.33 0.603 13.2     3.03     5.64
post  control     4.08 0.603 13.2     2.78     5.39
pre   control     4.25 0.603 13.2     2.95     5.55
fup   A           7.25 0.661 13.2     5.82     8.68
post  A           6.50 0.661 13.2     5.07     7.93
pre   A           5.00 0.661 13.2     3.57     6.43
fup   B           7.29 0.505 13.2     6.20     8.38
post  B           6.62 0.505 13.2     5.54     7.71
pre   B           4.17 0.505 13.2     3.08     5.26

Results are averaged over the levels of: gender, hour
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95 ``````
``````# run pairwise tests between the treatment groups within each phase
c_ <- contrast(em_, "pairwise", by = 'phase')
c_``````
``````phase = fup:
contrast    estimate    SE   df t.ratio p.value
control - A  -2.9167 0.895 13.2  -3.259  0.0157
control - B  -2.9583 0.787 13.2  -3.760  0.0061
A - B        -0.0417 0.832 13.2  -0.050  0.9986

phase = post:
contrast    estimate    SE   df t.ratio p.value
control - A  -2.4167 0.895 13.2  -2.700  0.0445
control - B  -2.5417 0.787 13.2  -3.230  0.0166
A - B        -0.1250 0.832 13.2  -0.150  0.9876

phase = pre:
contrast    estimate    SE   df t.ratio p.value
control - A  -0.7500 0.895 13.2  -0.838  0.6869
control - B   0.0833 0.787 13.2   0.106  0.9938
A - B         0.8333 0.832 13.2   1.002  0.5885

Results are averaged over the levels of: gender, hour
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates ``````
``````# extract the estimates
est_names <- c("fup:  control - A", "fup:  control - B", "fup:  A - B",
"post: control - A", "post: control - B", "post: A - B",
"pre:  control - A", "pre:  control - B", "pre:  A - B")
est_values <- summary(c_)\$estimate
names(est_values) <- est_names
est_values``````
``````fup:  control - A fup:  control - B       fup:  A - B post: control - A
-2.91666667       -2.95833333       -0.04166667       -2.41666667
post: control - B       post: A - B pre:  control - A pre:  control - B
-2.54166667       -0.12500000       -0.75000000        0.08333333
pre:  A - B
0.83333333 ``````

## 3. Run the bootstrap

Now let’s wrap this all in a function that accepts the fitted model as an argument:

``````treatment_phase_contrasts <- function(mod){
rg <- ref_grid(mod, mult.levs = rm_levels)

# get the expected means:
em_ <- emmeans(rg, ~ phase * treatment)

# run pairwise tests between the treatment groups within each phase
c_ <- contrast(em_, "pairwise", by = 'phase')

# extract the estimates
est_names <- c("fup:  control - A", "fup:  control - B", "fup:  A - B",
"post: control - A", "post: control - B", "post: A - B",
"pre:  control - A", "pre:  control - B", "pre:  A - B")
est_values <- summary(c_)\$estimate
names(est_values) <- est_names
est_values
}

# test it
treatment_phase_contrasts(fit_mixed)``````
``NOTE: Results may be misleading due to involvement in interactions``
``````fup:  control - A fup:  control - B       fup:  A - B post: control - A
-2.91666667       -2.95833333       -0.04166667       -2.41666667
post: control - B       post: A - B pre:  control - A pre:  control - B
-2.54166667       -0.12500000       -0.75000000        0.08333333
pre:  A - B
0.83333333 ``````

Finally, we will use `lme4::bootMer` to get the bootstrapped estimates!

``````treatment_phase_results <-
bootMer(fit_mixed, treatment_phase_contrasts, nsim = 50) # R = 599 at least``````
``NOTE: Results may be misleading due to involvement in interactions``
``summary(treatment_phase_results) # original vs. bootstrapped estimate (bootMed)``
``````
Number of bootstrap replications R = 50
original  bootBias  bootSE   bootMed
fup:  control - A -2.916667  0.017263 0.77841 -2.801902
fup:  control - B -2.958333 -0.017880 0.86119 -3.025705
fup:  A - B       -0.041667 -0.035143 0.98850 -0.066474
post: control - A -2.416667  0.031072 0.82654 -2.383370
post: control - B -2.541667 -0.024860 0.82351 -2.520263
post: A - B       -0.125000 -0.055932 1.03670 -0.216929
pre:  control - A -0.750000 -0.065397 0.73276 -0.851533
pre:  control - B  0.083333  0.024664 0.78592  0.111930
pre:  A - B        0.833333  0.090061 0.95015  0.994195``````
``confint(treatment_phase_results, type = "perc") # does include zero?``
``````                      2.5 %     97.5 %
fup:  control - A -5.062951 -1.2782764
fup:  control - B -4.985715 -1.0325666
fup:  A - B       -2.348035  2.1660820
post: control - A -4.451445 -0.5162071
post: control - B -4.840519 -1.1705024
post: A - B       -2.349137  2.3025369
pre:  control - A -2.427992  0.8830127
pre:  control - B -1.915388  1.7159931
pre:  A - B       -1.530049  2.7527436``````

Results indicate that the Control group is lower than both treatment groups in the post and fup (follow -up) phases.

If we wanted p-values, we could use this little function (based on this demo):

``````boot_pvalues <- function(x, side = c(0, -1, 1)) {
# Based on:
# https://blogs.sas.com/content/iml/2011/11/02/how-to-compute-p-values-for-a-bootstrap-distribution.html
side <- side
x <- as.data.frame(x\$t)

ps <- sapply(x, function(.x) {
s <- na.omit(.x)
s0 <- 0
N <- length(s)

if (side == 0) {
min((1 + sum(s >= s0)) / (N + 1),
(1 + sum(s <= s0)) / (N + 1)) * 2
} else if (side < 0) {
(1 + sum(s <= s0)) / (N + 1)
} else if (side > 0) {
(1 + sum(s >= s0)) / (N + 1)
}
})

setNames(ps,colnames(x))
}

boot_pvalues(treatment_phase_results)``````
``````fup:  control - A fup:  control - B       fup:  A - B post: control - A
0.03921569        0.03921569        0.94117647        0.03921569
post: control - B       post: A - B pre:  control - A pre:  control - B
0.03921569        0.74509804        0.23529412        0.94117647
pre:  A - B
0.27450980 ``````

These p-values can then be passed to `p.adjust()` for the p-value adjustment method of your choosing.

## Summary

I’ve demonstrated (again!) how to run permutation tests on main effects / interactions, with follow-up analysis using the bootstrap method. Using this code as a basis for any analysis you might have in mind gives you all the flexibility of `emmeans`, which supports many (many) models! 