1
$\begingroup$

I just did a series of 40 t-tests and then proceeded to use bonferroni correction for multiple testing and my P values reduced in size. Why does this happen? I was my impression that multiple tests correction would always result in an increase in p value size.

Before correction:

           substrate              p_value
19           glycan 0.000139091904182433
23 dermatan sulfate 0.000139091904182433
4            chitin 0.000294435140367691
22           xylose  0.00387472305660014
5       beta-glucan  0.00552400891530821
2         cellulose   0.0130881279666714

After correction:

          substrate      p_value
10           glucose 5.110415e-21
29          fructose 1.745709e-20
26            lignin 7.090204e-18
30 cyclomaltodextrin 3.569263e-10
31  lacto-N-tetraose 3.569263e-10
32       hyaluronate 3.569263e-10

Code used to generate the P values:

# loop labeling one substrate as "A" and every other substrate as "B" then doing a T test between then counts
sub_pvals = NULL
for(sub in unique(cazy_cata_melt$Substrate)){
  df= cazy_cata_melt
  df[df$Substrate != sub,]$Substrate = "B"
  df[df$Substrate == sub,]$Substrate = "A"
  input = cbind(substrate = sub, p_value = t.test(value ~ Substrate, data = df)[[3]][1])
  sub_pvals = rbind.data.frame(sub_pvals, input)
}

#correction for multiple testing
    sub_pvals$p_value = p.adjust(sub_pvals$p_value, method = "bonferroni", n = length(unique(cazy_cata_melt$Substrate)))

#ordering the dataframe
sub_pvals = sub_pvals[order(sub_pvals$p_value),]

Data available here: https://pastebin.com/vsbYGkQW

Improve this question
$\endgroup$
8
  • $\begingroup$ That should never happen. If you post the code you used to run the correction, it'll make it much easier to figure out what went wrong. And something did go wrong - not only are the p-values getting more significant, they're doing so by a factor of billions or more, which shouldn't happen in either direction with only 40 tests to correct for. This result would be strange even with the before/after values switched. $\endgroup$
    – Nuclear Hoagie
    Jul 24, 2020 at 14:50
  • $\begingroup$ Hi @NuclearWang I have added the code and will add the dataset also. EDIT: data added $\endgroup$
    – Lamma
    Jul 24, 2020 at 14:55
  • 1
    $\begingroup$ The before and after lists are of different substrates so...it's hard to say from your post. $\endgroup$
    – Matt Krause
    Jul 24, 2020 at 14:59
  • $\begingroup$ @MattKrause That is just because I have ordered to show smallest P values $\endgroup$
    – Lamma
    Jul 24, 2020 at 15:01
  • 1
    $\begingroup$ The correction shouldn't affect the ordering by p-values, so something is clearly wrong with the coding. $\endgroup$
    – EdM
    Jul 24, 2020 at 15:02

1 Answer 1

Reset to default
0
$\begingroup$

The underlying data are counts associated with each of a set of 40 substrates, so putting aside the coding problem (which isn't really on-topic here) the approach has two problems: t-tests aren't correct for count data, and serially comparing each substrate against the mean of all the other substrates isn't a correct way to do these tests.

Count data are best analyzed as Poisson or negative-binomial data. This is possible for example with the glm() function in R. In your case that would be set up similarly to an ANOVA, coding the substrates as levels of a single categorical predictor. The analysis would be performed with an underlying error distribution (needed to asses the significance of any differences) appropriate to count data, for which the normal distribution assumed by ANOVA and t-tests doesn't hold.

You start with the significance of the overall model. If the model as a whole isn't significant you stop and don't proceed to individual comparisons. If the model is significant overall, there are much better (and more powerful) ways to examine differences among the individual substrates; see this answer for an example with count data.

Improve this answer
$\endgroup$
4
  • $\begingroup$ How do i determine if a model is significant and then how do I decide is reporting estimated marginal means is the right thing to do? $\endgroup$
    – Lamma
    Aug 2, 2020 at 8:23
  • $\begingroup$ @Lamma if you run a model appropriate for count data (Poisson, quasi-Poisson, negative binomial) of the form counts ~ substrate with substrate a multi-category factor variable, you will get a p-value for the entire model that tests whether there are any significant differences among substrates with respect to counts. If that test shows significance, then examining estimated marginal means (e.g., via the emmeans package) provides principled ways to examine differences among the substrates. That pools information from all substrates in a way that's superior to repeated simple t-tests. $\endgroup$
    – EdM
    Aug 2, 2020 at 14:25
  • $\begingroup$ Thank you this makes sense :) I am still a bit stuck on how to interpret the emmeans though. I have looked over the package "basic of emmeans" it still don't really get hem $\endgroup$
    – Lamma
    Aug 2, 2020 at 15:35
  • $\begingroup$ @Lamma see the end of my answer here, based on this comment on the original version of that answer by emmeans author Russ Lenth: "Look at the documentation in emmeans for ‘contrast()‘ and eff.emmc. In short, specifying “eff” contrasts gives comparisons of each mean with the grand mean." I think that's what you want. $\endgroup$
    – EdM
    Aug 3, 2020 at 15:31

Your Answer

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

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