*! (SJ3-1: st0032)

capture program drop ggipower
program define ggipower

quietly {

clear

syntax using/

global using=`"`using'"'

version 7
set more off
capture set mem 16m
use `using'
count
local scenar=r(N)

cap which vars2
if _rc~=0 {
  di in red "YOU MUST INSTALL vars2.ado to run this"
  di "This function calulates the risk factor frequencies"
  exit
}


count if pg1<=0 | pg1>=1 | pg2<=0 | pg2>=1 | pd<=0 | pd>1 | rrg1<=0 | rrg2<=0| rrint<=0 |rrsibm<=0 | alpha_1 <=0 |alpha_1 >=1 |power<5 |power>99
if r(N)~=0 {
noisily {
display "one or more of the parameter values entered is outside permitted ranges"
display "please check and change your initial dataset, and re-save it."
display "type h geipower for a description of the input variables needed"
}
exit
}

count if inh1~="d" & inh1~="D" & inh1~="r" & inh1~="R"
if r(N)~=0 {
noisily {
display "inh1 can only have the values D or R (dominant or recessive)"
display "please check and change your initial dataset, and re-save it."
display "type h geipower for a description of the input variables needed"
}
exit
}

count if inh2~="d" & inh2~="D" & inh2~="r" & inh2~="R"
if r(N)~=0 {
noisily {
display "inh2 can only have the values D or R (dominant or recessive)"
display "please check and change your initial dataset, and re-save it."
display "type h geipower for a description of the input variables needed"
}
exit
}

count if ssize<30
if r(N)~=0 {
noisily {
display "you have entered one or more values for sample size that are less than 30"
display "this program is based on a large sample approximation and so results may not be accurate"
display
}
}

drop _all

vars2

local loop3=0
while `loop3'~=`scenar' {
local loop3=`loop3'+1

if mod(`loop3', 10)==0 {
noisily {
display "calculating observation `loop3' out of `scenar'"
}
}

local error1=pop_error in `loop3'
local error2=sib_error in `loop3'
use `using'
local ssize=ssize in `loop3'
local power=power in `loop3'
local power=`power'/100
local alpha_1=alpha_1 in `loop3'
local alpha_1=invchi2(1, (1-`alpha_1'))
if `error1'==0 & `error2'==0 {
local	a1=dom_a in `loop3'
local	b1=dom_b in `loop3'	
local	c1=dom_c in `loop3'	
local	d1=dom_d in `loop3'	
local	e1=dom_e in `loop3'	
local	f1=dom_f in `loop3'	
local	g1=dom_g in `loop3'	
local	h1=dom_h in `loop3'	
local	i1=dom_i in `loop3'	
local	j1=dom_j in `loop3'	
local	k1=dom_k in `loop3'	
local	l1=dom_l in `loop3'	
local	m1=dom_m in `loop3'	
local	n1=dom_n in `loop3'	
local	o1=dom_o in `loop3'	
local p1=dom_p in `loop3'
drop _all
set obs 16
gen id=_n
gen g0=0
replace g0=1 in 3/4
replace g0=1 in 7/8
replace g0=1 in 11/12
replace g0=1 in 15/16
gen e0=0
replace e0=1 in 2
replace e0=1 in 4
replace e0=1 in 6
replace e0=1 in 8
replace e0=1 in 10
replace e0=1 in 12
replace e0=1 in 14
replace e0=1 in 16
gen gei0=g0*e0
gen g1=0
replace g1=1 in 9/16
gen e1=0
replace e1=1 in 5/8
replace e1=1 in 13/16
gen gei1=g1*e1
gen weight=.
replace weight=`a1' in 1
replace weight=`b1' in 2
replace weight=`c1' in 3
replace weight=`d1' in 4
replace weight=`e1' in 5
replace weight=`f1' in 6
replace weight=`g1' in 7
replace weight=`h1' in 8
replace weight=`i1' in 9
replace weight=`j1' in 10
replace weight=`k1' in 11
replace weight=`l1' in 12
replace weight=`m1' in 13
replace weight=`n1' in 14
replace weight=`o1' in 15
replace weight=`p1' in 16
replace weight=round(weight*`ssize'*10000,1)
reshape long g e gei, i(id) j(case)
clogit case g e [w=weight], group(id) or
local temp1=e(ll)
clogit case g e gei [w=weight], group(id) or
matrix A=e(b)
local temp2=e(ll)
local temp3=2*(`temp2'-`temp1')/10000
local m_sib= 1-nchi(1, `temp3', `alpha_1')
local or_m_sib=exp(el(A, 1, 3))
local ss_m_sib=npnchi2(1, `alpha_1', (1-`power'))
local ss_m_sib=`ss_m_sib'*(`ssize')/`temp3'
drop _all
set obs 4
gen g=0
gen e=0
replace g=1 in 3/4
replace e=1 in 2
replace e=1 in 4
gen gei=g*e
gen aff=.
gen tot=.
local aff1=(`a1'+`b1'+`c1'+`d1')
local tot1=(`a1'+`b1'+`c1'+`d1')+(`a1'+`e1'+`i1'+`m1')
local aff2=(`e1'+`f1'+`g1'+`h1')
local tot2=(`e1'+`f1'+`g1'+`h1')+(`b1'+`f1'+`j1'+`n1')
local aff3=(`i1'+`j1'+`k1'+`l1')
local tot3=(`i1'+`j1'+`k1'+`l1')+(`c1'+`g1'+`k1'+`o1')
local aff4=(`m1'+`n1'+`o1'+`p1')
local tot4=(`m1'+`n1'+`o1'+`p1'+`d1'+`h1'+`l1'+`p1')
replace aff=`aff1' in 1
replace tot=`tot1' in 1
replace aff=`aff2' in 2
replace tot=`tot2' in 2
replace aff=`aff3' in 3
replace tot=`tot3' in 3
replace aff=`aff4' in 4
replace tot=`tot4' in 4
replace aff=round(aff*`ssize'*10000, 1)
replace tot=round(tot*`ssize'*10000, 1)
blogit aff tot g e, or
local temp5=e(ll)
blogit aff tot g e gei, or
local temp6=e(ll)
local temp7=2*(`temp6'-`temp5')/10000
local u_sib= 1-nchi(1, `temp7', `alpha_1')
local ss_u_sib=npnchi2(1, `alpha_1', (1-`power'))
local ss_u_sib=`ss_u_sib'*(`ssize')/`temp7'
matrix A=e(b)
local or_u_sib=exp(el(A, 1, 3))

drop _all
use `using'
local aff1=p_g0e0_d in `loop3'
local tot1=p_g0e0_not_d in `loop3'
local tot1=`tot1'+`aff1'
local aff2=p_g0e1_d in `loop3'
local tot2=p_g0e1_not_d in `loop3'
local tot2=`tot2'+`aff2'
local aff3=p_g1e0_d in `loop3'
local tot3=p_g1e0_not_d in `loop3'
local tot3=`tot3'+`aff3'
local aff4=p_g1e1_d in `loop3'
local tot4=p_g1e1_not_d in `loop3'
local tot4=`tot4'+`aff4'
drop _all
set obs 4
gen g=0
gen e=0
replace g=1 in 3/4
replace e=1 in 2
replace e=1 in 4
gen gei=g*e
gen aff=.
gen tot=.
replace aff=`aff1' in 1
replace tot=`tot1' in 1
replace aff=`aff2' in 2
replace tot=`tot2' in 2
replace aff=`aff3' in 3
replace tot=`tot3' in 3
replace aff=`aff4' in 4
replace tot=`tot4' in 4
replace aff=round(aff*`ssize'*10000, 1)
replace tot=round(tot*`ssize'*10000, 1)
blogit aff tot g e, or
local temp5=e(ll)
blogit aff tot g e gei, or
local temp6=e(ll)
local temp7=2*(`temp6'-`temp5')/10000
local cc_u= 1-nchi(1, `temp7', `alpha_1')
local ss_cc_u=npnchi2(1, `alpha_1', (1-`power'))
local ss_cc_u=`ss_cc_u'*(`ssize')/`temp7'
matrix A=e(b)
local or_cc_u=exp(el(A, 1, 3))


drop _all
use `using'
local aff1=dom_p_g0e0_soc_d in `loop3'
local tot1=p_g0e0_not_d in `loop3'
local tot1=`tot1'+`aff1'
local aff2=dom_p_g0e1_soc_d in `loop3'
local tot2=p_g0e1_not_d in `loop3'
local tot2=`tot2'+`aff2'
local aff3=dom_p_g1e0_soc_d in `loop3'
local tot3=p_g1e0_not_d in `loop3'
local tot3=`tot3'+`aff3'
local aff4=dom_p_g1e1_soc_d in `loop3'
local tot4=p_g1e1_not_d in `loop3'
local tot4=`tot4'+`aff4'
drop _all
set obs 4
gen g=0
gen e=0
replace g=1 in 3/4
replace e=1 in 2
replace e=1 in 4
gen gei=g*e
gen aff=.
gen tot=.
replace aff=`aff1' in 1
replace tot=`tot1' in 1
replace aff=`aff2' in 2
replace tot=`tot2' in 2
replace aff=`aff3' in 3
replace tot=`tot3' in 3
replace aff=`aff4' in 4
replace tot=`tot4' in 4
replace aff=round(aff*`ssize'*10000, 1)	
replace tot=round(tot*`ssize'*10000, 1)
blogit aff tot g e, or
local temp5=e(ll)
blogit aff tot g e gei, or
local temp6=e(ll)
local temp7=2*(`temp6'-`temp5')/10000
local asp= 1-nchi(1, `temp7', `alpha_1')
local ss_asp=npnchi2(1, `alpha_1', (1-`power'))
local ss_asp=`ss_asp'*(`ssize')/`temp7'
matrix A=e(b)
local or_asp=exp(el(A, 1, 3))


drop _all
use `using'
local aff1=dom_p_g0e0_soc_d in `loop3'
local tot1=dom_p_g0e0_soc_nd in `loop3'
local tot1=`tot1'+`aff1'
local aff2=dom_p_g0e1_soc_d in `loop3'
local tot2=dom_p_g0e1_soc_nd in `loop3'
local tot2=`tot2'+`aff2'
local aff3=dom_p_g1e0_soc_d in `loop3'
local tot3=dom_p_g1e0_soc_nd in `loop3'
local tot3=`tot3'+`aff3'
local aff4=dom_p_g1e1_soc_d in `loop3'
local tot4=dom_p_g1e1_soc_nd in `loop3'
local tot4=`tot4'+`aff4'
drop _all
set obs 4
gen g=0
gen e=0
replace g=1 in 3/4
replace e=1 in 2
replace e=1 in 4
gen gei=g*e
gen aff=.
gen tot=.
replace aff=`aff1' in 1
replace tot=`tot1' in 1
replace aff=`aff2' in 2
replace tot=`tot2' in 2
replace aff=`aff3' in 3
replace tot=`tot3' in 3
replace aff=`aff4' in 4
replace tot=`tot4' in 4
replace aff=round(aff*`ssize'*10000, 1)
replace tot=round(tot*`ssize'*10000, 1)
blogit aff tot g e, or
local temp5=e(ll)
blogit aff tot g e gei, or
local temp6=e(ll)
local temp7=2*(`temp6'-`temp5')/10000
local coh= 1-nchi(1, `temp7', `alpha_1')
local ss_coh=npnchi2(1, `alpha_1', (1-`power'))
local ss_coh=`ss_coh'*(`ssize')/`temp7'
matrix A=e(b)
local or_coh=exp(el(A, 1, 3))


drop _all
use `using'


display `loop3'
replace m_sib=`m_sib' in `loop3'
replace u_sib=`u_sib' in `loop3'
replace cc_u=`cc_u' in `loop3'
replace asp=`asp' in `loop3'
replace coh=`coh' in `loop3'
replace or_m_sib=`or_m_sib' in `loop3'
replace or_u_sib=`or_u_sib' in `loop3'
replace or_cc_u=`or_cc_u' in `loop3'
replace or_asp=`or_asp' in `loop3'
replace or_coh=`or_coh' in `loop3'
replace ss_m_sib=`ss_m_sib' in `loop3'
replace ss_u_sib=`ss_u_sib' in `loop3'
replace ss_cc_u=`ss_cc_u' in `loop3'
replace ss_asp=`ss_asp' in `loop3'
replace ss_coh=`ss_coh' in `loop3'
save `using', replace
}

else if `error2'==float(1) & `error1'==float(0) {
drop _all
use `using'
local aff1=p_g0e0_d in `loop3'
local tot1=p_g0e0_not_d in `loop3'
local tot1=`tot1'+`aff1'
local aff2=p_g0e1_d in `loop3'
local tot2=p_g0e1_not_d in `loop3'
local tot2=`tot2'+`aff2'
local aff3=p_g1e0_d in `loop3'
local tot3=p_g1e0_not_d in `loop3'
local tot3=`tot3'+`aff3'
local aff4=p_g1e1_d in `loop3'
local tot4=p_g1e1_not_d in `loop3'
local tot4=`tot4'+`aff4'
drop _all
set obs 4
gen g=0
gen e=0
replace g=1 in 3/4
replace e=1 in 2
replace e=1 in 4
gen gei=g*e
gen aff=.
gen tot=.
replace aff=`aff1' in 1
replace tot=`tot1' in 1
replace aff=`aff2' in 2
replace tot=`tot2' in 2
replace aff=`aff3' in 3
replace tot=`tot3' in 3
replace aff=`aff4' in 4
replace tot=`tot4' in 4
replace aff=round(aff*`ssize'*100, 1)
replace tot=round(tot*`ssize'*100, 1)
blogit aff tot g e, or
local temp5=e(ll)
blogit aff tot g e gei, or
local temp6=e(ll)
local temp7=2*(`temp6'-`temp5')/100
local cc_u= 1-nchi(1, `temp7', `alpha_1')
local ss_cc_u=npnchi2(1, `alpha_1', (1-`power'))
local ss_cc_u=`ss_cc_u'*(`ssize')/`temp7'
matrix A=e(b)
local or_cc_u=exp(el(A, 1, 3))
drop _all
use `using'
replace cc_u=`cc_u' in `loop3'
replace or_cc_u=`or_cc_u' in `loop3'
replace ss_cc_u=`ss_cc_u' in `loop3'
save `using', replace
}	
else {
}

}


drop _all
use `using'



replace m_sib=round(m_sib*100, 0.01)
replace u_sib=round(u_sib*100, 0.01)
replace cc_u=round(cc_u*100, 0.01)
replace asp=round(asp*100, 0.01)
replace coh=round(coh*100, 0.01)
replace or_m_sib=round(or_m_sib, 0.01)
replace or_u_sib=round(or_u_sib, 0.01)
replace or_cc_u=round(or_cc_u, 0.01)
replace or_asp=round(or_asp, 0.01)
replace or_coh=round(or_coh, 0.01)
replace ss_m_sib=round(ss_m_sib, 1)
replace ss_u_sib=round(ss_u_sib, 1)
replace ss_cc_u=round(ss_cc_u, 1)
replace ss_asp=round(ss_asp, 1)
replace ss_coh=round(ss_coh, 1)
save `using', replace

keep pg pe inh1 inh2 rrg0e1 rrg1e0 rrint orge pd rrsibm m_sib u_sib  cc_u asp coh or_m_sib or_u_sib  or_cc_u or_asp or_coh ssize power ss_m_sib  ss_u_sib  ss_cc_u ss_asp  ss_coh  alpha_1 pop_error sib_error

rename pg pg1
rename pe pg2
rename rrg1e0 rrg1
rename rrg0e1 rrg2

label var pg1 "frequency genotype 1"
label var pg2 "frequency genotype 2"
label var rrg1 "relative risk of disease genotype 1 only"
label var rrg2 "relative risk of disease genotype 2 only"
label var rrint "interaction relative risk"
label var orge "OR(GE)"
label var rrsibm "RR Sib(M)"
label var alpha_1 "significance level"
label var pd "disease frequency"
label var m_sib "Power (%) sib control, matched analysis"
label var u_sib "Power (%) sib control, unmatched analysis"
label var cc_u "Power (%) unmatched case-control study"
label var asp "Power (%) case from ASP"
label var coh "Power (%) 'cohort of sibs of cases'"
label var or_m_sib "Interaction odds ratio, sibling control (matched)"
label var or_u_sib "Interaction odds ratio, sibling control (unmatched)"
label var or_cc_u "Interaction odds ratio, population based"
label var or_asp "Interaction odds ratio, case from ASP"
label var or_coh "Interaction odds ratio, 'cohort of sibs of cases'"
label var ss_m_sib "sample size, matched sib design"
label var ss_u_sib "sample size, unmatched sib design"
label var ss_cc_u "sample size, population based case control study"
label var ss_asp "sample size, case from ASP"
label var ss_coh "sample size, cohort of sibs of cases"
save `using', replace

}

if `scenar'==1 {
foreach temp of varlist _all {
local `temp'=`temp'
}

local alpha_1=round(`alpha_1', 0.01)


di
di as txt "parameter file:" _col(30) as res "`using'"
di in gr "observations:" in ye _col(30) 1
di
di in gr "-------------------------------------------------------------------------------
di in gr "PARAMETERS
di in gr "-------------------------------------------------------------------------------
di
di in gr "Power to detect interaction (%), interaction odds ratios and required sample 
di in gr "sizes have been calculated from an exemplary dataset for the following
di in gr "population risk factor frequencies and effects:
di
display in gr "susceptibility genotype 1 frequency: (pg1)" in ye _col(60) pg1
display in gr "susceptibility genotype 1 inheritance: (inh1)" in ye _col(60) inh1
display in gr "susceptibility genotype 2 frequency: (pg2)" in ye _col(60) pg2
display in gr "susceptibility genotype 2 inheritance: (inh2)" in ye _col(60) inh2
display in gr "genetic risk factor 1 main effect: (rrg1)" in ye _col(60) rrg1
display in gr "genetic risk factor 1 main effect: (rrg2)" in ye _col(60) rrg2
display in gr "interaction relative risk: (rrint)" in ye _col(60) rrint
display in gr "disease prevalence: (pd)" in ye _col(60) pd
display in gr "relative risk of disease to sibs of cases: (rrsibm)" in ye _col(60) rrsibm
display
di in gr "-------------------------------------------------------------------------------
di in gr "POWER
di in gr "-------------------------------------------------------------------------------
di
di in gr "Power to detect an interaction for a sample size of `ssize' cases and `ssize'
di in gr "controls with a two-sided significance level (alpha=`alpha_1') 
di
display in gr "1) matched sib design, power (%):"  in ye _col(70)  m_sib
display in gr "2) unmatched sib design, power (%):"  in ye _col(70)  u_sib
display in gr "3) case from ASP, population control, power (%):"  in ye _col(70)  asp
display in gr "4) unmatched design, case and control both have affected sibling,"
di in gr _col(5) "power (%):"   in ye _col(70)   coh
display in gr "5) population based case-control study, power (%):"  in ye _col(70)  cc_u
display
di in gr "-------------------------------------------------------------------------------
di in gr "INTERACTION ODDS RATIOS
di in gr "-------------------------------------------------------------------------------
di
di in gr "The interaction odds ratio calculated from the exemplary dataset:
di
display in gr "1) matched sib design, interaction odds ratio:" in ye _col(70) or_m_sib
display in gr "2) unmatched sib design, interaction odds ratio:" in ye _col(70) or_u_sib
display in gr "3) case from ASP, population control, interaction odds ratio:" in ye _col(70) or_asp
display in gr "4) unmatched design, case and control both have affected sibling,"
di in gr _col(5) "interaction odds ratio:" in ye _col(70) or_coh
display in gr "5) population based case-control study:" in ye _col(70) or_cc_u
di
di in gr "-------------------------------------------------------------------------------
di in gr "REQUIRED SAMPLE SIZE
di in gr "-------------------------------------------------------------------------------
di
di in gr "for a power of `power'% and a two-sided significance level (alpha=`alpha_1') the 
di in gr "required number of cases (with equal numbers of controls)
di
display in gr "1) matched sib design:" in ye _col(70)  ss_m_sib
display in gr "2) unmatched sib design:" in ye _col(70)  ss_u_sib
display in gr "3) case from ASP, population control:" in ye _col(70)  ss_asp
display in gr "4) unmatched design, case and control both have affected sibling:" in ye _col(70)  ss_coh
display in gr "5) population based case-control study:" in ye _col(70)  ss_cc_u
display
di in gr "-------------------------------------------------------------------------------
di in gr "NOTES
di in gr "-------------------------------------------------------------------------------
di
di in gr "The power, required sample sizes and calculated interaction odds ratio for the 
di in gr "these designs have been saved into file - " in ye "`using'"
di in gr "Type d for details of the output variables
if pop_error==1 {
display
display in gr "All results are missing because the parameters specified in " in ye "`using'"
display in gr "lead to disease probabilities greater than one;"
display in gr "try specifying a lower disease frequency, or smaller effects."
}
else if pop_error==0 & sib_error==1 {
display
display "No results are given for the sibling designs because the values in "
display  in ye "`using'" in gr " lead to the probability of disease in siblings of cases"
di in gr "being greater than one."
display "try using a lower disease frequency or a smaller value for rrsibm"
}
drop sib_error pop_error

quietly {
save `using', replace
}

}
else {
di
di in gr "parameter file:" in ye _col(30) "`using'"
di in gr "observations:" in ye _col(30) `scenar'

display

di in gr "-------------------------------------------------------------------------------
di in gr "NOTES
di in gr "-------------------------------------------------------------------------------

di 
display "The power, required sample sizes and calculated interaction odds ratio
di in gr "for these designs have been saved into file - " in ye "`using'"
display in gr  "Type d for details of the output variables"
di
di in gr "If you want output written to the results window please just use one
di in gr "observation in the initial dataset


quietly {
sum sib_error
}
if r(mean)~=0 {
display
di in gr "One or more or the sets of parameters (observations) that you entered, resulted
di in gr "in disease probabilities being greater than one in some or all situations.  For
di in gr "these observations the results are posted as missing and variables pop_error 
di in gr "(population error) or sib_error (error in the sibling parameters will equal 1"
di in gr "rather than 0"
}
else {
drop pop_error sib_error


quietly {
save `using', replace
}

}
}
end