[Stata] Multiple Hypothesis Testing: wyoung and rwolf2
Social work and health studies often test one exposure against several outcomes: access to care, use of services, insurance, distress. Each test has its own 5 percent chance of a false positive, so with six outcomes the chance that at least one null effect looks significant is well above 5 percent. A family-wise error rate correction adjusts the p-values so that this chance stays at 5 percent across the whole set of tests.
Bonferroni and Holm corrections are easy to compute by hand but conservative, because they ignore the correlation between outcomes. Two resampling methods use that correlation and lose less power: the Westfall-Young stepdown procedure and the Romano-Wolf stepdown procedure. In Stata they are available as wyoung by Julian Reif, used in Jones, Molitor, and Reif (2019), and as rwolf and rwolf2 by Damian Clarke, described in Clarke, Romano, and Wolf (2020).
The example tests whether adults with limited English proficiency differ from other adults on six outcomes.
Installation
ssc install wyoung
ssc install rwolf2The current wyoung (2.0, December 2024) requires Stata 13 or newer. rwolf2 installs alongside the original rwolf and has a more flexible syntax, one parenthesized model per hypothesis.
Data
The examples use the 2024 California Health Interview Survey (CHIS) adult public use file, which has 24,810 adult respondents. The CHIS public use files are free and available through a brief application on the UCLA Center for Health Policy Research site. Every post in this series builds its variables with the same short do-file, 00_chis_prep.do:
use "chis_adult_2024.dta", clear
* outcomes
gen k6 = distress if distress >= 0 // K6 score, past month (0-24)
gen nousual = usual == 2 // no usual source of care
gen delay = ah22 == 1 // delayed or did not get medical care, past 12 months
gen erv = er == 1 // any ER visit, past 12 months
gen uninsured = ins == 2
gen needhelp = af81 == 1 if af81 > 0 // needed help for emotional/mental or substance problem
gen fairpoor = ab1 >= 4 // fair or poor self-rated health
gen spd = dstrs12 == 1 // likely serious psychological distress, past year
* groups
gen lep = spk_eng == 3 // speaks English not well or not at all
gen foreign = ah33new == 2
gen female = srsex == 2
* covariates
gen age = srage_p1
recode ahedc_p1 (1 2 = 1 "Less than high school") (3 = 2 "High school") ///
(4 5 6 = 3 "Some college") (7 8 9 = 4 "Bachelor's or more"), gen(educ)
gen pov = povll
label define povllf 4 "300%+ FPL", modify
label values pov povllf
recode ombsrr_p1 (1 = 1 "Latino") (2 = 2 "White") (3 = 3 "Black") (5 = 4 "Asian") ///
(4 6 = 5 "Other or multiracial"), gen(race)
* CHIS design: final weight plus 80 jackknife replicate weights
svyset [pw=rakedw0], jkrweight(rakedw1-rakedw80, multiplier(1)) vce(jackknife) mseLimited English proficiency (lep) means the respondent reported speaking English not well or not at all; 1,520 respondents (6.1 percent) meet that definition. Commands that accept the svy prefix use the jackknife replicate weights; commands that do not accept it use the final weight rakedw0 as a probability weight, which gives the right point estimates and robust standard errors that ignore the replicate design.
Six outcomes, one exposure
The outcomes are having no usual source of care, delaying or forgoing medical care, any emergency room visit, being uninsured, needing help for a mental health or substance use problem, and likely serious psychological distress in the past year. Each is a 0/1 variable, and each model is a weighted linear probability model, so the coefficient on lep is a difference in percentage points.
local ys nousual delay erv uninsured needhelp spd
drop if missing(needhelp)
tab educ, gen(ed)
tab pov, gen(pv)
local x "age female ed2 ed3 ed4 pv2 pv3 pv4"
foreach y of local ys {
regress `y' lep `x' [pw=rakedw0]
}The rows for lep from the six regressions:
nousual lep | .0548357 .0196216 2.79 0.005 .0163761 .0932953
delay lep | -.0208476 .0168985 -1.23 0.217 -.0539696 .0122743
erv lep | -.0369416 .0194525 -1.90 0.058 -.0750697 .0011865
uninsured lep | .0743223 .0173396 4.29 0.000 .0403358 .1083089
needhelp lep | -.1037591 .0166768 -6.22 0.000 -.1364467 -.0710714
spd lep | -.0701126 .0138896 -5.05 0.000 -.0973371 -.0428882Holding age, gender, education, and income constant, adults with limited English proficiency were 5.5 percentage points more likely to have no usual source of care and 7.4 points more likely to be uninsured, and they were 10.4 points less likely to report needing help for a mental health or substance use problem and 7.0 points less likely to meet the screen for serious psychological distress. The unadjusted p-values for delayed care (0.217) and emergency room visits (0.058) are above 0.05.
Westfall-Young adjusted p-values
wyoung takes the list of outcomes, a model template in which OUTCOMEVAR stands for each outcome in turn, and the variable whose coefficients form the family of hypotheses.
wyoung `ys', cmd(regress OUTCOMEVAR lep `x' [pw=rakedw0]) ///
familyp(lep) bootstraps(500) seed(20260920)| familyp | coef | stderr | p | pwyoung | pbonf | psidak |
nousual | .05483569 | .01962163 | .0051995 | .026 | .01559851 | .01551755 |
delay | -.02084764 | .01689846 | .21732617 | .266 | .21732617 | .21732617 |
erv | -.03694159 | .01945252 | .05756761 | .094 | .11513522 | .11182119 |
uninsured | .07432235 | .01733956 | .00001824 | 0 | .00007295 | .00007295 |
needhelp | -.10375906 | .01667684 | 4.995e-10 | 0 | 2.997e-09 | 2.997e-09 |
spd | -.07011263 | .01388962 | 4.500e-07 | 0 | 2.250e-06 | 2.250e-06 |The output, shortened here to one row per outcome, gives the unadjusted p-value (p), the Westfall-Young adjusted p-value (pwyoung), and Bonferroni-Holm (pbonf) and Sidak-Holm (psidak) adjusted p-values for comparison. A pwyoung of 0 means none of the 500 bootstrap replications produced a test statistic as large as the observed one; report it as p < 0.002 rather than 0. Four of the six associations remain significant after adjustment. The result for no usual source of care moves from 0.005 to 0.026. The emergency room result, which was close to 0.05 before adjustment, is 0.094 after it.
Romano-Wolf adjusted p-values
rwolf2 takes one model per hypothesis and the independent variable to correct in each model, separated by commas.
rwolf2 (regress nousual lep `x' [pw=rakedw0]) (regress delay lep `x' [pw=rakedw0]) ///
(regress erv lep `x' [pw=rakedw0]) (regress uninsured lep `x' [pw=rakedw0]) ///
(regress needhelp lep `x' [pw=rakedw0]) (regress spd lep `x' [pw=rakedw0]), ///
indepvars(lep, lep, lep, lep, lep, lep) reps(500) seed(20260920) holmRomano-Wolf step-down adjusted p-values
Number of resamples: 500
------------------------------------------------------------------------------
| Model Resample Romano-Wolf Holm
| p-value p-value p-value p-value
----------------+-------------------------------------------------------------
nousual |
lep | 0.0052 0.0120 0.0279 0.0359
delay |
lep | 0.2173 0.2675 0.2675 0.2675
erv |
lep | 0.0576 0.0379 0.0958 0.0758
uninsured |
lep | 0.0000 0.0020 0.0020 0.0080
needhelp |
lep | 0.0000 0.0020 0.0020 0.0120
spd |
lep | 0.0000 0.0020 0.0020 0.0100
------------------------------------------------------------------------------The Romano-Wolf p-values are close to the Westfall-Young ones, and the conclusions are the same. By default rwolf2 adds one to both the count and the number of replications, so the smallest possible p-value with 500 replications is 1/501, which rounds to 0.0020 for three outcomes. The Holm column is larger than the Romano-Wolf column for most outcomes, which is the gain in power from accounting for the correlation between outcomes.
Choices to report
Define the family before looking at the results. The correction depends on which hypotheses are grouped together, so the family should be set in the analysis plan, usually all primary outcomes for one exposure.
Set the number of replications. Use several thousand replications for the final analysis; 500 keeps this example quick. Set a seed so the adjusted p-values can be reproduced.
State what the resampling ignores. Both commands resample individual respondents. With weighted survey data, this treats the weight as fixed and ignores the CHIS replicate-weight design, so the adjusted p-values are approximate. Both commands have cluster() and strata() options for designs where those are known.
References
Clarke, D., Romano, J. P., & Wolf, M. (2020). The Romano–Wolf multiple-hypothesis correction in Stata. The Stata Journal, 20(4), 812–843. https://doi.org/10.1177/1536867X20976314
Jones, D., Molitor, D., & Reif, J. (2019). What do workplace wellness programs do? Evidence from the Illinois Workplace Wellness Study. The Quarterly Journal of Economics, 134(4), 1747–1791. https://doi.org/10.1093/qje/qjz023
Romano, J. P., & Wolf, M. (2005). Stepwise multiple testing as formalized data snooping. Econometrica, 73(4), 1237–1282. https://doi.org/10.1111/j.1468-0262.2005.00615.x
