Two stage power calculation

two-stage-power.do

// One run of simulation

cap prog drop powerfinder
prog def powerfinder , rclass
args alloc

  // Set up 1000 hypotheses
  clear
  set obs 1000
  gen id = _n
  gen effect = rnormal() * (runiform() < 0.05) // 5% have an effect

  // Evaluate a study for each hypothesis
  gen n = 100000 * `alloc' / _N
  expand n

  gen realization = effect + rnormal()
  collapse (mean) realization effect n (semean) se = realization , by(id)

  // T-statistics and significance
  gen t = realization / se
  gen sig = abs(t) > 2

  gen false = effect == 0 & sig == 1
  sort false

  su sig if effect != 0
    return scalar s1power = r(mean)
  su false if sig != 0
    return scalar s1false = r(mean)

  // Second stage
  keep if sig
  keep id effect
  gen n = round( 100000 * (1-`alloc') / (_N) , 1)

  // Re-study
  expand n

  gen realization = effect + rnormal()
  collapse (mean) realization effect n (semean) se = realization , by(id)

  // T-statistics and significance
  gen t = realization / se
  gen sig = abs(t) > 2

  gen false = effect == 0 & sig == 1
  sort false

  su sig if effect != 0
    return scalar s2power = r(mean)
  su false if sig != 0
    return scalar s2false = r(mean)


end

// Test

powerfinder 0.5
return list

// Sim

clear
tempfile results
save `results' , emptyok

forvalues i = 0.05(0.05)0.95 {
  simulate ///
    s1power = r(s1power) s1false = r(s1false) s2power = r(s2power) s2false = r(s2false) ///
  , reps(1000) : powerfinder `i'
  gen alloc = `i'
  append using `results'
  save `results' , replace

}

save "/users/bbdaniels/desktop/temp.dta" , replace
use "/users/bbdaniels/desktop/temp.dta" , clear

collapse (mean) s1power s1false s2power s2false , by(alloc) fast

tw (lowess s1power alloc) (lowess s1false alloc) (lowess s2power alloc) (lowess s2false alloc) ///
  , legend(on order(1 "Stage 1 Power" 2 "Stage 1 False Pos" 3 "Stage 2 Power" 4 "Stage 2 False Pos")) ///
  xtit("Allocation to Stage 1 (100k Obs, 1k Hypotheses)") xscale(log)

// End