0

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

llewmills
  • 2,959
  • 3
  • 31
  • 58

1 Answers1

1

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")
Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • Brilliant thank you so much. I never would have worked that out on my own. – llewmills Sep 24 '21 at 16:05
  • 1
    I think it's in the FAQs vignette. – Russ Lenth Sep 24 '21 at 17:29
  • 1
    Turns out it isn't in FAQs, so I added it. Will be in next update to **emmeans** – Russ Lenth Sep 25 '21 at 02:07
  • @RussLenth could you elaborate on this topic? This seems counterintuitive to me. We are estimating marginal means from the model, then comparing multiple of those marginal means. Why would this not be considered multiple testing, regardless of the groupings? – creutzml Jun 15 '22 at 22:44
  • Well, different people make different choices. But I think this is pretty conventionally done. And in the case of pairwise comparisons of means in each of more than one groups, the Tukey adjustment is available *only* when you adjust each group separately. Otherwise the mathematics of the Tukey adjustment don't apply. – Russ Lenth Jun 15 '22 at 23:41
  • Adding a bit more. Context is obviously important -- is it exploratory, are we in a regulatory setting, what are the risks, etc. I think it may even depend on the group size. When it is 2, as in this example, we have four comparisons and it seems we could be asking, "in which cond is there an mcType effect?" -- which would call for a multiplicity adjustment on the set of 4 tests. If we reverse the factors, though, there are 6 pairwise comparisons of cond with each mcType. And it feels like we are more likely to view those as two separate families of comparisons rather than one family of 12. – Russ Lenth Jun 16 '22 at 00:23