DeclareDesign Community

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