# Correcting p-values with the Holm Correction

In a diagnosis I’m conducting that involves more than one estimator, I have several estimators whose p-values I would like to adjust with the Holm Correction. This correction involves conducting a calculation across the p values of more than one estimator.

I presume that the best place to implement this is in the declare_diagnosand stage? But do diagnosands have access to the whole set of all estimates?

Presumably another option would be to report the mean p-value as an additional diagnosand and then conduct the Holm correction on a subset of estimators returned by diagnose_design?

1 Like

Hi Nathan,

Great question! Thanks for this. There are a few ways of doing it, but I think the method below is easiest. There is some weird code in there due to a problem with `dplyr::mutate()` stripping parameters from the `simulations` data.frame (see here, we’re working on fixing). But basically, it leverages the fact that `diagnose_design()` can also just take a data.frame and use that for diagnosis.

``````library(tidyverse)
multiple_comparisons <-
declare_population(
N = 100,
X1 = rnorm(N),
X2 = rnorm(N),
X3 = rnorm(N),
X4 = rnorm(N)
) +
declare_estimand(true_effect = 0) +
declare_estimator(X1 ~ X2, model = lm_robust, estimand = "true_effect", label = "X1 on X2") +
declare_estimator(X2 ~ X3, model = lm_robust, estimand = "true_effect", label = "X2 on X3") +
declare_estimator(X3 ~ X4, model = lm_robust, estimand = "true_effect", label = "X3 on X4") +
declare_estimator(X1 ~ X4, model = lm_robust, estimand = "true_effect", label = "X1 on X4")

simulations <-
simulate_design(multiple_comparisons, sims = 20) %>%
group_by(sim_ID) %>%
mutate(p.value = p.adjust(p.value, method = "holm")) %>%
as.data.frame()

attr(simulations,"parameters") <- data.frame(design_label = "multiple_comparisons")

diagnose_design(simulations)
``````
1 Like

Note that, if you use a reasonable number of simulations (e.g., `simulate_design(multiple_comparisons, sims = 500) %>%`), you’ll see that the power for all estimates is 2%, implying the rate of false positives is now 2% even when using an alpha of .05. In other words, the testwise rate of false positives is correctly being adjusted downwards.

1 Like

Also, this is going outside DD a bit (you could do it within DeclareDesign using a custom diagnosand), but if you want the FWER that’s pretty easy too:

``````simulations %>%
group_by(sim_ID) %>%
summarize(any_rejection = any(p.value < .05)) %>%
with(., mean(any_rejection))
``````
1 Like

Generally I would recommend adjusting p values inside a custom estimator handler - https://github.com/DeclareDesign/DeclareDesign/blob/nfultz/postestimation/vignettes/postestimation.Rmd#L225 for example. This lets you build the adjustment into the design itself instead of having to mess with simulations post hoc.

1 Like

Hi Jasper, Neal, I just realized that I failed to respond to this thread blush. Both of your responses have been super helpful, and I’m planning to come back with a followup and report of how things went in a couple of weeks.

1 Like