*! (SJ3-1: st0032)

capture program drop tdt_ggipower
program define tdt_ggipower

quietly {


syntax using/

set type double

clear

use `using'

keep pg1 pg2 rrg1 rrg2 rrint pd ssize power alpha_1 

gen dd_power=.
gen dd_ss=.
gen dd_or=.
gen dr_power=.
gen dr_ss=.
gen dr_or=.
gen rr_power=.
gen rr_ss=.
gen rr_or=.

gen pop_error=0
gen p_d_g0e0=pd/(1*(1-pg1)*(1-pg2)+rrg1*pg1*(1-pg2)+rrg2*pg2*(1-pg1)+rrg1*rrg2*rrint*pg1*pg2)
gen p_d_g0e1=p_d_g0e0*rrg1
gen p_d_g1e0=p_d_g0e0*rrg2
gen p_d_g1e1=p_d_g0e0*rrg1*rrg2*rrint

replace pop_error=1 if p_d_g0e0>1 | p_d_g0e1>1 | p_d_g1e0>1 | p_d_g1e1>1 | p_d_g0e0<0 | p_d_g0e1<0 | p_d_g1e0<0 | p_d_g1e1<0

drop p_d_g0e0 p_d_g0e1 p_d_g1e0 p_d_g1e1

save `using', replace

count if pg1<=0 | pg1>=1 | pg2<=0 | pg2>=1 | pd<=0 | pd>1 | rrg1<=0 | rrg2<=0| rrint<=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 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
}
}


count
local scenar=r(N)

local loop=0

while `loop'~=`scenar' {
local loop=`loop'+1

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

* population parameters for case parent exemplary dataset
local pg1=pg1 in `loop'
local pg2=pg2 in `loop'
local rrg1=rrg1 in `loop'
local rrg2=rrg2 in `loop'
local rrint=rrint in `loop'
local pd=pd in `loop'
local alpha_1=alpha_1 in `loop'
local sample_size=ssize in `loop'
local power=power in `loop'

drop _all

* assume independence of loci
* alternatively simply specify the four possible combination frequencies

tempfile gene2 gene1
local g00=(1-`pg1')*(1-`pg2')
local g01=(1-`pg1')*`pg2'
local g10=`pg1'*(1-`pg2')
local g11=`pg1'*`pg2'

local rr00=1
local rr01=`rrg2'
local rr10=`rrg1'
local rr11=`rrg1'*`rrg2'*`rrint'

local sumproduct=`g00'*`rr00'+`g01'*`rr01'+`g10'*`rr10'+`g11'*`rr11'

local p_g_00_d=`g00'*`rr00'/`sumproduct'
local p_g_01_d=`g01'*`rr01'/`sumproduct'
local p_g_10_d=`g10'*`rr10'/`sumproduct'
local p_g_11_d=`g11'*`rr11'/`sumproduct'

local af_g1_rec=`pg1'^0.5
local af_g1_dom=1-(1-`pg1')^0.5
local af_g2_rec=`pg2'^0.5
local af_g2_dom=1-(1-`pg2')^0.5

local p_00_rr=(1-`af_g1_rec')^2 * (1-`af_g2_rec')^2
local p_01_rr=(1-`af_g1_rec')^2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_02_rr=(1-`af_g1_rec')^2 * `af_g2_rec'^2
local p_10_rr=(1-`af_g1_rec') * `af_g1_rec' * 2 * (1-`af_g2_rec')^2
local p_11_rr=(1-`af_g1_rec') * `af_g1_rec' * 2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_12_rr=(1-`af_g1_rec') * `af_g1_rec' * 2 * `af_g2_rec'^2
local p_20_rr=`af_g1_rec'^2 * (1-`af_g2_rec')^2
local p_21_rr=`af_g1_rec'^2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_22_rr=`af_g1_rec'^2 * `af_g2_rec'^2

local p_00_dr=(1-`af_g1_dom')^2 * (1-`af_g2_rec')^2
local p_01_dr=(1-`af_g1_dom')^2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_02_dr=(1-`af_g1_dom')^2 * `af_g2_rec'^2
local p_10_dr=(1-`af_g1_dom') * `af_g1_dom' * 2 * (1-`af_g2_rec')^2
local p_11_dr=(1-`af_g1_dom') * `af_g1_dom' * 2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_12_dr=(1-`af_g1_dom') * `af_g1_dom' * 2 * `af_g2_rec'^2
local p_20_dr=`af_g1_dom'^2 * (1-`af_g2_rec')^2
local p_21_dr=`af_g1_dom'^2 * (1-`af_g2_rec') * `af_g2_rec' * 2
local p_22_dr=`af_g1_dom'^2 * `af_g2_rec'^2

local p_00_dd=(1-`af_g1_dom')^2 * (1-`af_g2_dom')^2
local p_01_dd=(1-`af_g1_dom')^2 * (1-`af_g2_dom') * `af_g2_dom' * 2
local p_02_dd=(1-`af_g1_dom')^2 * `af_g2_dom'^2
local p_10_dd=(1-`af_g1_dom') * `af_g1_dom' * 2 * (1-`af_g2_dom')^2
local p_11_dd=(1-`af_g1_dom') * `af_g1_dom' * 2 * (1-`af_g2_dom') * `af_g2_dom' * 2
local p_12_dd=(1-`af_g1_dom') * `af_g1_dom' * 2 * `af_g2_dom'^2
local p_20_dd=`af_g1_dom'^2 * (1-`af_g2_dom')^2
local p_21_dd=`af_g1_dom'^2 * (1-`af_g2_dom') * `af_g2_dom' * 2
local p_22_dd=`af_g1_dom'^2 * `af_g2_dom'^2

local rr_00_rr=`rr00'
local rr_01_rr=`rr00'
local rr_02_rr=`rr01'
local rr_10_rr=`rr00'
local rr_11_rr=`rr00'
local rr_12_rr=`rr01'
local rr_20_rr=`rr10'
local rr_21_rr=`rr10'
local rr_22_rr=`rr11'

local sumproduct_rr=`rr_00_rr' * `p_00_rr' + `rr_01_rr' * `p_01_rr' + `rr_02_rr' * `p_02_rr' + `rr_10_rr' * `p_10_rr' + `rr_11_rr' * `p_11_rr' + `rr_12_rr' * `p_12_rr' + `rr_20_rr' * `p_20_rr' + `rr_21_rr' * `p_21_rr' + `rr_22_rr' * `p_22_rr'

local rr_00_dr=`rr00'
local rr_01_dr=`rr00'
local rr_02_dr=`rr01'
local rr_10_dr=`rr10'
local rr_11_dr=`rr10'
local rr_12_dr=`rr11'
local rr_20_dr=`rr10'
local rr_21_dr=`rr10'
local rr_22_dr=`rr11'

local sumproduct_dr=`rr_00_dr' * `p_00_dr' + `rr_01_dr' * `p_01_dr' + `rr_02_dr' * `p_02_dr' + `rr_10_dr' * `p_10_dr' + `rr_11_dr' * `p_11_dr' + `rr_12_dr' * `p_12_dr' + `rr_20_dr' * `p_20_dr' + `rr_21_dr' * `p_21_dr' + `rr_22_dr' * `p_22_dr'

local rr_00_dd=`rr00'
local rr_01_dd=`rr01'
local rr_02_dd=`rr01'
local rr_10_dd=`rr10'
local rr_11_dd=`rr11'
local rr_12_dd=`rr11'
local rr_20_dd=`rr10'
local rr_21_dd=`rr11'
local rr_22_dd=`rr11'

local sumproduct_dd=`rr_00_dd' * `p_00_dd' + `rr_01_dd' * `p_01_dd' + `rr_02_dd' * `p_02_dd' + `rr_10_dd' * `p_10_dd' + `rr_11_dd' * `p_11_dd' + `rr_12_dd' * `p_12_dd' + `rr_20_dd' * `p_20_dd' + `rr_21_dd' * `p_21_dd' + `rr_22_dd' * `p_22_dd'

local p_00_rr_d=`rr_00_rr' * `p_00_rr' / `sumproduct_rr'
local p_01_rr_d=`rr_01_rr' * `p_01_rr' / `sumproduct_rr'
local p_02_rr_d=`rr_02_rr' * `p_02_rr' / `sumproduct_rr'
local p_10_rr_d=`rr_10_rr' * `p_10_rr' / `sumproduct_rr'
local p_11_rr_d=`rr_11_rr' * `p_11_rr' / `sumproduct_rr'
local p_12_rr_d=`rr_12_rr' * `p_12_rr' / `sumproduct_rr'
local p_20_rr_d=`rr_20_rr' * `p_20_rr' / `sumproduct_rr'
local p_21_rr_d=`rr_21_rr' * `p_21_rr' / `sumproduct_rr'
local p_22_rr_d=`rr_22_rr' * `p_22_rr' / `sumproduct_rr'

local p_00_dr_d=`rr_00_dr' * `p_00_dr' / `sumproduct_dr'
local p_01_dr_d=`rr_01_dr' * `p_01_dr' / `sumproduct_dr'
local p_02_dr_d=`rr_02_dr' * `p_02_dr' / `sumproduct_dr'
local p_10_dr_d=`rr_10_dr' * `p_10_dr' / `sumproduct_dr'
local p_11_dr_d=`rr_11_dr' * `p_11_dr' / `sumproduct_dr'
local p_12_dr_d=`rr_12_dr' * `p_12_dr' / `sumproduct_dr'
local p_20_dr_d=`rr_20_dr' * `p_20_dr' / `sumproduct_dr'
local p_21_dr_d=`rr_21_dr' * `p_21_dr' / `sumproduct_dr'
local p_22_dr_d=`rr_22_dr' * `p_22_dr' / `sumproduct_dr'

local p_00_dd_d=`rr_00_dd' * `p_00_dd' / `sumproduct_dd'
local p_01_dd_d=`rr_01_dd' * `p_01_dd' / `sumproduct_dd'
local p_02_dd_d=`rr_02_dd' * `p_02_dd' / `sumproduct_dd'
local p_10_dd_d=`rr_10_dd' * `p_10_dd' / `sumproduct_dd'
local p_11_dd_d=`rr_11_dd' * `p_11_dd' / `sumproduct_dd'
local p_12_dd_d=`rr_12_dd' * `p_12_dd' / `sumproduct_dd'
local p_20_dd_d=`rr_20_dd' * `p_20_dd' / `sumproduct_dd'
local p_21_dd_d=`rr_21_dd' * `p_21_dd' / `sumproduct_dd'
local p_22_dd_d=`rr_22_dd' * `p_22_dd' / `sumproduct_dd'


set obs 3
gen case_g1=_n-1
gen parent1_g1=_n-1
gen parent2_g1=_n-1
fillin  case_g1 parent1_g1 parent2_g1
drop _fillin
count
gen p_g1_case_parents=.
replace  p_g1_case_parents=0 if ( parent1_g1==2 & case_g1==0)|( parent2_g1==2 &  case_g1==0)
replace  p_g1_case_parents=0 if ( parent1_g1==0 & case_g1==2)|( parent2_g1==0 &  case_g1==2)
replace  p_g1_case_parents=0 if  case_g1>( parent1_g1+ parent2_g1)
replace  p_g1_case_parents=0 if ( parent1_g1==2 & case_g1==0)|( parent2_g1==2 &  case_g1==0)
replace  p_g1_case_parents=0 if ( parent1_g1==0 & case_g1==2)|( parent2_g1==0 &  case_g1==2)
replace  p_g1_case_parents=1 if case_g1==0 &  parent1_g1==0 & parent2_g1==0
replace  p_g1_case_parents=1 if case_g1==2 &  parent1_g1==2 & parent2_g1==2

replace p_g1_case_parents = 0.5 in 2
replace p_g1_case_parents = 0.5 in 4
replace p_g1_case_parents = 0.25 in 5
replace p_g1_case_parents = 0.5 in 11
replace p_g1_case_parents = 1 in 12
replace p_g1_case_parents = 0.5 in 13
replace p_g1_case_parents = 0.5 in 14
replace p_g1_case_parents = 0.5 in 15
replace p_g1_case_parents = 1 in 16
replace p_g1_case_parents = 0.5 in 17
replace p_g1_case_parents = 0 in 18
replace p_g1_case_parents = 0.25 in 23
replace p_g1_case_parents = 0.5 in 24
replace p_g1_case_parents = 0.5 in 26

gen g1_untrans=( parent1_g1+ parent2_g1)- case_g1
replace  g1_untrans=. if  p_g1_case_parents==0
keep if g1_untrans~=.
*drop g1_untrans


gen byte var6 = 0 in 1
gen byte var7 = 0 in 1
gen byte var8 = 0 in 1
replace var6 = 0 in 2
replace var7 = 1 in 2
replace var8 = 1 in 2
replace var6 = 0 in 3
replace var7 = 1 in 3
replace var8 = 1 in 3
replace var6 = 1 in 4
replace var7 = 1 in 4
replace var8 = 2 in 4
replace var6 = 0 in 5
replace var7 = 0 in 5
replace var8 = 1 in 5
replace var6 = 1 in 6
replace var7 = 1 in 6
replace var8 = 1 in 6
replace var6 = 0 in 7
replace var7 = 0 in 7
replace var8 = 1 in 7
replace var6 = 0 in 8
replace var7 = 2 in 8
replace var7 = 1 in 8
replace var8 = 2 in 8
replace var6 = 1 in 9
replace var7 = 2 in 9
replace var8 = 2 in 9
replace var6 = 1 in 10
replace var7 = 1 in 10
replace var8 = 1 in 10
replace var6 = 1 in 11
replace var7 = 2 in 11
replace var8 = 2 in 11
replace var6 = 0 in 12
replace var7 = 1 in 12
replace var8 = 1 in 12
replace var6 = 1 in 13
replace var7 = 1 in 13
replace var8 = 2 in 13
replace var6 = 1 in 14
replace var7 = 1 in 14
replace var8 = 2 in 14
replace var6 = 2 in 15
replace var7 = 2 in 15
replace var8 = 2 in 15

rename var6 g1_control1
rename var7 g1_control2
rename var8 g1_control3

save `gene1', replace

gen case_g2= case_g1
gen parent1_g2=parent1_g1
gen parent2_g2=parent2_g1
gen p_g2_case_parents=p_g1_case_parents
gen g2_control1=g1_control1
gen g2_control2=g1_control2
gen g2_control3=g1_control3
gen g2_untrans=g1_untrans

drop case_g1 parent1_g1 parent2_g1 p_g1_case_parents g1_control1 g1_control2 g1_control3 g1_untrans
keep if g2_untrans~=.

save `gene2', replace

cross using `gene1'

gen p_caseg1g2_parents= p_g1_case_parents* p_g2_case_parents

drop p_g1_case_parents p_g2_case_parents

tempvar p_parents_g1_d p_parents_g2_d p_parents_g1_r p_parents_g2_r

gen `p_parents_g1_d'=0 if g1_untrans==.
replace `p_parents_g1_d'=(1-`af_g1_dom')^2 if g1_untrans==0
replace `p_parents_g1_d'=`af_g1_dom'*(1-`af_g1_dom') if g1_untrans==1
replace `p_parents_g1_d'=`af_g1_dom'^2 if g1_untrans==2
replace `p_parents_g1_d'=2*`af_g1_dom'*(1-`af_g1_dom') if case_g1==1 & parent1_g1==1 & parent2_g1==1

gen `p_parents_g2_d'=0 if g2_untrans==.
replace `p_parents_g2_d'=(1-`af_g2_dom')^2 if g2_untrans==0
replace `p_parents_g2_d'=`af_g2_dom'*(1-`af_g2_dom') if g2_untrans==1
replace `p_parents_g2_d'=`af_g2_dom'^2 if g2_untrans==2
replace `p_parents_g2_d'=2*`af_g2_dom'*(1-`af_g2_dom') if case_g2==1 & parent1_g2==1 & parent2_g2==1

gen `p_parents_g1_r'=0 if g1_untrans==.
replace `p_parents_g1_r'=(1-`af_g1_rec')^2 if g1_untrans==0
replace `p_parents_g1_r'=`af_g1_rec'*(1-`af_g1_rec') if g1_untrans==1
replace `p_parents_g1_r'=`af_g1_rec'^2 if g1_untrans==2
replace `p_parents_g1_r'=2*`af_g1_rec'*(1-`af_g1_rec') if case_g1==1 & parent1_g1==1 & parent2_g1==1

gen `p_parents_g2_r'=0 if g2_untrans==.
replace `p_parents_g2_r'=(1-`af_g2_rec')^2 if g2_untrans==0
replace `p_parents_g2_r'=`af_g2_rec'*(1-`af_g2_rec') if g2_untrans==1
replace `p_parents_g2_r'=`af_g2_rec'^2 if g2_untrans==2
replace `p_parents_g2_r'=2*`af_g2_rec'*(1-`af_g2_rec') if case_g2==1 & parent1_g2==1 & parent2_g2==1

gen p_par_child_dd=`p_parents_g1_d'*`p_parents_g2_d'
gen p_par_child_dr=`p_parents_g1_d'*`p_parents_g2_r'
gen p_par_child_rr=`p_parents_g1_r'*`p_parents_g2_r'

gen p_case_dd=`p_00_dd_d' if case_g1==0 & case_g2==0
replace p_case_dd=`p_01_dd_d' if case_g1==0 & case_g2==1
replace p_case_dd=`p_02_dd_d' if case_g1==0 & case_g2==2
replace p_case_dd=`p_10_dd_d' if case_g1==1 & case_g2==0
replace p_case_dd=`p_11_dd_d' if case_g1==1 & case_g2==1
replace p_case_dd=`p_12_dd_d' if case_g1==1 & case_g2==2
replace p_case_dd=`p_20_dd_d' if case_g1==2 & case_g2==0
replace p_case_dd=`p_21_dd_d' if case_g1==2 & case_g2==1
replace p_case_dd=`p_22_dd_d' if case_g1==2 & case_g2==2

gen p_case_dr=`p_00_dr_d' if case_g1==0 & case_g2==0
replace p_case_dr=`p_01_dr_d' if case_g1==0 & case_g2==1
replace p_case_dr=`p_02_dr_d' if case_g1==0 & case_g2==2
replace p_case_dr=`p_10_dr_d' if case_g1==1 & case_g2==0
replace p_case_dr=`p_11_dr_d' if case_g1==1 & case_g2==1
replace p_case_dr=`p_12_dr_d' if case_g1==1 & case_g2==2
replace p_case_dr=`p_20_dr_d' if case_g1==2 & case_g2==0
replace p_case_dr=`p_21_dr_d' if case_g1==2 & case_g2==1
replace p_case_dr=`p_22_dr_d' if case_g1==2 & case_g2==2

gen p_case_rr=`p_00_rr_d' if case_g1==0 & case_g2==0
replace p_case_rr=`p_01_rr_d' if case_g1==0 & case_g2==1
replace p_case_rr=`p_02_rr_d' if case_g1==0 & case_g2==2
replace p_case_rr=`p_10_rr_d' if case_g1==1 & case_g2==0
replace p_case_rr=`p_11_rr_d' if case_g1==1 & case_g2==1
replace p_case_rr=`p_12_rr_d' if case_g1==1 & case_g2==2
replace p_case_rr=`p_20_rr_d' if case_g1==2 & case_g2==0
replace p_case_rr=`p_21_rr_d' if case_g1==2 & case_g2==1
replace p_case_rr=`p_22_rr_d' if case_g1==2 & case_g2==2

replace p_par_child_dd=p_par_child_dd/2 if case_g1==1
replace p_par_child_dd=p_par_child_dd/2 if case_g2==1

replace p_par_child_dr=p_par_child_dr/2 if case_g1==1
replace p_par_child_dr=p_par_child_dr/2 if case_g2==1

replace p_par_child_rr=p_par_child_rr/2 if case_g1==1
replace p_par_child_rr=p_par_child_rr/2 if case_g2==1

gen p_case_and_par_dd=p_case_dd*p_par_child_dd
gen p_case_and_par_dr=p_case_dr*p_par_child_dr
gen p_case_and_par_rr=p_case_rr*p_par_child_rr

drop if p_case_and_par_dd==0
drop  g2_untrans g1_untrans p_caseg1g2_parents p_par_child_dd p_par_child_rr p_par_child_dr p_case_dd p_case_dr p_case_rr

drop parent1_g2 parent2_g2 parent1_g1 parent2_g1

collapse (sum)  p_case_and_par_dd p_case_and_par_dr p_case_and_par_rr, by( case_g2 g2_control1 g2_control2 g2_control3 case_g1 g1_control1 g1_control2 g1_control3)

gen id=_n
rename case_g1 g1case
rename case_g2 g2case
rename g1_control1 g1control1
rename g1_control2 g1control2
rename g1_control3 g1control3
rename g2_control1 g2control1
rename g2_control2 g2control2
rename g2_control3 g2control3

gen g1control4=g1control1
gen g1control5=g1control2
gen g1control6=g1control3
gen g1control7=g1control1
gen g1control8=g1control2
gen g1control9=g1control3
gen g1control10=g1control1
gen g1control11=g1control2
gen g1control12=g1control3
gen g1control13=g1case
gen g1control14=g1case
gen g1control15=g1case


gen g2control4=g2control2
gen g2control5=g2control2
gen g2control6=g2control2
gen g2control7=g2control3
gen g2control8=g2control3
gen g2control9=g2control3
gen g2control10=g2case
gen g2control11=g2case
gen g2control12=g2case
gen g2control13=g2control1
gen g2control14=g2control2
gen g2control15=g2control3

replace g2control2=g2control1
replace g2control3=g2control1

reshape long g1 g2, i(id) j(caseness) string

gen aff=1 if caseness=="case"
replace aff=0 if aff==.
drop caseness

gen weight_dd=round(p_case_and_par_dd*10000000, 1)
gen weight_dr=round(p_case_and_par_dr*10000000, 1)
gen weight_rr=round(p_case_and_par_rr*10000000, 1)

gen gen_locus1_d=1 if g1>0
replace gen_locus1_d=0 if gen_locus1_d==.

gen gen_locus2_d=1 if g2>0
replace gen_locus2_d=0 if gen_locus2_d==.

gen gen_locus1_r=1 if g1>1
replace gen_locus1_r=0 if gen_locus1_r==.

gen gen_locus2_r=1 if g2>1
replace gen_locus2_r=0 if gen_locus2_r==.

gen ggi_dd=gen_locus1_d*gen_locus2_d
gen ggi_dr=gen_locus1_d*gen_locus2_r
gen ggi_rr=gen_locus1_r*gen_locus2_r

clogit aff gen_locus1_r gen_locus2_r ggi_rr [fweight=weight_rr], group(id) or
lrtest, saving(0)
matrix a =e(b)
local rr_or=exp(a[1,3])
clogit aff gen_locus1_r gen_locus2_r [fweight=weight_rr], group(id) or
lrtest
local alpha_1=invchi2(1, (1-`alpha_1'))
local temp7=`sample_size'*r(chi2)/10000000
local tdt= 1-nchi(1, `temp7', `alpha_1')
local rrtdt=`tdt'*100
local power = `power'/100
local ss_tdt=npnchi2(1, `alpha_1', (1-`power'))
local rrss_tdt=`ss_tdt'*`sample_size'/`temp7'

clogit aff gen_locus1_d gen_locus2_d ggi_dd [fweight=weight_dd], group(id) or
lrtest, saving(0)
matrix a =e(b)
local dd_or=exp(a[1,3])
clogit aff gen_locus1_d gen_locus2_d [fweight=weight_dd], group(id) or
lrtest
local temp7=`sample_size'*r(chi2)/10000000
local tdt= 1-nchi(1, `temp7', `alpha_1')
local ddtdt=`tdt'*100
local ss_tdt=npnchi2(1, `alpha_1', (1-`power'))
local ddss_tdt=`ss_tdt'*`sample_size'/`temp7'

clogit aff gen_locus1_d gen_locus2_r ggi_dr [fweight=weight_dr], group(id) or
lrtest, saving(0)
matrix a =e(b)
local dr_or=exp(a[1,3])
clogit aff gen_locus1_d gen_locus2_r [fweight=weight_dr], group(id) or
lrtest
local temp7=`sample_size'*r(chi2)/10000000
local tdt= 1-nchi(1, `temp7', `alpha_1')
local drtdt=`tdt'*100
local ss_tdt=npnchi2(1, `alpha_1', (1-`power'))
local drss_tdt=`ss_tdt'*`sample_size'/`temp7'

drop _all

use `using'
replace dd_power=`ddtdt' in `loop'
replace dr_power=`drtdt' in `loop'
replace rr_power=`rrtdt' in `loop'
replace dd_ss=`ddss_tdt' in `loop'
replace dr_ss=`drss_tdt' in `loop'
replace rr_ss=`rrss_tdt' in `loop'
replace dd_or=`dd_or' in `loop'
replace dr_or=`dr_or' in `loop'
replace rr_or=`rr_or' in `loop'
save `using', replace

}

replace dd_power=round(dd_power, 0.01)
replace dr_power=round(dr_power, 0.01)
replace rr_power=round(rr_power, 0.01)
replace dd_or=round(dd_or, 0.01)
replace dr_or=round(dr_or, 0.01)
replace rr_or=round(rr_or, 0.01)
replace dd_ss=round(dd_ss, 1)
replace dr_ss=round(dr_ss, 1)
replace rr_ss=round(rr_ss, 1)

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 alpha_1 "significance level"
label var pd "disease frequency"
label var dd_power "case-parent power (%) to detect gene-gene interaction (inheritance=DD)"
label var dr_power "case-parent power (%) to detect gene-gene interaction (inheritance=DR)"
label var rr_power "case-parent power (%) to detect gene-gene interaction (inheritance=RR)"
label var dd_ss "required sample size (inheritance=DD)"
label var dr_ss "required sample size (inheritance=DR)"
label var rr_ss "required sample size (inheritance=RR)"
label var dd_or "calculated interaction odds ratio (inheritance=DD)"
label var dr_or "calculated interaction odds ratio (inheritance=DR)"
label var rr_or "calculated interaction odds ratio (inheritance=RR)"

replace rr_power=. if pop_error==1
replace dd_power=. if pop_error==1
replace dr_power=. if pop_error==1
replace rr_or=. if pop_error==1
replace rr_ss=. if pop_error==1
replace dd_or=. if pop_error==1
replace dd_ss=. if pop_error==1
replace dr_or=. if pop_error==1
replace dr_ss=. if pop_error==1


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(70) pg1
display in gr "susceptibility genotype 2 frequency (pg2):" in ye _col(70) pg2
display in gr "genetic risk factor 1 main effect (rrg1)" in ye _col(70) rrg1
display in gr "genetic risk factor 2 main effect (rrg2):" in ye _col(70) rrg2
display in gr "interaction relative risk (rrint):" in ye _col(70) rrint
display in gr "disease prevalence (pd):" in ye _col(70) pd
di
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
di in gr "with a two-sided significance level (alpha=`alpha_1') 
di
display "for a case-parent design with " in ye `ssize' in gr " cases"
display in gr "power (%) both genetic risk factors dominant inheritance:" in ye _col(70) round(dd_power, 0.01)
display in gr "power (%) both genetic risk factors recessive inheritance:" in ye _col(70) round(rr_power, 0.01)
display in gr "power (%) g1 dominant, g2 recessive:" in ye _col(70) round(dr_power, 0.01)
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
di in gr "dominant inheritance:" in ye _col(50) dd_or
di in gr "recessive inheritance:" in ye _col(50) rr_or
display in gr "g1 dominant, g2 recessive inheritance:" in ye _col(50) dr_or
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 two parental controls per case)
di
di in gr "dominant inheritance " in ye _col(50) dd_ss
di in gr "recessive inheritance " in ye _col(50) rr_ss
display in gr "g1 dominant, g2 recessive inheritance:" in ye _col(50) dr_ss
di
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 "case-parent design 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."
}
drop 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 the case-parent design 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 pop_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 variable pop_error 
di in gr "(population error) will equal 1 rather than 0
}
else {
drop pop_error 

quietly {
save `using', replace
}
}
}



end