*! v.1.0.0 Matteo Bottai 13mar2020
cap program drop emagnification
program emagnification, rclass
version 14
syntax namelist , *
if strpos("proportion","`namelist'")==1 {
syntax namelist , p0(numlist >0 <1) or(numlist >0) ///
	n0(numlist integer >0) n1(numlist integer >0) ///
	[ pctile(numlist integer >0 <100) ifactor(numlist integer >0 <100) nsim(numlist max=1 integer >0) ///
	level(numlist max=1 >0 <1) onesided EXACT seed(string) log clean format(string) ]
if "`nsim'"=="" {
	local nsim 10
	}
if "`seed'"=="" {
	local seed = c(rngstate)
	}
if "`level'"=="" {
	local level .05
	}
if "`pctile'"=="" {
	local pctile = "10 50 90"
	}
if "`format'"=="" {
	local format %4.3f
	}
if "`exact'"=="" {
	local test "chi2"
	local retp = cond("`onesided'"=="","r(p)","r(p)/2")
	}
	else {
	local test "exact"
	local retp = cond("`onesided'"=="","r(p_exact)","r(p1_exact)")
	}
preserve
qui drop _all
set seed `seed'
qui set obs `=(`nsim')*(`: word count `p0'')*(`: word count `or'')*(`: word count `n0'')*(`: word count `n1'')'
qui gen p0 = .
qui gen p1 = .
qui gen true_or = .
qui gen n0 = .
qui gen n1 = .
qui gen sim = .
qui gen or = .
qui gen p = .
local c = 0
local s = 0
qui foreach p_0 of local p0 {
foreach o of local or {
local p_1 = (`o'*`p_0'/(1-`p_0'))/(1+(`o'*`p_0'/(1-`p_0')))
foreach n_0 of local n0 {
foreach n_1 of local n1 {
local ++s
if "`log'"!="" {
	noi di _new as txt "Scenario " as res `s' as txt ": p0 = " as res `p_0' as txt ", or = " as res `o' as txt ", n0 = " as res `n_0' as txt ", n1 = " as res `n_1'
	noi di as txt "Completed: " _cont
	}
cap set obs `=`n_0'+`n_1''
gen x = _n>`n_0' in 1/`=`n_0'+`n_1''
forv i = 1/`nsim' {
if "`log'"!="" & mod(`i',round(`nsim'/10))==0 {
	noi di as txt int(`i'/round(`nsim'/10))*10 "% " _cont
	}
local ++c
replace p0 = `p_0' in `c'
replace p1 = `p_1' in `c'
replace true_or = `o' in `c'
replace n0 = `n_0' in `c'
replace n1 = `n_1' in `c'
replace sim = `i' in `c'
gen y = (runiform()<(cond(x==0,`p_0',`p_1')))
tab y x, matcell(f) `test'
replace or = f[1,1]*f[2,2]/f[1,2]/f[2,1] in `c'
replace p = `retp' in `c'
drop y
}
drop x
}
}
}
}
// qui replace or = c(maxfloat) if sim!=. & p!=. & or==.
qui drop if sim==. | p==.
qui count
if r(N)==0 {
	noi di as err "No valid simulations"
	exit
	}
qui gen if_or = or/true_or
qui egen group = group(p0 true_or n0 n1)
qui egen valid = count(p), by(group)
qui egen power = mean((p<`level') & ((or>1)==(true_or>1) | "`onesided'"=="")), by(group)
qui keep if (p<`level') & ((or>1)==(true_or>1) | "`onesided'"=="")
qui count
if r(N)==0 {
	noi di as err "No significant p-value; increase the number of simulations"
	exit
	}
qui foreach p of local pctile {
	egen p`p' = pctile(or), p(`p') by(group)
	format `format' p`p'
	}
qui foreach p of local ifactor {
	egen if_p`p' = pctile(if_or), p(`p') by(group)
	format `format' if_p`p'
	}
qui egen tag = tag(group)
qui keep if tag==1
drop or p sim group tag if_or
di _newline(2) as txt "The tests are " as res cond("`onesided'"=="", "two-sided", "one-sided") as txt " with level = " as res `level'
list , noobs sepby(p0) `clean'
}
else if strpos("rate","`namelist'")==1 {
syntax namelist , r0(numlist >0) rr(numlist >0) ///
	n0(numlist integer >0) n1(numlist integer >0) ///
	[ pctile(numlist integer >0 <100) ifactor(numlist integer >0 <100) nsim(numlist max=1 integer >0) ///
	level(numlist max=1 >0 <1) onesided EXACT seed(string) log clean format(string) ]
if "`nsim'"=="" {
	local nsim 10
	}
if "`seed'"=="" {
	local seed = c(rngstate)
	}
if "`level'"=="" {
	local level = .05
	}
if "`pctile'"=="" {
	local pctile = "10 50 90"
	}
if "`format'"=="" {
	local format %4.3f
	}
preserve
qui drop _all
set seed `seed'
qui set obs `=(`nsim')*(`: word count `r0'')*(`: word count `rr'')*(`: word count `n0'')*(`: word count `n1'')'
qui gen r0 = .
qui gen r1 = .
qui gen true_rr = .
qui gen n0 = .
qui gen n1 = .
qui gen sim = .
qui gen rr = .
qui gen p = .
local c = 0
local s = 0
qui foreach r_0 of local r0 {
foreach r of local rr {
local r_1 = `r'*`r_0'
foreach n_0 of local n0 {
foreach n_1 of local n1 {
local ++s
if "`log'"!="" {
	noi di _new as txt "Scenario " as res `s' as txt ": r0 = " as res `r_0' as txt ", rr = " as res `r' as txt ", n0 = " as res `n_0' as txt ", n1 = " as res `n_1'
	noi di as txt "Completed: " _cont
	}
cap set obs `=`n_0'+`n_1''
gen x = _n>`n_0' in 1/`=`n_0'+`n_1''
forv i = 1/`nsim' {
if "`log'"!="" & mod(`i',round(`nsim'/10))==0 {
	noi di as txt int(`i'/round(`nsim'/10))*10 "% " _cont
	}
local ++c
qui replace r0 = `r_0' in `c'
qui replace r1 = `r_1' in `c'
qui replace true_rr = `r' in `c'
qui replace n0 = `n_0' in `c'
qui replace n1 = `n_1' in `c'
qui replace sim = `i' in `c'
gen y = rpoisson(cond(x==0,`r_0',`r_1'))
cap poisson y x
if _rc==0 {
	replace rr = exp(_b[x]) in `c'
	if "`exact'"=="" {
		replace p = normal(-abs(_b[x]/_se[x]))*2/(1+("`onesided'"!="")) in `c'
		}
		else {
		cap expoisson y x
		if _rc==0 {
			mat p = e(p_sufficient)
			replace p = p[1,1]/(1+("`onesided'"!="")) in `c'
			}
		}
	}
drop y
}
drop x
}
}
}
}
qui drop if sim==. | p==.
if r(N)==0 {
	noi di as err "No valid simulations"
	exit
	}
qui gen if_rr = rr/true_rr
qui egen group = group(r0 true_rr n0 n1)
qui egen valid = count(p), by(group)
qui egen power = mean((p<`level') & (sign(rr-1)==sign(true_rr-1))), by(group)
qui keep if (p<`level') & (sign(rr-1)==sign(true_rr-1))
qui count
if r(N)==0 {
	noi di as err "No significant p-value; increase the number of simulations"
	exit
	}
qui foreach p of local pctile {
	egen p`p' = pctile(rr), p(`p') by(group)
	format `format' p`p'
	}
qui foreach p of local ifactor {
	egen if_p`p' = pctile(if_rr), p(`p') by(group)
	format `format' if_p`p'
	}
qui egen tag = tag(group)
qui keep if tag==1
drop rr p sim group tag if_rr
di _newline(2) as txt "The tests are " as res cond("`onesided'"=="", "two-sided", "one-sided") as txt " with level = " as res `level'
list , noobs sepby(r0) `clean'
}
else {
	di as err "Only 'proportion' and 'rate' are allowed"
	exit
	}
tempname table
mkmat _all, matrix(`table')
return matrix table `table'
return scalar level = `level'
return local seed `seed'
return local cmdline "emagnification `0'"
restore
end