// owen ozier
// 2012.05
// quick demo of checking a simple power calculation in stata


// PART 1.1: PARAMETERS OF SIMULATION
local loopmax=1000
set seed 123

// PART 1.2: PARAMETERS OF CALCULATION
local power=0.8
local effectsize=0.1
local sderror=0.4



clear


// PART 2.1: USING SAMPSI TO GET THE ANSWER
sampsi 0 `effectsize', sd1(`sderror') power(`power')
di "Sampsi says we need: " r(N_1) " obs per treatment."
local N1=r(N_1)
local doubleN1=2*r(N_1)

// PART 2.2: USING A VERSION OF THE FORMULA TO GET THE ANSWER
local byhandroughly = ceil(2 * ((invnorm(`power')+invnorm(0.975)) / (`effectsize'/`sderror') )^2)
di "The simple version of the formula (normal rather than T) says we need: `byhandroughly' obs per treatment."

// PART 2.3: USING A SIMULATION TO CONFIRM THE ANSWER
// a bit of housekeeping first:
set obs `doubleN1'
g n=_n
g t=_n>`N1'
local rejectcount=0
local shouldbeabout=`power'*`loopmax'
// then the simulation loop
forvalues i=1/`loopmax' {
 quietly {
  cap drop e
  g e=invnorm(uniform())*`sderror'
  cap drop y
  g y=e+t*`effectsize'
  reg y t
  if 2*ttail(e(df_r),abs(_b[t]/_se[t]))<0.05 {
    local rejectcount=`rejectcount'+1
  }
 } // quietly
} // forvalues
di "Rejected null: `rejectcount' out of `loopmax' with effect `effectsize' and sample size `doubleN1'"
di "(Should be about `shouldbeabout' rejections, if calculations were correct.)"