1
$\begingroup$

I ran a multilevel binary logistic regression / generalized linear mixed-effects model in R, and then ran the following code to get post-hoc tests for a significant A x B interaction where A is a binary variable and B has 4 categories.

emms <- emmeans(model, ~ A | B)
con <- contrast(emms, interaction = "pairwise")
pairs(con, by = NULL)

I got the following output: 

contrast                estimate        SE  df z.ratio p.value
 0 - 1,X - 0 - 1,Y    1.2841638  0.4307531 Inf   2.981  0.0152
 0 - 1,X - 0 - 1,Z    0.7205958  0.4547173 Inf   1.585  0.3873
 0 - 1,X - 0 - 1,Q    0.6951286  0.4497946 Inf   1.545  0.4103
 0 - 1,Y - 0 - 1,Z   -0.5635680  0.3647502 Inf  -1.545  0.4105
 0 - 1,Y - 0 - 1,Q   -0.5890352  0.3605842 Inf  -1.634  0.3597
 0 - 1,Z - 0 - 1,Q   -0.0254672  0.3908556 Inf  -0.065  0.9999

I tried to get confidence intervals for the X v Y estimate (1.2841638) but failed.

Is it okay to do it this way?

lower_bound <- exp(1.2841638 - (1.96*0.4307531))
upper_bound <- exp(1.2841638 + (1.96*0.4307531))

Is it okay to report this? I'm reporting the confint() results for most other parameters (terms that come out of the model, and not out of emmeans post-hoc stuff) and I know that looks at slightly different confidence intervals, but I'm not sure how to get those a) manually or b) with a function out of this emmeans object.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Two things: 1. I think you want regular pairwise comparisons, not interaction contrasts. Skip that con step and do pairs(emm) directly. 2. Add type = “response”) to the emmeans call and the results will be back-transformed. OK, also 3. Look at the vignettes that come with emmeans for suggestions and examples. $\endgroup$
    – Russ Lenth
    Commented Oct 4, 2018 at 0:09
  • $\begingroup$ Thanks! I did try that first but that gives me a list of (8*7)/2=28 comparisons while what I'm interested in is how the difference between the two levels of A differs between groups (6 comparisons, like above). (Similarly to what was going on here, but with a less complex interaction: stats.stackexchange.com/questions/355611/…) $\endgroup$
    – MGy
    Commented Oct 4, 2018 at 8:44
  • $\begingroup$ To put it in different terms: I assumed that each line here in the output I'm reporting corresponds to the A x B interaction, where the levels of B are restricted to two levels at a time. $\endgroup$
    – MGy
    Commented Oct 4, 2018 at 8:48
  • $\begingroup$ As constructed, emms has B as a by variable. If you don’t mess with it, pairs(emms) will produce separate comparisons of A at each B $\endgroup$
    – Russ Lenth
    Commented Oct 4, 2018 at 12:12
  • $\begingroup$ Also look at the row labels in your output. Just because you have the number of results you expected doesn’t mean they are the right results, and in this case you have a set of 6 comparisons of comparisons of cell means. $\endgroup$
    – Russ Lenth
    Commented Oct 4, 2018 at 12:18

1 Answer 1

Reset to default
3
$\begingroup$

Here's a parallel example:

> warp = transform(warpbreaks, A=wool, B=tension)
> warp.lm = lm(log(breaks) ~ A*B, data = warp)
> emms = emmeans(warp.lm, ~ A | B)

Look at the estimates:

> emms
B = L:
 A   emmean        SE df lower.CL upper.CL
 A 3.717945 0.1246647 48 3.467290 3.968601
 B 3.282378 0.1246647 48 3.031723 3.533034

B = M:
 A   emmean        SE df lower.CL upper.CL
 A 3.116750 0.1246647 48 2.866094 3.367405
 B 3.309327 0.1246647 48 3.058671 3.559982

B = H:
 A   emmean        SE df lower.CL upper.CL
 A 3.117623 0.1246647 48 2.866967 3.368278
 B 2.904152 0.1246647 48 2.653496 3.154807

Results are given on the log (not the response) scale. 
Confidence level used: 0.95

Look at these results, back-transformed:

> summary(emms, type = "response")
B = L:
 A response       SE df lower.CL upper.CL
 A 41.17969 5.133656 48 32.04977 52.91043
 B 26.63906 3.320951 48 20.73293 34.22765

B = M:
 A response       SE df lower.CL upper.CL
 A 22.57289 2.814043 48 17.56827 29.00316
 B 27.36669 3.411661 48 21.29924 35.16256

B = H:
 A response       SE df lower.CL upper.CL
 A 22.59260 2.816501 48 17.58361 29.02849
 B 18.24975 2.275101 48 14.20361 23.44851

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

Obtain confidence intervals for comparisons within each by group, back-transformed:

> confint(pairs(emms), type = "response")
B = L:
 contrast     ratio        SE df  lower.CL upper.CL
 A / B    1.5458390 0.2725354 48 1.0844650 2.203500

B = M:
 contrast     ratio        SE df  lower.CL upper.CL
 A / B    0.8248307 0.1454198 48 0.5786502 1.175746

B = H:
 contrast     ratio        SE df  lower.CL upper.CL
 A / B    1.2379675 0.2182568 48 0.8684814 1.764648

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

Note that with a log transformation or link, back-transformed differences become ratios. Similarly, with a logit link, the comparisons will back-transform to odds ratios. For other transformations or links, there is not a sensible way to back-transform a comparison or contrast.

Finally, compare the A differences among levels of B, again back-transformed to ratios:

> confint(contrast(emms, interaction = "pairwise", by = NULL), 
+         type = "response")
 A_pairwise B_pairwise     ratio        SE df  lower.CL upper.CL
 A / B      L / M      1.8741288 0.4672755 48 1.1352279 3.093968
 A / B      L / H      1.2486911 0.3113355 48 0.7563776 2.061443
 A / B      M / H      0.6662782 0.1661228 48 0.4035889 1.099948

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

The first entry, for example, estimates the A/B ratio with B = L, divided by the A/B ratio with B = M. (It is unfortunate that A has levels A and B)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.