Introduction

Chapter 11 talks about the problem of multiple tests. These examples will be fairly short but will highlight a few ways to adjust for multiple tests.

Adjusting for Multiple Tests

  1. Let’s start by loading the tidyverse package and the rlm package.
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1.9000 ──
## ✔ ggplot2 2.2.1.9000      ✔ purrr   0.2.5      
## ✔ tibble  1.4.2.9004      ✔ dplyr   0.7.99.9000
## ✔ tidyr   0.8.1           ✔ stringr 1.3.1      
## ✔ readr   1.2.0           ✔ forcats 0.3.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(rlm)
## ── rlm 0.1.0 ───────────────────────────────────────────────────────────────────────────── learn more at tysonbarrett.com ──
## ✔ rlm attached
## ✔ No potential conflicts found
  1. Import gss into R.
data(gss)
  1. Let’s do a multiple with degree, race, and education predicting income (income06). Our first adjustment is using Bonferroni’s adjustment. To do so, we will grab the p-values from the summary() and feed that to p.adjust().
fit <- gss %>%
  lm(income06 ~ degree + race + educ, data = .) %>%
  summary()
  
p_bonf <- fit %>% 
  coef() %>%
  data.frame %>%
  .[[4]] %>%
  p.adjust("bonferroni")
fit %>%
  coef() %>%
  data.frame %>%
  .[,1:3] %>%
  cbind(p_bonf)
##                    Estimate Std..Error   t.value       p_bonf
## (Intercept)       7885.8124  4776.5311  1.650950 7.914133e-01
## degreeAssociates 29800.3711  4442.5264  6.707978 2.118946e-10
## degreeBachelors  45455.5562  4227.3359 10.752766 2.912338e-25
## degreeGraduate   57663.9272  5123.0925 11.255687 1.573189e-27
## degreeHS         20296.5103  3067.3139  6.617031 3.875839e-10
## raceOther        13466.5777  4026.0555  3.344856 6.723918e-03
## raceWhite        16327.4677  2740.7819  5.957230 2.471381e-08
## educ               706.0966   366.1085  1.928654 4.314717e-01
  1. Are there significant mean differences between the non-reference levels and the reference level regarding income for either degree or race now that we have adjusted for multiple tests? What about the effect of education?
  2. Let’s try a different approach; namely, let’s try "fdr" (the false discovery rate). This approach controls the number of false discoveries (as the name inplies). This approach is less conservative than Bonferroni’s.
p_fdr <- fit %>% 
  coef() %>%
  data.frame %>%
  .[[4]] %>%
  p.adjust("fdr")
fit %>%
  coef() %>%
  data.frame %>%
  .[,1:3] %>%
  cbind(p_fdr)
##                    Estimate Std..Error   t.value        p_fdr
## (Intercept)       7885.8124  4776.5311  1.650950 9.892666e-02
## degreeAssociates 29800.3711  4442.5264  6.707978 7.063154e-11
## degreeBachelors  45455.5562  4227.3359 10.752766 1.456169e-25
## degreeGraduate   57663.9272  5123.0925 11.255687 1.573189e-27
## degreeHS         20296.5103  3067.3139  6.617031 9.689598e-11
## raceOther        13466.5777  4026.0555  3.344856 1.120653e-03
## raceWhite        16327.4677  2740.7819  5.957230 4.942763e-09
## educ               706.0966   366.1085  1.928654 6.163882e-02
  1. Do any conclusions change based on using the false discovery rate rather than Bonferroni’s? Which change? Why?

Conclusion

Adjusting for multiple comparisons is often an important consideration. In R, there are several approaches that can be quickly applied. Although important to adjust, possibly more important is the consideration of which tests to run and the likelihood of the effect prior to making any tests.