// ECON 626 FALL 2019
// L9. Multiple Hypothesis Testing
// A1. Implementing Anderson's Benjamini-Hochberg q-values
clear all
set seed 12345
cd "E:\Dropbox\econ-626-2019\lectures\L9 Multiple Testing\activities"
// 1. Generate a list of k=20 p-values, including several below 0.05 that are
// quite close together and many that would normally be rejected (absent
// adjustment). How many hypotheses would be rejected?
set obs 20
gen pval = runiform()/10
sort pval
gen p_reject = pval<=0.05
// 1a. Calculate Bonferroni-adjusted p-values - how many would be rejected now?
count // number of hypotheses being tested
local k = r(N)
gen bval = min(pval*`k',1)
gen b_reject = bval<=0.05
// 1b. Calculate Benjamini-Hochberg q-values as follows:
// - sort p-values and calculate ranks (smallest to largest)
// - calculate preliminary adjusted q-values by multiplying p by rank/k
// - Whenever the q-value for hypothesis j is above that of hypothesis j+1 (etc.)
// adjust accordingly
count // number of hypotheses being tested
local k = r(N)
sort pval
gen rank = _n
gen temp_q = pval * (`k'/rank)
gen qval = temp_q // will be replaced as needed
local j = `k' - 1
sum temp_q in `k'
local thisqval = r(mean)
forvalues i = 1/`j' {
local temprank = `k' - `i'
sum temp_q in `temprank'
local tempqval = r(mean)
if `tempqval' > `thisqval' {
replace qval = `thisqval' in `temprank'
}
if `tempqval' < `thisqval' {
local thisqval = `tempqval'
}
}
gen q_reject = qval<=0.05
// 1c. How many hypotheses are rejected under each of the three approaches
// (no adjustments, Bonferroni, and Benjamini-Hochberg)?
// 1d. Compare the q-values you calculated to those you obtain using
// Michael Anderson's do file (anderson_qvalues.do) - are they similar?