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
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")