*! (SJ3-1: st0032)

capture program drop tdt_geipower
program define tdt_geipower

quietly {

clear

syntax using/

set type double

use `using'
keep pg pe rrg1e0 rrg0e1 rrint ssize power alpha_1 pd

count if pg<=0 | pg>=1 | pe<=0 | pe>=1 | pd<=0 | pd>1 | rrg1e0<=0 | rrg0e1<=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 co_power 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
}
}

gen r_power=.
gen r_ss=.
gen r_or=.
gen d_power=.
gen d_ss=.
gen d_or=.

save `using', replace

gen p_g0e0=(1-pe)*(1-pg)
gen p_g0e1=(1-pg)*pe
gen p_g1e0=pg*(1-pe)
gen p_g1e1=pg*pe
gen rrg0e0=1
gen p_d_g0e0=pd/(p_g0e0* rrg0e0+ p_g0e1* rrg0e1+ p_g1e0* rrg1e0+ rrg1e0* rrg0e1*rrint* p_g1e1)
gen p_d_g0e1= p_d_g0e0* rrg0e1
gen p_d_g1e0= p_d_g0e0* rrg1e0
gen p_d_g1e1= p_d_g0e0* rrg0e1* rrg1e0* rrint

gen pop_error=0
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_g0e0 p_g0e1 p_g1e0 p_g1e1 rrg0e0 p_d_g0e0 p_d_g0e1 p_d_g1e0 p_d_g1e1

save `using', replace

local loop=0

count
local scenar=r(N)
while `loop'~=`scenar' {
local loop=`loop'+1

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

drop _all
use `using'


* population parameters for case parent exemplary dataset
local pg=pg in `loop'
local pe=pe in `loop'
local rrg=rrg1e0 in `loop'
local rre=rrg0e1 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'

local af_rec=`pg'^(0.5)
local af_dom=1-(1-`pg')^(0.5)

display `af_rec'
display `af_dom'


local p_g0e0_r= (1-`af_rec')^2 * (1-`pe')
local p_g1e0_r= 2 * `af_rec' * (1-`af_rec') * (1-`pe')
local p_g2e0_r= `af_rec'^2 * (1-`pe')
local p_g0e1_r= (1-`af_rec')^2 * `pe'
local p_g1e1_r= 2 * `af_rec' * (1-`af_rec') * `pe'
local p_g2e1_r= `af_rec'^2 * `pe'

local p_g0e0_d= (1-`af_dom')^2 * (1-`pe')
local p_g1e0_d= 2 * `af_dom' * (1-`af_dom') * (1-`pe')
local p_g2e0_d= `af_dom'^2 * (1-`pe')
local p_g0e1_d= (1-`af_dom')^2 * `pe'
local p_g1e1_d= 2 * `af_dom' * (1-`af_dom') * `pe'
local p_g2e1_d= `af_dom'^2 * `pe'

local rr_g0e0_r=1
local rr_g1e0_r=1
local rr_g2e0_r=`rrg'
local rr_g0e1_r=`rre'
local rr_g1e1_r=`rre'
local rr_g2e1_r=`rrg'*`rre'*`rrint'

local rr_g0e0_d=1
local rr_g1e0_d=`rrg'
local rr_g2e0_d=`rrg'
local rr_g0e1_d=`rre'
local rr_g1e1_d=`rrg'*`rre'*`rrint'
local rr_g2e1_d=`rrg'*`rre'*`rrint'

local sumproduct_r=`p_g0e0_r'*`rr_g0e0_r'+`p_g1e0_r'*`rr_g1e0_r'+`p_g2e0_r'*`rr_g2e0_r'+`p_g0e1_r'*`rr_g0e1_r'+`p_g1e1_r'*`rr_g1e1_r'+`p_g2e1_r'*`rr_g2e1_r'
local sumproduct_d=`p_g0e0_d'*`rr_g0e0_d'+`p_g1e0_d'*`rr_g1e0_d'+`p_g2e0_d'*`rr_g2e0_d'+`p_g0e1_d'*`rr_g0e1_d'+`p_g1e1_d'*`rr_g1e1_d'+`p_g2e1_d'*`rr_g2e1_d'

local p_g0e0_d_r=`p_g0e0_r'*`rr_g0e0_r'/`sumproduct_r'
local p_g1e0_d_r=`p_g1e0_r'*`rr_g1e0_r'/`sumproduct_r'
local p_g2e0_d_r=`p_g2e0_r'*`rr_g2e0_r'/`sumproduct_r'
local p_g0e1_d_r=`p_g0e1_r'*`rr_g0e1_r'/`sumproduct_r'
local p_g1e1_d_r=`p_g1e1_r'*`rr_g1e1_r'/`sumproduct_r'
local p_g2e1_d_r=`p_g2e1_r'*`rr_g2e1_r'/`sumproduct_r'

local p_g0e0_d_d=`p_g0e0_d'*`rr_g0e0_d'/`sumproduct_d'
local p_g1e0_d_d=`p_g1e0_d'*`rr_g1e0_d'/`sumproduct_d'
local p_g2e0_d_d=`p_g2e0_d'*`rr_g2e0_d'/`sumproduct_d'
local p_g0e1_d_d=`p_g0e1_d'*`rr_g0e1_d'/`sumproduct_d'
local p_g1e1_d_d=`p_g1e1_d'*`rr_g1e1_d'/`sumproduct_d'
local p_g2e1_d_d=`p_g2e1_d'*`rr_g2e1_d'/`sumproduct_d'

drop _all
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

expand 2
gen e_case=1 in 1/15
replace e_case=0 if e_case==.
gen e_control1=1 in 1/15
replace e_control1=0 if e_control1==.
gen e_control2=1 in 1/15
replace e_control2=0 if e_control2==.
gen e_control3=1 in 1/15
replace e_control3=0 if e_control3==.

renpfix g1 g

rename case_g1 g_case

gen p_case_r=`p_g0e0_d_r' if g_case==0 & e_case==0
replace p_case_r=`p_g1e0_d_r' if g_case==1 & e_case==0
replace p_case_r=`p_g2e0_d_r' if g_case==2 & e_case==0
replace p_case_r=`p_g0e1_d_r' if g_case==0 & e_case==1
replace p_case_r=`p_g1e1_d_r' if g_case==1 & e_case==1
replace p_case_r=`p_g2e1_d_r' if g_case==2 & e_case==1

gen p_case_d=`p_g0e0_d_d' if g_case==0 & e_case==0
replace p_case_d=`p_g1e0_d_d' if g_case==1 & e_case==0
replace p_case_d=`p_g2e0_d_d' if g_case==2 & e_case==0
replace p_case_d=`p_g0e1_d_d' if g_case==0 & e_case==1
replace p_case_d=`p_g1e1_d_d' if g_case==1 & e_case==1
replace p_case_d=`p_g2e1_d_d' if g_case==2 & e_case==1

gen p_parents_g_r=0 if g_untrans==.
replace p_parents_g_r=(1-`af_rec')^2 if g_untrans==0
replace p_parents_g_r=`af_rec'*(1-`af_rec') if g_untrans==1
replace p_parents_g_r=`af_rec'^2 if g_untrans==2
replace p_parents_g_r=2*`af_rec'*(1-`af_rec') if g_case==1 & parent1_g1==1 & parent2_g1==1

gen p_parents_g_d=0 if g_untrans==.
replace p_parents_g_d=(1-`af_dom')^2 if g_untrans==0
replace p_parents_g_d=`af_dom'*(1-`af_dom') if g_untrans==1
replace p_parents_g_d=`af_dom'^2 if g_untrans==2
replace p_parents_g_d=2*`af_dom'*(1-`af_dom') if g_case==1 & parent1_g1==1 & parent2_g1==1

replace p_parents_g_d=p_parents_g_d/2 if g_case==1
replace p_parents_g_r=p_parents_g_r/2 if g_case==1

gen p_par_child_r=p_case_r* p_parents_g_r
gen p_par_child_d=p_case_d* p_parents_g_d

drop parent1_g1 parent2_g1 p_g1_case_parents g_untrans p_case_r p_case_d p_parents_g_r p_parents_g_d

gen id=_n
reshape long g e, i(id) j(caseness) string

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

gen g_d=1 if  g>0
replace g_d=0 if g_d~=1
gen g_r=1 if  g>1
replace g_r=0 if g_r~=1
gen gei_r= g_r* e
gen gei_d= g_d* e

gen weight_r=round(1000000* p_par_child_r,1)
gen weight_d=round(1000000* p_par_child_d,1)

local alpha_1=invchi2(1, (1-`alpha_1'))
local power = `power'/100

clogit aff  g_r gei_r [fweight=weight_r], group(id) or
matrix a =e(b)
local r_or=exp(a[1,2])
lrtest, saving(0)
clogit aff  g_r [fweight=weight_r], group(id) or
lrtest
local temp7=`sample_size'*r(chi2)/1000000
local r_tdt= 1-nchi(1, `temp7', `alpha_1')
local r_tdt=`r_tdt'*100
local r_ss_tdt=npnchi2(1, `alpha_1', (1-`power'))
local r_ss_tdt=`r_ss_tdt'*`sample_size'/`temp7'

clogit aff  g_d gei_d [fweight=weight_d], group(id) or
lrtest, saving(0)
matrix a=e(b)
local d_or=exp(a[1,2])
clogit aff  g_d [fweight=weight_d], group(id) or
lrtest
local temp7=`sample_size'*r(chi2)/1000000
local tdt= 1-nchi(1, `temp7', `alpha_1')
local tdt=`tdt'*100
local ss_tdt=npnchi2(1, `alpha_1', (1-`power'))
local ss_tdt=`ss_tdt'*`sample_size'/`temp7'

drop _all
use `using'
replace r_power=`r_tdt' in `loop'
replace r_ss=`r_ss_tdt' in `loop'
replace d_power=`tdt' in `loop'
replace d_ss=`ss_tdt' in `loop'
replace r_or=`r_or' in `loop'
replace d_or=`d_or' in `loop'

save `using', replace


}

label var pg "P(G)"
label var pe "P(E)"
label var rrg0e1 "RR(G0E1)"
label var rrg1e0 "RR(G1E0)"
label var rrint "RR(INT)"
label var alpha_1 "significance level"
label var pd "disease frequency"
label var d_power "Power (%) parent control (dominant inheritance)"
label var r_power "Power (%) parent control (recessive inheritance)"
label var r_or "Interaction odds ratio, parent control (recessive inheritance)"
label var d_or "Interaction odds ratio, parent control (dominant inheritance)"
label var d_ss "sample size, parent control (dominant inheritance)"
label var r_ss "sample size, parent control (recessive inheritance)"

replace d_power=round(d_power, 0.01)
replace r_power=round(r_power, 0.01)
replace d_or=round(d_or, 0.01)
replace r_or=round(r_or, 0.01)
replace d_ss=round(d_ss, 1)
replace r_ss=round(r_ss,1)

replace d_power=. if pop_error==1
replace r_power=. if pop_error==1
replace d_or=. if pop_error==1
replace r_or=. if pop_error==1
replace d_ss=. if pop_error==1
replace r_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
di in gr "suscpetibility genotype frequency (pg):" in ye _col(70) pg
di in gr "environmental risk factor frequency (pe):" in ye _col(70) pe
di in gr "genetic risk factor main effect (rrg1e0):" in ye _col(70) rrg1e0
di in gr "environmental risk factor main effect (rrg0e1):" in ye _col(70) rrg0e1
di in gr "interaction relative risk (rrint):" in ye _col(70) rrint
di 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 (%) dominant inheritance:" in ye _col(35) round(d_power, 0.01)
display in gr "power (%) recessive inheritance:" in ye _col(35) round(r_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(35) d_or
di in gr "recessive inheritance " in ye _col(35) r_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(35) d_ss
di in gr "recessive inheritance " in ye _col(35) r_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