[SOLVED] When and how to apply error correction in emmeans package

When and how to apply error correction in emmeans package

I originally posted this on cross--validated but I think it might be more appropriate for SO since it's purely about software syntax.

This is a follow-up question to this post. I ran a multinomial logistic regression examining the difference in log-odds of respondents indicating they treated a range of different medical conditions (pain, sleep, mental-health/substance use (mhsu) and all other conditions (allOther)) with either licit or illicit medical cannabis.

Here is the toy data

df <- tibble(mcType = factor(rep(c("licit", "illicit"),
                                 times = c(534,1207))),
             cond = factor(c(rep(c("pain","mhsu","allOther","sleep"), 
                                 times = c(280,141,82,31)),
                             rep(c("pain","mhsu","allOther","sleep"), 
                                 times = c(491,360,208,148))),
                           levels = c("pain","sleep","mhsu","allOther")))

And the proportions of each type of condition reported for each type of cannabis

mcType  cond         n   tot  perc
<fct>   <fct>    <int> <int> <dbl>
1 illicit pain       491  1207 40.7 
2 illicit sleep      148  1207 12.3 
3 illicit mhsu       360  1207 29.8 
4 illicit allOther   208  1207 17.2 
5 licit   pain       280   534 52.4 
6 licit   sleep       31   534  5.81
7 licit   mhsu       141   534 26.4 
8 licit   allOther    82   534 15.4

To see whether there were differences in the relative proportion of respondents indicating each type of condition based on the type of cannabis they report using I ran a multinomial logistic regression using multinom() in the nnet package. Output below,

library(nnet)
summary(mm <- multinom(cond ~ mcType,
                       data = df))


# output
Coefficients:
  (Intercept) mcTypelicit
sleep     -1.1992431  -1.0014884
mhsu      -0.3103369  -0.3756443
allOther  -0.8589398  -0.3691759

Std. Errors:
  (Intercept) mcTypelicit
sleep     0.09377333   0.2112368
mhsu      0.06938587   0.1244098
allOther  0.08273132   0.1503720

Residual Deviance: 4327.814 
AIC: 4339.814

The I ran tests of simple effects, using the emmeans package. In this blog post the author suggests that the emmeans package applies error correction by default, but that you can control this via the adjust = argument.

# testing effect of mc type at each level of condition. first create emmeans object
library(emmeans)
(em_mcTypeByCond <- emmeans(object = mm,
                            specs = ~mcType|cond,
                            adjust = "bonferroni"))

# output  
cond = pain:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.4068 0.01414  6   0.3648   0.4488
 licit   0.5243 0.02161  6   0.4602   0.5885

cond = sleep:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.1226 0.00944  6   0.0946   0.1506
 licit   0.0581 0.01012  6   0.0280   0.0881

cond = mhsu:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.2983 0.01317  6   0.2592   0.3374
 licit   0.2641 0.01908  6   0.2074   0.3207

cond = allOther:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.1723 0.01087  6   0.1401   0.2046
 licit   0.1535 0.01560  6   0.1072   0.1999

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates

The problem is that I don't seem to be able to choose any other method of error correction (e.g. "BH", "fdr", "westfall", "holm"). I am not sure if it is because I am applying the correction at the wrong step, i.e. before I apply any tests.

So I tried applying the adjust argument within the pairs() function (testing the difference in probability of each condition between the two types of cannabis)

(mcTypeByCond_test <- pairs(em_mcTypeByCond,
                            adjust = "bonferroni"))

cond = pain:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit  -0.1175 0.0258  6 -4.551  0.0039 

cond = sleep:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0646 0.0138  6  4.665  0.0034 

cond = mhsu:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0342 0.0232  6  1.476  0.1905 

cond = allOther:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0188 0.0190  6  0.987  0.3616

But as you can see this does not provide any message telling me what type of error correction was applied (I assume none, and tried several different methods). Also I want to control error across all four pairwise comparisons.

So I need to know how and at what stage I need to make the arguments specifying adjustment of p-values.

Any help much appreciated

Solution

P-value adjustments are applied to each by group, and there is only one comparison - hence no multiplicity - in each group. And no annotation about adjustments is shown when no adjustments are made.

To apply an adjustment to all the results, you need to remove the by variable from consideration when displaying the results:

summary(pairs(...), by = NULL, adjust = "bonf")