library(lme4)
# read in data;
# see get_emergence_data.R
d <- readRDS(file = "emergence.rds")
# fit logistic regression model with experiment as random effect
m <- glmer(emergence ~ day * treatment + (1|experiment),
data = d, family = binomial(link = "probit"),
weights = total)Examples of emmeans
Intro
In this document I show how to use {emmeans} after fitting a model to make various comparisons.
Fit a model
The interaction certainly seems significant.
car::Anova(m)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: emergence
Chisq Df Pr(>Chisq)
day 1596.97 1 < 2.2e-16 ***
treatment 1628.81 3 < 2.2e-16 ***
day:treatment 724.07 3 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We can visualize the model using {ggeffects}.
library(ggeffects)
ggeffect(m, terms = c("day[0:16]", "treatment"),
bias_correction = TRUE) |> plot()
Use {emmeans} to investigate model
Some examples of using emmeans to quantify relationships we see in the plot.
1. Compare Cl to Combo at day 16
Specify pairwise ~ treatment | day to make all pairwise comparisons, but them set at = list(day = 16, treatment = c("Cl", "Combo")) to restrict comparisons to Cl and Combo at Day 16. We set regrid = "response" to make comparison in terms of probabilities instead of probits.
library(emmeans)
emmeans(m, specs = pairwise ~ treatment | day,
at = list(day = 16, treatment = c("Cl", "Combo")),
regrid = "response") $emmeans
day = 16:
treatment prob SE df asymp.LCL asymp.UCL
Cl 0.415 0.108 Inf 0.20321 0.628
Combo 0.117 0.055 Inf 0.00939 0.225
Confidence level used: 0.95
$contrasts
day = 16:
contrast estimate SE df z.ratio p.value
Cl - Combo 0.298 0.0563 Inf 5.296 <.0001The output says that Cl is 0.415 at Day 16 and that Combo is 0.117 at Day 16. The difference between them (under contrasts) is about 0.298. That difference appears to reliably positive (ie, “significant”). We can get confidence intervals by piping the result into confint().
emmeans(m, specs = pairwise ~ treatment | day,
at = list(day = 16, treatment = c("Cl", "Combo")),
regrid = "response") |>
confint()$emmeans
day = 16:
treatment prob SE df asymp.LCL asymp.UCL
Cl 0.415 0.108 Inf 0.20321 0.628
Combo 0.117 0.055 Inf 0.00939 0.225
Confidence level used: 0.95
$contrasts
day = 16:
contrast estimate SE df asymp.LCL asymp.UCL
Cl - Combo 0.298 0.0563 Inf 0.188 0.409
Confidence level used: 0.95 The 95% CI is reported as [0.19,0.41]. The difference in Cl - Combo appears to be at least 0.19, perhaps as high as 0.41. All data in the CI range is plausible. The estimated difference of 0.298 is but one value. Instead of stating with certainty the difference is 0.298 we might say the difference appears to be in the range of [0.19, 0.41].
2. Compare Day 10 to Day 16 for all treatments
Is there a difference in emergence between Day 10 and Day 16 within each treatment? The syntax revpairwise ~ day | treatment says compare days within treatments in “reverse” order, which in this case means day 16 - day 10. By setting at = list(day = c(10, 16)) we make predictions at just day 10 and 16. This returns predicted probabilities at days 10 and 16 as well as the difference in probabilities between those days. Each difference is also tested against the null of no difference between days.
emmeans(m, specs = revpairwise ~ day | treatment,
at = list(day = c(10, 16)),
regrid = "response")$emmeans
treatment = Ag:
day prob SE df asymp.LCL asymp.UCL
10 0.1048 0.0501 Inf 0.00658 0.203
16 0.1925 0.0762 Inf 0.04317 0.342
treatment = Cl:
day prob SE df asymp.LCL asymp.UCL
10 0.1998 0.0771 Inf 0.04859 0.351
16 0.4154 0.1080 Inf 0.20321 0.628
treatment = Combo:
day prob SE df asymp.LCL asymp.UCL
10 0.0643 0.0347 Inf -0.00382 0.132
16 0.1172 0.0550 Inf 0.00939 0.225
treatment = Control:
day prob SE df asymp.LCL asymp.UCL
10 0.5326 0.1100 Inf 0.31768 0.747
16 0.9475 0.0299 Inf 0.88890 1.006
Confidence level used: 0.95
$contrasts
treatment = Ag:
contrast estimate SE df z.ratio p.value
day16 - day10 0.0876 0.0271 Inf 3.239 0.0012
treatment = Cl:
contrast estimate SE df z.ratio p.value
day16 - day10 0.2157 0.0327 Inf 6.589 <.0001
treatment = Combo:
contrast estimate SE df z.ratio p.value
day16 - day10 0.0529 0.0211 Inf 2.509 0.0121
treatment = Control:
contrast estimate SE df z.ratio p.value
day16 - day10 0.4149 0.0801 Inf 5.180 <.0001Again we can pipe into confint() for confidence intervals.
emmeans(m, specs = revpairwise ~ day | treatment,
at = list(day = c(10, 16)),
regrid = "response") |>
confint()$emmeans
treatment = Ag:
day prob SE df asymp.LCL asymp.UCL
10 0.1048 0.0501 Inf 0.00658 0.203
16 0.1925 0.0762 Inf 0.04317 0.342
treatment = Cl:
day prob SE df asymp.LCL asymp.UCL
10 0.1998 0.0771 Inf 0.04859 0.351
16 0.4154 0.1080 Inf 0.20321 0.628
treatment = Combo:
day prob SE df asymp.LCL asymp.UCL
10 0.0643 0.0347 Inf -0.00382 0.132
16 0.1172 0.0550 Inf 0.00939 0.225
treatment = Control:
day prob SE df asymp.LCL asymp.UCL
10 0.5326 0.1100 Inf 0.31768 0.747
16 0.9475 0.0299 Inf 0.88890 1.006
Confidence level used: 0.95
$contrasts
treatment = Ag:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.0876 0.0271 Inf 0.0346 0.1407
treatment = Cl:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.2157 0.0327 Inf 0.1515 0.2798
treatment = Combo:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.0529 0.0211 Inf 0.0116 0.0942
treatment = Control:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.4149 0.0801 Inf 0.2579 0.5719
Confidence level used: 0.95 3. Contrasts of contrasts
In the previous comparison, Ag and Combo had a difference on 0.0876 and 0.0529, respectively, between days 10 and 16. Are those differences significantly different?
This requires a little extra work. First we use emmeans as before without revpairwise and save the result.
emg <- emmeans(m, specs = ~ day | treatment,
at = list(day = c(10, 16)),
regrid = "response")
emgtreatment = Ag:
day prob SE df asymp.LCL asymp.UCL
10 0.1048 0.0501 Inf 0.00658 0.203
16 0.1925 0.0762 Inf 0.04317 0.342
treatment = Cl:
day prob SE df asymp.LCL asymp.UCL
10 0.1998 0.0771 Inf 0.04859 0.351
16 0.4154 0.1080 Inf 0.20321 0.628
treatment = Combo:
day prob SE df asymp.LCL asymp.UCL
10 0.0643 0.0347 Inf -0.00382 0.132
16 0.1172 0.0550 Inf 0.00939 0.225
treatment = Control:
day prob SE df asymp.LCL asymp.UCL
10 0.5326 0.1100 Inf 0.31768 0.747
16 0.9475 0.0299 Inf 0.88890 1.006
Confidence level used: 0.95 Now we use the pairs() function to create the first set of contrasts by “day”. This is where the previous example ended
pairs(emg, simple = "day", reverse = TRUE) |> confint()treatment = Ag:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.0876 0.0271 Inf 0.0346 0.1407
treatment = Cl:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.2157 0.0327 Inf 0.1515 0.2798
treatment = Combo:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.0529 0.0211 Inf 0.0116 0.0942
treatment = Control:
contrast estimate SE df asymp.LCL asymp.UCL
day16 - day10 0.4149 0.0801 Inf 0.2579 0.5719
Confidence level used: 0.95 Now we use pairs() again to compare those contrasts, except this time we set by = NULL.
pairs(pairs(emg, simple = "day", reverse = TRUE), by = NULL) contrast estimate SE df z.ratio
(day16 - day10 Ag) - (day16 - day10 Cl) -0.1280 0.0166 Inf -7.714
(day16 - day10 Ag) - (day16 - day10 Combo) 0.0347 0.0138 Inf 2.510
(day16 - day10 Ag) - (day16 - day10 Control) -0.3273 0.1060 Inf -3.096
(day16 - day10 Cl) - (day16 - day10 Combo) 0.1628 0.0183 Inf 8.888
(day16 - day10 Cl) - (day16 - day10 Control) -0.1993 0.1110 Inf -1.796
(day16 - day10 Combo) - (day16 - day10 Control) -0.3620 0.0999 Inf -3.625
p.value
<.0001
0.0583
0.0106
<.0001
0.2754
0.0016
P value adjustment: tukey method for comparing a family of 4 estimates In row 2 we see the difference in differences between Day 10 and Day 16 for Ag vs Combo is marginally significant at 0.0583. There is possibly some evidence of difference but we probably shouldn’t rule out the possibility that a different sample could yield different results.