2
$\begingroup$

I am performing quite some binomial/multinomial models for my thesis. After doing the emmeans statement, I used the contrast statement to compare the emmeans of the different groups. But, I start to think about a few things:

  • My data is sampled as such that the groups don't have an equal size
  • Therefore the SE is not equal for each estimate. If my reasoning is correct, the variance between the two groups is not the same. I didn't found how the contrast test works, but if it actually just performs a t-test, then the assumptions for a t-test are violated. How can I adjust for this?

Or am I making things more difficult than they actually are?

Example:

I want to compare the distribution of 4 age classes between two methods. HD and DD, HD contain 322 individuals and DD 1711 individuals. The null hypothesis is that there is no difference in age distribution between the two methods. The Anova test on the model was significant. Therefore I want to know which age classes differ the most between the two methods.

the piece of code (Age class consists of a 0,1,2,3+ age class, method consists of the HD and DD group):

a <- multinom(`Age class` ~  `Method`  , data = goodyears)
emmeans = emmeans(a,~ `Age class`  |`Method`, mode = "prob")
z =contrast(emmeans, "pairwise", simple = "each", combine = TRUE, adjust = "mvt")

For example, the emmeans for the 1y age class in the HD group is 0.32 (95%CL [0.29,0.34]) and the emmeans for the 1y age class in the DD group is 0.24 (95%CL [0.18,0.30]). The contrast test is insignificant (logical as the intervals overlap). Is there a way to deal with the rather large difference in group size? Statistically the smaller group size caused a larger confidence interval which made the effect less significant. Or do I better accept the result as it is, and state in a discussion that the insignificance of the effect might be due to the limited sample size of one group?

$\endgroup$
4
  • $\begingroup$ Please edit the question to provide more information about the model a that you are using for the basis of the emmeans post-modeling analysis. Without that information it's hard to know if the call to emmeans actually accomplishes what you want. $\endgroup$
    – EdM
    Mar 19, 2023 at 2:39
  • $\begingroup$ Hi I provided the code, and some additional information on the variables, I hope it helps to understand the question. $\endgroup$
    – S_vdp
    Mar 19, 2023 at 7:57
  • $\begingroup$ Are you interested in the contrast between 1y age class in the HD group and the 2y age class in the DD group? Can you explain what you mean by "I want to compare the distribution of 4 age classes between the two methods." $\endgroup$
    – dipetkov
    Mar 19, 2023 at 15:16
  • $\begingroup$ Excuse me, it is a typo. I am interested in the contrast between 1y age class in HD group and the 1y age class in the DD group. $\endgroup$
    – S_vdp
    Mar 19, 2023 at 15:42

2 Answers 2

1
$\begingroup$

The emmeans package does not require equal sample sizes, so to that extent you are worrying too much. It uses the variance-covariance matrix of coefficient estimates to compare different scenarios after the initial model has been fit.

Based on that matrix, the formula for a weighted sum of correlated variables provides the variance estimate for any sum/difference of coefficient estimates involved in a comparison of interest. The assumption for significance testing is generally that the coefficient estimates follow a multivariate normal or multivariate t distribution, an assumption that underlies much statistical analysis.

The sample size enters through its effects on those variance/covariance estimates. A coefficient variance estimate typically decreases with the number of observations, so in your case there is probably more precision in estimates for the DD than in the HD class.

Where you might be making things harder for yourself is in the number of comparisons you perform. Results in emmeans are corrected for multiple comparisons. The more comparisons you perform, the harder it is to show that any single comparison is "significant." I'm not sure how many comparisons are provided by your call to the contrast() function, but I don't think that you need to perform more than an HD versus DD comparison within each of the 4 age classes to accomplish what you want, or 4 comparisons total.

Finally, don't jump to an assumption that overlap of confidence intervals necessarily means no significant difference between groups. See this answer for a more nuanced explanation.

$\endgroup$
1
$\begingroup$

Is there a way to deal with the rather large difference in group size? Statistically the smaller group size caused a larger confidence interval which made the effect less significant.

You have less information (= fewer samples) about the HD group and about the distribution of age classes within this group. So you should expect to estimate the HD parameters with less accuracy (= wider confidence intervals). To get more precise estimates, you would have to collect more data.

That said, you may be over-adjusting by computing more comparisons than you are actually interested in.

For illustration purposes, I've generated some fake data according to your specifications. There are two methods, DD and HD, and four age classes 0, 1, 2 and 3+.

Here is the output of contrast as you use it:

emm <- emmeans(fit, ~ `Age class` | `Method`, mode = "prob")
contrast(emm, "pairwise", simple = "each", combine = TRUE, adjust = "mvt")
#>  Method Age class contrast estimate     SE df t.ratio p.value
#>  DD     .         0 - 1    -0.03507 0.0170  6  -2.065  0.4246
#>  DD     .         0 - 2    -0.02747 0.0169  6  -1.630  0.6360
#>  DD     .         0 - (3+) -0.01870 0.0167  6  -1.119  0.8811
#>  DD     .         1 - 2     0.00760 0.0175  6   0.435  0.9986
#>  DD     .         1 - (3+)  0.01636 0.0173  6   0.945  0.9377
#>  DD     .         2 - (3+)  0.00877 0.0172  6   0.510  0.9968
#>  HD     .         0 - 1     0.05280 0.0477  6   1.106  0.8861
#>  HD     .         0 - 2     0.22360 0.0400  6   5.590  0.0106
#>  HD     .         0 - (3+)  0.30125 0.0351  6   8.581  0.0011
#>  HD     .         1 - 2     0.17080 0.0387  6   4.409  0.0332
#>  HD     .         1 - (3+)  0.24845 0.0340  6   7.299  0.0028
#>  HD     .         2 - (3+)  0.07765 0.0283  6   2.743  0.2039
#>  .      0         DD - HD  -0.16472 0.0291  6  -5.666  0.0101
#>  .      1         DD - HD  -0.07686 0.0285  6  -2.697  0.2156
#>  .      2         DD - HD   0.08635 0.0235  6   3.677  0.0720
#>  .      3+        DD - HD   0.15523 0.0193  6   8.054  0.0017
#> 
#> P value adjustment: mvt method for 16 tests

Why adjust for all possible 16 tests? If you are not interested in the contrast between 1y age class in the HD group and the 2y age class in the DD group and the rest of the comparisons across both ages and methods, you are over-adjusting.

Here is how to compare the two methods at each age class.

# NB: Switch the order in the `condition on` statement.
emm <- emmeans(fit, ~ `Method` | `Age class`, mode = "prob")
contrast(emm, "pairwise", adjust = "mvt")
#> Age class = 0:
#>  contrast estimate     SE df t.ratio p.value
#>  DD - HD   -0.1647 0.0291  6  -5.666  0.0013
#> 
#> Age class = 1:
#>  contrast estimate     SE df t.ratio p.value
#>  DD - HD   -0.0769 0.0285  6  -2.697  0.0357
#> 
#> Age class = 2:
#>  contrast estimate     SE df t.ratio p.value
#>  DD - HD    0.0863 0.0235  6   3.677  0.0104
#> 
#> Age class = 3+:
#>  contrast estimate     SE df t.ratio p.value
#>  DD - HD    0.1552 0.0193  6   8.054  0.0002

Note that now we get no statement about number of adjustments. Since for each age class emmeans calculates a single pairwise comparison, it applies no adjustment to the p-values.

The question if and how to adjust for multiple comparisons of interest is trickier than the fact we shouldn't calculate and adjust for comparisons of no interest. Since you are studying the difference in age distribution between the two methods (rather than, say, focusing on babies and toddlers separately), I think you should do multiple comparisons adjustment.

So finally here is how to compare the two methods at each age class and adjust for making four comparisons across the four age classes.

emm %>%
  contrast("pairwise", by = "Age class") %>%
  summary(by = NULL, adjust = "mvt")
#>  contrast Age class estimate     SE df t.ratio p.value
#>  DD - HD  0          -0.1647 0.0291  6  -5.666  0.0038
#>  DD - HD  1          -0.0769 0.0285  6  -2.697  0.1054
#>  DD - HD  2           0.0863 0.0235  6   3.677  0.0324
#>  DD - HD  3+          0.1552 0.0193  6   8.054  0.0006
#> 
#> P value adjustment: mvt method for 4 tests
$\endgroup$
1
  • $\begingroup$ Nice demonstration of the problem of more comparisons than you need! (+1) $\endgroup$
    – EdM
    Mar 19, 2023 at 21:03

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.