#delim cr
version 8.0
/*
 Run examples and create .eps figures
 for the Stata Journal paper on -powercal- entitled:
 "Power calculations for generalized linear models and more"
 (executable in Stata version 8).
 This program uses the -powercal-, -parmest-, -somersd- and -explist- packages
 (downloadable from SSC).
*/

set scheme sj

/*
 Plot lognormal tail ratios against lognormal CV
*/
clear
set obs 129
local cvmin 0.25
local cvmax 4
explist -2(1)2, base(2)
local xlabs "`r(explist)'"
gene double cv=exp( log(`cvmin') + (log(`cvmax')-log(`cvmin'))*(_n-1)/(_N-1) )
gene double sdlog=sqrt(log(cv*cv + 1))
lab var cv "Coefficient of variation"
lab var sdlog "SD of logs"
foreach Q of numlist 5 10 15 20 25 {
  scal q=`Q'/100
  gene double r_`Q'=exp(-2*sdlog*invnorm(q))
  lab var r_`Q' "`Q' percent tail ratio"
}
explist 0(1)8, base(2)
local ylabs "`r(explist)'"
line r_* cv, xscale(log) xlab(`xlabs', grid) yscale(log) ylab(`ylabs', grid) ///
  ytitle("Tail ratio (bottom of upper tail/top of lower tail)")
more
list

/*
 Example 1.1
*/
clear
scal cv=0.5
scal sdlog=sqrt(log(cv*cv + 1))
scal r20=exp(-2*sdlog*invnorm(0.2))
disp _n as text "Coefficient of variation: " as result cv ///
  _n as text "SD of logs: " as result sdlog ///
  _n as text "20% tail ratio: " as result r20
set obs 100
gene logratio=log(2)*(_n/_N)
lab var logratio "Log GM ratio"
gene gmratio=exp(logratio)
lab var gmratio "GM ratio"
powercal power, alpha(0.01) delta(logratio) sdinf(sdlog*sqrt(2)) ///
  nunit(50) tdf(98)
line power gmratio, ///
  yscale(range(0 1)) ///
  ylab(0(0.05)1, grid gmin gmax angle(0)) yline(0.8 0.9, lpattern(shortdash)) ///
  xscale(log range(1 2)) xlab(1(0.1)2, grid gmin gmax)
more
list
graph export figseq1.eps, replace

/*
 Example 1.2
*/
clear
scal cv=0.5
scal sdlog=sqrt(log(cv*cv + 1))
scal r20=exp(-2*sdlog*invnorm(0.2))
disp _n as text "Coefficient of variation: " as result cv ///
  _n as text "SD of logs: " as result sdlog ///
  _n as text "20% tail ratio: " as result r20
set obs 100
gene npergp=_n
lab var npergp "Number per group"
powercal logratio, power(0.9) alpha(0.01) sdinf(sdlog*sqrt(2)) ///
  nunit(npergp) tdf(2*(npergp-1))
gene hiratio=exp(logratio)
gene loratio=exp(-logratio)
lab var hiratio "Detectable GM ratio >1"
lab var loratio "Detectable GM ratio <1"
line hiratio loratio npergp if _n>=5, ///
  ylabel(, angle(0) grid gmin gmax) yline(1, lpattern(shortdash)) ///
  ytitle("Detectable GM ratio") ///
  xlab(0(10)100, grid gmin gmax)
more
list
graph export figseq2.eps, replace

/*
 Example 2
*/
clear
scal conprev=0.25
scal conodds=conprev/(1-conprev)
disp _n as text "Expected control prevalence: " as result conprev ///
  _n as text "Expected control odds: " as result conodds
set obs 101
gene logor=log(1.25)+(log(5)-log(1.25))*(_n-1)/(_N-1)
gene or=exp(logor)
gene caseodds=conodds*or
gene caseprev=caseodds/(1+caseodds)
gene sdinflor=sqrt( ///
  1/caseprev + 1/(1-caseprev) + (1/2)*( 1/conprev + 1/(1-conprev) ) ///
  )
lab var logor "Log odds ratio"
lab var or "Odds ratio"
lab var caseodds "Case exposure odds"
lab var caseprev "Case exposure prevalence"
lab var sdinflor "SD of influence for log OR"
desc
* Detectable OR by number of cases *
powercal ncases01, power(0.9) alpha(0.01) delta(logor) sdinf(sdinflor)
powercal ncases001, power(0.9) alpha(0.001) delta(logor) sdinf(sdinflor)
lab var ncases01 "Minimum cases (alpha=0.01)"
lab var ncases001 "Minimum cases (alpha=0.001)"
line or ncases01 if ncases01<=2000 || line or ncases001 if ncases001<=2000 , ///
  yscale(log range(1 5)) ylabel(1 1.5 2(1)5, angle(0) grid gmin gmax) ///
  xlab(, grid gmin gmax) xtitle("Minimum number of cases") ///
  legend(label(1 "Alpha=0.01") label(2 "Alpha=0.001"))
more
graph export figseq3.eps, replace
* Significance level by odds ratio *
powercal alpha50, power(0.50) delta(logor) sdinf(sdinflor) nunit(100)
powercal alpha90, power(0.90) delta(logor) sdinf(sdinflor) nunit(100)
line alpha50 or || line alpha90 or, ///
  yscale(log range(1e-9 1) reverse) ///
  ylab(1 0.05 1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-9, ///
    angle(0) format(%8.2g) grid gmin gmax) ///
  yline(0.05 0.01 0.001, lpattern(shortdash)) ///
  xscale(log) xlab(1 1.25 1.5 2(1)5, grid gmin gmax) ///
  legend(label(1 "Power=0.50") label(2 "Power=0.90"))
more
graph export figseq4.eps, replace
list

/*
 Example 3
*/
sysuse bpwide, clear
gene byte male=1-sex
lab var male "Male patient"
parmby "somersd male bp_before, tr(z) td", norestore ///
  escal(N) rename(es_1 N)
list N parm estimate stderr min* max* p, clean noobs
scal sdinf=stderr[1]*sqrt(N[1])
disp _n as text "SD of influence function for z-transformed Somers' D: " ///
  as result sdinf
drop _all
set obs 1000
gene int npat=_n
lab var npat "Number of patients"
foreach X in 05 01 001 0001 {
  powercal detz`X', power(0.9) alpha(0.`X') sdinf(sdinf) ///
    nunit(npat) tdf(npat-1)
  gene detd`X'=exp(2*detz`X')
  replace detd`X'=(detd`X'-1)/(detd`X'+1)
  lab var detz`X' "z-transformed Somers' D (P<=0.`X')"
  lab var detd`X' "Somers' D (P<=0.`X')"
}
line detd* npat, ///
 xlab(0(100)1000, grid gmin gmax) ///
 ylab(0(0.05)1, angle(0) grid gmin gmax) ///
 ytitle("Detectable Somers' D")
more
graph export figseq5.eps, replace
foreach X in 10 15 20 30 {
  scal z=0.5*log((1+0.`X')/(1-0.`X'))
  powercal alpha`X', power(0.9) delta(z) sdinf(sdinf) ///
    nunit(npat) tdf(npat-1)
  lab var alpha`X' "Alpha (Somers' D = 0.`X')"
  format alpha`X' %8.2g
}
line alpha* npat, ///
  xlab(0(100)1000, grid gmin gmax) ///
  yscale(log reverse) ///
  ylab(1 0.05 0.01 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-9 1e-10 1e-11, ///
    format(%8.2g) angle(0) grid gmin gmax) ///
  ytitle("Minimum alpha")
more
graph export figseq6.eps, replace
list npat detd* alpha*

exit