*! (SJ3-1: st0032)

capture program drop gei_matching
program define gei_matching
quietly {
version 7
set more off
syntax using/
drop _all
capture set mem 16m
set type double
use `using'
capture gen rrg0e0=1
capture gen orge=1
keep pg pe rrg0e1 rrg1e0 rrint pd ssize power alpha_1 sensg sense specg spece orge rrg0e0 dm
count if pg<=0 | pg>=1 | pe<=0 | pe>=1 | pd<=0 | pd>1 | rrg1e0<=0 | rrg0e1<=0| rrint<=0 | orge<=0| alpha_1 <=0 |alpha_1 >=1 |power<5 |power>99 | spece<0 |spece>1|sense<0|sense>1|specg<0 |specg>1|sensg<0|sensg>1
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 gei_matching 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 cmge_power=.
gen cmge_ss=.
gen cmge_or=.
gen cme_power=.
gen cme_ss=.
gen cme_or=.
gen flme_power=.
gen flme_ss=.
gen flme_or=.
gen flmse_power=.
gen flmse_ss=.
gen flmse_or=.
gen flmsge_power=.
gen flmsge_ss=.
gen flmsge_or=.

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 orge "OR(GE)"
label var alpha_1 "significance level"
label var pd "disease frequency"
label var sense "sensitivity of exposure surrogate"
label var sensg "sensitivity of genotype surrogate"
label var spece "specificity of exposure surrogate"
label var specg "specificity of genotype surrogate"
label var cmge_power "power (%) counter-matching on surrogate G and E"
label var cmge_ss "required sample size, counter-matching on surrogate G and E"
label var cmge_or "interaction odds ratio, counter-matching on surrogate G and E"
label var cme_power "power (%) counter-matching on surrogate E"
label var cme_ss "required sample size, counter-matching on surrogate E"
label var cme_or "interaction odds ratio, counter-matching on surrogate E"
label var flme_power "power (%) flexible matching on E"
label var flme_ss "required sample size, flexible-matching on E"
label var flme_or "interaction odds ratio, flexible-matching on E"
label var flmse_power "power (%) flexible matching on surrogate E"
label var flmse_ss "required sample size, flexible-matching on surrogate E"
label var flmse_or "interaction odds ratio, flexible-matching on surrogate E"
label var flmsge_power "power (%) flexible matching on surrogate G and E"
label var flmsge_ss "required sample size, flexible-matching on surrogate G and E"
label var flmsge_or "interaction odds ratio, flexible-matching on surrogate G and E"

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 a=orge-1
gen b=-2*orge+orge*pg+orge*pe-pg-pe+1
gen c=orge*(1-pe-pg+pe*pg)

replace p_g0e0=(-1*b-(b^2-(4*a*c))^0.5)/(2*a) if orge~=1
replace p_g0e1=1-pg-p_g0e0 if orge~=1
replace p_g1e0=1-pe-p_g0e0 if orge~=1
replace p_g1e1=1-p_g0e0-p_g0e1-p_g1e0 if orge~=1

drop a b c

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_g0e0<0 | p_d_g0e1>1 | p_d_g0e1<0 | p_d_g1e0>1 | p_d_g1e0<0 | p_d_g1e1>1 | p_d_g1e1<0 

gen p_g0e0_d= p_g0e0*  p_d_g0e0/ pd
gen p_g0e1_d= p_g0e1*  p_d_g0e1/ pd
gen p_g1e0_d= p_g1e0*  p_d_g1e0/ pd
gen p_g1e1_d= p_g1e1*  p_d_g1e1/ pd
gen p_g0e0_not_d= p_g0e0*(1- p_d_g0e0)/(1- pd)
gen p_g0e1_not_d= p_g0e1*(1- p_d_g0e1)/(1- pd)
gen p_g1e0_not_d= p_g1e0*(1- p_d_g1e0)/(1- pd)
gen p_g1e1_not_d= p_g1e1*(1- p_d_g1e1)/(1- pd)

gen pe_cases=p_g0e1_d+p_g1e1_d
gen pe_controls=p_g0e1_not_d+p_g1e1_not_d
gen pe_matched=dm*(pe_cases-pe)+pe
gen p_g0e0_not_d_fm= p_g0e0_not_d*(1-pe_matched)/(1-pe_controls)
gen p_g0e1_not_d_fm= p_g0e1_not_d*pe_matched/pe_controls
gen p_g1e0_not_d_fm= p_g1e0_not_d*(1-pe_matched)/(1-pe_controls)
gen p_g1e1_not_d_fm= p_g1e1_not_d*pe_matched/pe_controls
drop pe_cases pe_controls

gen spg=sensg*pg+(1-specg)*(1-pg)
gen spe=sense*pe+(1-spece)*(1-pe)

gen casessurrg0e0actualg0e0=p_g0e0_d*specg*spece
gen casessurrg0e1actualg0e0=p_g0e0_d*specg*(1-spece)
gen casessurrg1e0actualg0e0=p_g0e0_d*(1-specg)*spece
gen casessurrg1e1actualg0e0=p_g0e0_d*(1-specg)*(1-spece)
gen casessurrg0e0actualg0e1=p_g0e1_d*specg*(1-sense)
gen casessurrg0e1actualg0e1=p_g0e1_d*specg*sense
gen casessurrg1e0actualg0e1=p_g0e1_d*(1-specg)*(1-sense)
gen casessurrg1e1actualg0e1=p_g0e1_d*(1-specg)*sense
gen casessurrg0e0actualg1e0=p_g1e0_d*(1-sensg)*spece
gen casessurrg0e1actualg1e0=p_g1e0_d*(1-sensg)*(1-spece)
gen casessurrg1e0actualg1e0=p_g1e0_d*sensg*spece
gen casessurrg1e1actualg1e0=p_g1e0_d*sensg*(1-spece)
gen casessurrg0e0actualg1e1=p_g1e1_d*(1-sensg)*(1-sense)
gen casessurrg0e1actualg1e1=p_g1e1_d*(1-sensg)*sense
gen casessurrg1e0actualg1e1=p_g1e1_d*sensg*(1-sense)
gen casessurrg1e1actualg1e1=p_g1e1_d*sensg*sense

gen controlsurrg0e0actualg0e0=(p_g0e0_not_d*specg*spece)
gen controlsurrg0e1actualg0e0=(p_g0e0_not_d*specg*(1-spece))
gen controlsurrg1e0actualg0e0=(p_g0e0_not_d*(1-specg)*spece)
gen controlsurrg1e1actualg0e0=(p_g0e0_not_d*(1-specg)*(1-spece))
gen controlsurrg0e0actualg0e1=(p_g0e1_not_d*specg*(1-sense))
gen controlsurrg0e1actualg0e1=(p_g0e1_not_d*specg*sense)
gen controlsurrg1e0actualg0e1=(p_g0e1_not_d*(1-specg)*(1-sense))
gen controlsurrg1e1actualg0e1=(p_g0e1_not_d*(1-specg)*sense)
gen controlsurrg0e0actualg1e0=(p_g1e0_not_d*(1-sensg)*spece)
gen controlsurrg0e1actualg1e0=(p_g1e0_not_d*(1-sensg)*(1-spece))
gen controlsurrg1e0actualg1e0=(p_g1e0_not_d*sensg*spece)
gen controlsurrg1e1actualg1e0=(p_g1e0_not_d*sensg*(1-spece))
gen controlsurrg0e0actualg1e1=(p_g1e1_not_d*(1-sensg)*(1-sense))
gen controlsurrg0e1actualg1e1=(p_g1e1_not_d*(1-sensg)*sense)
gen controlsurrg1e0actualg1e1=(p_g1e1_not_d*sensg*(1-sense))
gen controlsurrg1e1actualg1e1=(p_g1e1_not_d*sensg*sense)

gen controlsurrg0e0=controlsurrg0e0actualg0e0+controlsurrg0e0actualg0e1+controlsurrg0e0actualg1e0+controlsurrg0e0actualg1e1
gen controlsurrg0e1=controlsurrg0e1actualg0e0+controlsurrg0e1actualg0e1+controlsurrg0e1actualg1e0+controlsurrg0e1actualg1e1
gen controlsurrg1e0=controlsurrg1e0actualg0e0+controlsurrg1e0actualg0e1+controlsurrg1e0actualg1e0+controlsurrg1e0actualg1e1
gen controlsurrg1e1=controlsurrg1e1actualg0e0+controlsurrg1e1actualg0e1+controlsurrg1e1actualg1e0+controlsurrg1e1actualg1e1

gen casessurrg0e0=casessurrg0e0actualg0e0+casessurrg0e0actualg0e1+casessurrg0e0actualg1e0+casessurrg0e0actualg1e1
gen casessurrg0e1=casessurrg0e1actualg0e0+casessurrg0e1actualg0e1+casessurrg0e1actualg1e0+casessurrg0e1actualg1e1
gen casessurrg1e0=casessurrg1e0actualg0e0+casessurrg1e0actualg0e1+casessurrg1e0actualg1e0+casessurrg1e0actualg1e1
gen casessurrg1e1=casessurrg1e1actualg0e0+casessurrg1e1actualg0e1+casessurrg1e1actualg1e0+casessurrg1e1actualg1e1

gen surrg0e0=casessurrg0e0*pd+controlsurrg0e0*(1-pd)
gen surrg0e1=casessurrg0e1*pd+controlsurrg0e1*(1-pd)
gen surrg1e0=casessurrg1e0*pd+controlsurrg1e0*(1-pd)
gen surrg1e1=casessurrg1e1*pd+controlsurrg1e1*(1-pd)

replace controlsurrg0e0actualg0e0=(p_g0e0_not_d*specg*spece)/controlsurrg0e0
replace controlsurrg0e1actualg0e0=(p_g0e0_not_d*specg*(1-spece))/controlsurrg0e1
replace controlsurrg1e0actualg0e0=(p_g0e0_not_d*(1-specg)*spece)/controlsurrg1e0
replace controlsurrg1e1actualg0e0=(p_g0e0_not_d*(1-specg)*(1-spece))/controlsurrg1e1
replace controlsurrg0e0actualg0e1=(p_g0e1_not_d*specg*(1-sense))/controlsurrg0e0
replace controlsurrg0e1actualg0e1=(p_g0e1_not_d*specg*sense)/controlsurrg0e1
replace controlsurrg1e0actualg0e1=(p_g0e1_not_d*(1-specg)*(1-sense))/controlsurrg1e0
replace controlsurrg1e1actualg0e1=(p_g0e1_not_d*(1-specg)*sense)/controlsurrg1e1
replace controlsurrg0e0actualg1e0=(p_g1e0_not_d*(1-sensg)*spece)/controlsurrg0e0
replace controlsurrg0e1actualg1e0=(p_g1e0_not_d*(1-sensg)*(1-spece))/controlsurrg0e1
replace controlsurrg1e0actualg1e0=(p_g1e0_not_d*sensg*spece)/controlsurrg1e0
replace controlsurrg1e1actualg1e0=(p_g1e0_not_d*sensg*(1-spece))/controlsurrg1e1
replace controlsurrg0e0actualg1e1=(p_g1e1_not_d*(1-sensg)*(1-sense))/controlsurrg0e0
replace controlsurrg0e1actualg1e1=(p_g1e1_not_d*(1-sensg)*sense)/controlsurrg0e1
replace controlsurrg1e0actualg1e1=(p_g1e1_not_d*sensg*(1-sense))/controlsurrg1e0
replace controlsurrg1e1actualg1e1=(p_g1e1_not_d*sensg*sense)/controlsurrg1e1

gen casesurre0actualg0e0=spece* p_g0e0_d
gen casesurre1actualg0e0=(1-spece)* p_g0e0_d
gen casesurre0actualg0e1=(1-sense)* p_g0e1_d
gen casesurre1actualg0e1=sense* p_g0e1_d
gen casesurre0actualg1e0=spece* p_g1e0_d
gen casesurre1actualg1e0=(1-spece)* p_g1e0_d
gen casesurre0actualg1e1=(1-sense)* p_g1e1_d
gen casesurre1actualg1e1=sense* p_g1e1_d

gen controlsurre0actualg0e0=spece* p_g0e0_not_d
gen controlsurre1actualg0e0=(1-spece)* p_g0e0_not_d
gen controlsurre0actualg0e1=(1-sense)* p_g0e1_not_d
gen controlsurre1actualg0e1=sense* p_g0e1_not_d
gen controlsurre0actualg1e0=spece* p_g1e0_not_d
gen controlsurre1actualg1e0=(1-spece)* p_g1e0_not_d
gen controlsurre0actualg1e1=(1-sense)* p_g1e1_not_d
gen controlsurre1actualg1e1=sense* p_g1e1_not_d

egen controlsurre0=rsum(controlsurre0actualg0e0 controlsurre0actualg0e1 controlsurre0actualg1e0 controlsurre0actualg1e1)
egen controlsurre1=rsum(controlsurre1actualg0e0 controlsurre1actualg0e1 controlsurre1actualg1e0 controlsurre1actualg1e1)

replace controlsurre0actualg0e0=spece* p_g0e0_not_d/controlsurre0
replace controlsurre1actualg0e0=(1-spece)* p_g0e0_not_d/controlsurre1
replace controlsurre0actualg0e1=(1-sense)* p_g0e1_not_d/controlsurre0
replace controlsurre1actualg0e1=sense* p_g0e1_not_d/controlsurre1
replace controlsurre0actualg1e0=spece* p_g1e0_not_d/controlsurre0
replace controlsurre1actualg1e0=(1-spece)* p_g1e0_not_d/controlsurre1
replace controlsurre0actualg1e1=(1-sense)* p_g1e1_not_d/controlsurre0
replace controlsurre1actualg1e1=sense* p_g1e1_not_d/controlsurre1

egen casesurre0=rsum(casesurre0actualg0e0 casesurre0actualg0e1 casesurre0actualg1e0 casesurre0actualg1e1)
egen casesurre1=rsum(casesurre1actualg0e0 casesurre1actualg0e1 casesurre1actualg1e0 casesurre1actualg1e1)

gen surre0=casesurre0*pd+controlsurre0*(1-pd)
gen surre1=casesurre1*pd+controlsurre1*(1-pd)

gen pe_matchedsurr=dm*(casesurre1-spe)+spe
gen p_g0e0_not_d_fmse= controlsurre0actualg0e0*(1-pe_matchedsurr) + controlsurre1actualg0e0*pe_matchedsurr
gen p_g0e1_not_d_fmse= controlsurre0actualg0e1*(1-pe_matchedsurr) + controlsurre1actualg0e1*pe_matchedsurr
gen p_g1e0_not_d_fmse= controlsurre0actualg1e0*(1-pe_matchedsurr) + controlsurre1actualg1e0*pe_matchedsurr
gen p_g1e1_not_d_fmse= controlsurre0actualg1e1*(1-pe_matchedsurr) + controlsurre1actualg1e1*pe_matchedsurr

gen pg_matchedsurr=dm*((casessurrg1e0+casessurrg1e1)-spg)+spg
gen p_g0e0_not_d_fmsge=controlsurrg0e0actualg0e0*(1-pe_matchedsurr)*(1-pg_matchedsurr)+controlsurrg0e1actualg0e0*pe_matchedsurr*(1-pg_matchedsurr)+controlsurrg1e0actualg0e0*(1-pe_matchedsurr)*pg_matchedsurr+controlsurrg1e1actualg0e0*pe_matchedsurr*pg_matchedsurr
gen p_g0e1_not_d_fmsge=controlsurrg0e0actualg0e1*(1-pe_matchedsurr)*(1-pg_matchedsurr)+controlsurrg0e1actualg0e1*pe_matchedsurr*(1-pg_matchedsurr)+controlsurrg1e0actualg0e1*(1-pe_matchedsurr)*pg_matchedsurr+controlsurrg1e1actualg0e1*pe_matchedsurr*pg_matchedsurr
gen p_g1e0_not_d_fmsge=controlsurrg0e0actualg1e0*(1-pe_matchedsurr)*(1-pg_matchedsurr)+controlsurrg0e1actualg1e0*pe_matchedsurr*(1-pg_matchedsurr)+controlsurrg1e0actualg1e0*(1-pe_matchedsurr)*pg_matchedsurr+controlsurrg1e1actualg1e0*pe_matchedsurr*pg_matchedsurr
gen p_g1e1_not_d_fmsge=controlsurrg0e0actualg1e1*(1-pe_matchedsurr)*(1-pg_matchedsurr)+controlsurrg0e1actualg1e1*pe_matchedsurr*(1-pg_matchedsurr)+controlsurrg1e0actualg1e1*(1-pe_matchedsurr)*pg_matchedsurr+controlsurrg1e1actualg1e1*pe_matchedsurr*pg_matchedsurr

gen match_error=0
replace match_error=1 if pe_matched>1 | pg_matchedsurr>1 | pe_matchedsurr>1|pe_matched<0 | pg_matchedsurr<0 | pe_matchedsurr<0
drop pe_matchedsurr pg_matchedsurr pe_matched


save `using', replace


local loop2=0

count
local no_in_data=r(N)

while `loop2'~=`no_in_data' {
local loop2=`loop2'+1

noisily {
display "calculating observation `loop2' of `no_in_data'"
}

use `using'

local alpha_1=alpha_1 in `loop2'
local power=power in `loop2'
local ssize=ssize in `loop2'

local casessurrg0e0actualg0e0= casessurrg0e0actualg0e0 in `loop2'
local casessurrg0e1actualg0e0= casessurrg0e1actualg0e0 in `loop2'
local casessurrg1e0actualg0e0= casessurrg1e0actualg0e0 in `loop2'
local casessurrg1e1actualg0e0= casessurrg1e1actualg0e0 in `loop2'
local casessurrg0e0actualg0e1= casessurrg0e0actualg0e1 in `loop2'
local casessurrg0e1actualg0e1= casessurrg0e1actualg0e1 in `loop2'
local casessurrg1e0actualg0e1= casessurrg1e0actualg0e1 in `loop2'
local casessurrg1e1actualg0e1= casessurrg1e1actualg0e1 in `loop2'
local casessurrg0e0actualg1e0= casessurrg0e0actualg1e0 in `loop2'
local casessurrg0e1actualg1e0= casessurrg0e1actualg1e0 in `loop2'
local casessurrg1e0actualg1e0= casessurrg1e0actualg1e0 in `loop2'
local casessurrg1e1actualg1e0= casessurrg1e1actualg1e0 in `loop2'
local casessurrg0e0actualg1e1= casessurrg0e0actualg1e1 in `loop2'
local casessurrg0e1actualg1e1= casessurrg0e1actualg1e1 in `loop2'
local casessurrg1e0actualg1e1= casessurrg1e0actualg1e1 in `loop2'
local casessurrg1e1actualg1e1= casessurrg1e1actualg1e1 in `loop2'

local surrg0e0=surrg0e0 in `loop2'
local surrg0e1=surrg0e1 in `loop2'
local surrg1e0=surrg1e0 in `loop2'
local surrg1e1=surrg1e1 in `loop2'

local controlsurrg0e0=controlsurrg0e0 in `loop2'
local controlsurrg0e1=controlsurrg0e1 in `loop2'
local controlsurrg1e0=controlsurrg1e0 in `loop2'
local controlsurrg1e1=controlsurrg1e1 in `loop2'

local controlssurrg0e0actualg0e0= controlsurrg0e0actualg0e0 in `loop2'
local controlssurrg0e1actualg0e0= controlsurrg0e1actualg0e0 in `loop2'
local controlssurrg1e0actualg0e0= controlsurrg1e0actualg0e0 in `loop2'
local controlssurrg1e1actualg0e0= controlsurrg1e1actualg0e0 in `loop2'
local controlssurrg0e0actualg0e1= controlsurrg0e0actualg0e1 in `loop2'
local controlssurrg0e1actualg0e1= controlsurrg0e1actualg0e1 in `loop2'
local controlssurrg1e0actualg0e1= controlsurrg1e0actualg0e1 in `loop2'
local controlssurrg1e1actualg0e1= controlsurrg1e1actualg0e1 in `loop2'
local controlssurrg0e0actualg1e0= controlsurrg0e0actualg1e0 in `loop2'
local controlssurrg0e1actualg1e0= controlsurrg0e1actualg1e0 in `loop2'
local controlssurrg1e0actualg1e0= controlsurrg1e0actualg1e0 in `loop2'
local controlssurrg1e1actualg1e0= controlsurrg1e1actualg1e0 in `loop2'
local controlssurrg0e0actualg1e1= controlsurrg0e0actualg1e1 in `loop2'
local controlssurrg0e1actualg1e1= controlsurrg0e1actualg1e1 in `loop2'
local controlssurrg1e0actualg1e1= controlsurrg1e0actualg1e1 in `loop2'
local controlssurrg1e1actualg1e1= controlsurrg1e1actualg1e1 in `loop2'

local casessurre0actualg0e0=casesurre0actualg0e0 in `loop2'
local casessurre1actualg0e0=casesurre1actualg0e0 in `loop2'
local casessurre0actualg0e1=casesurre0actualg0e1 in `loop2'
local casessurre1actualg0e1=casesurre1actualg0e1 in `loop2'
local casessurre0actualg1e0=casesurre0actualg1e0 in `loop2'
local casessurre1actualg1e0=casesurre1actualg1e0 in `loop2'
local casessurre0actualg1e1=casesurre0actualg1e1 in `loop2'
local casessurre1actualg1e1=casesurre1actualg1e1 in `loop2'

local controlssurre0actualg0e0=controlsurre0actualg0e0 in `loop2'
local controlssurre1actualg0e0=controlsurre1actualg0e0 in `loop2'
local controlssurre0actualg0e1=controlsurre0actualg0e1 in `loop2'
local controlssurre1actualg0e1=controlsurre1actualg0e1 in `loop2'
local controlssurre0actualg1e0=controlsurre0actualg1e0 in `loop2'
local controlssurre1actualg1e0=controlsurre1actualg1e0 in `loop2'
local controlssurre0actualg1e1=controlsurre0actualg1e1 in `loop2'
local controlssurre1actualg1e1=controlsurre1actualg1e1 in `loop2'

local surre1=surre1 in `loop2'
local surre0=surre0 in `loop2'

local p_g0e0_d=p_g0e0_d in `loop2'
local p_g0e1_d=p_g0e1_d in `loop2'
local p_g1e0_d=p_g1e0_d in `loop2'
local p_g1e1_d=p_g1e1_d in `loop2'

local p_g0e0_not_d_fm=p_g0e0_not_d_fm in `loop2'
local p_g0e1_not_d_fm=p_g0e1_not_d_fm in `loop2'
local p_g1e0_not_d_fm=p_g1e0_not_d_fm in `loop2'
local p_g1e1_not_d_fm=p_g1e1_not_d_fm in `loop2'

local p_g0e0_not_d_fmse=p_g0e0_not_d_fmse in `loop2'
local p_g0e1_not_d_fmse=p_g0e1_not_d_fmse in `loop2'
local p_g1e0_not_d_fmse=p_g1e0_not_d_fmse in `loop2'
local p_g1e1_not_d_fmse=p_g1e1_not_d_fmse in `loop2'

local p_g0e0_not_d_fmsge=p_g0e0_not_d_fmsge in `loop2'
local p_g0e1_not_d_fmsge=p_g0e1_not_d_fmsge in `loop2'
local p_g1e0_not_d_fmsge=p_g1e0_not_d_fmsge in `loop2'
local p_g1e1_not_d_fmsge=p_g1e1_not_d_fmsge in `loop2'

drop _all
set obs 1024
gsurrcase                
gactualcase                   
esurrcase                   
eactualcase                    
gsurrcontrol1                   
gactualcontrol1                 
esurrcontrol1                
eactualcontrol1                  
gsurrcontrol2                  
gactualcontrol2                
esurrcontrol2                   
eactualcontrol2                 
gsurrcontrol3                     
gactualcontrol3                   
esurrcontrol3                    
eactualcontrol3   

gen str1 strgcasesurr=string(gsurrcase)
gen str1 strgcaseactual=string(gactualcase)
gen str1 strecasesurr=string(esurrcase)
gen str1 strecaseactual=string(eactualcase)
gen str1 strgcontrol1surr=string(gsurrcontrol1)
gen str1 strgcontrol1actual=string(gactualcontrol1)
gen str1 strecontrol1surr=string(esurrcontrol1)
gen str1 strecontrol1actual=string(eactualcontrol1)
gen str1 strgcontrol2surr=string(gsurrcontrol2)
gen str1 strgcontrol2actual=string(gactualcontrol2)
gen str1 strecontrol2surr=string(esurrcontrol2)
gen str1 strecontrol2actual=string(eactualcontrol2)
gen str1 strgcontrol3surr=string(gsurrcontrol3)
gen str1 strgcontrol3actual=string(gactualcontrol3)
gen str1 strecontrol3surr=string(esurrcontrol3)
gen str1 strecontrol3actual=string(eactualcontrol3)

gen str32 case="casessurrg"+ strgcasesurr+"e"+ strecasesurr+"actualg"+strgcaseactual+"e"+ strecaseactual

gen str32 control1="controlssurrg"+ strgcontrol1surr+"e"+ strecontrol1surr+"actualg"+strgcontrol1actual+"e"+ strecontrol1actual
gen str32 control2="controlssurrg"+ strgcontrol2surr+"e"+ strecontrol2surr+"actualg"+strgcontrol2actual+"e"+ strecontrol2actual
gen str32 control3="controlssurrg"+ strgcontrol3surr+"e"+ strecontrol3surr+"actualg"+strgcontrol3actual+"e"+ strecontrol3actual

drop  strgcasesurr strgcaseactual strecasesurr strecaseactual strgcontrol1surr strgcontrol1actual strecontrol1surr  strecontrol1actual strgcontrol2surr strgcontrol2actual strecontrol2surr strecontrol2actual strgcontrol3surr strgcontrol3actual strecontrol3surr strecontrol3actual

gen p_actcase_actcontrols=. 

count
local datasetsize=r(N)

local loop1=0
while `loop1'~=`datasetsize'{
local loop1=`loop1'+1
local temp1=case in `loop1'
local temp2=control1 in `loop1'
local temp3=control2 in `loop1'
local temp4=control3 in `loop1'
local temp5=``temp1''*``temp2''*``temp3''*``temp4''
replace p_actcase_actcontrols=`temp5' in `loop1'
}

drop  case control1 control2 control3
gen id=_n
reshape long gsurr gactual esurr eactual, i(id) j(caseness) string
gen offset=.
replace offset=`surrg0e0' if gsurr==0 & esurr==0
replace offset=`surrg0e1' if gsurr==0 & esurr==1
replace offset=`surrg1e0' if gsurr==1 & esurr==0
replace offset=`surrg1e1' if gsurr==1 & esurr==1
gen gei=gactual*eactual
gen no_sets= round(p_actcase_actcontrols*10^9, 1)
gen aff=1 if caseness=="case"
replace aff=0 if aff~=1
gen offset2=ln(offset)
clogit  aff  gactual eactual gei [fweight= no_sets], group(id) or offset(offset2)
lrtest, saving(0)
matrix a=e(b)
local cmge_or=exp(a[1,3])
clogit  aff  gactual eactual [fweight= no_sets], group(id) or offset(offset2)
lrtest
local lrtest=r(chi2) 
local chi_sl=invchi2(1, (1-`alpha_1'))
local ncp=`lrtest'*`ssize'/10^9
local cmge_power=(1-nchi2(1, `ncp', `chi_sl'))*100
local rec_ncp=npnchi2(1, `chi_sl', (1-`power'/100))
local cmge_ss=10^9*(`rec_ncp'/`lrtest')

drop _all
set obs 2

gen gactualcase1=_n-1
gen esurrcase1=_n-1
gen eactualcase1=_n-1
gen gactualcontrol1=_n-1
gen esurrcontrol1=_n-1
gen eactualcontrol1=_n-1
gen gactualcase2=_n-1
gen esurrcase2=_n-1
gen eactualcase2=_n-1
gen gactualcontrol2=_n-1
gen esurrcontrol2=_n-1
gen eactualcontrol2=_n-1

fillin  gactualcase1 esurrcase1 eactualcase1 gactualcontrol1 esurrcontrol1 eactualcontrol1 gactualcase2 esurrcase2 eactualcase2 gactualcontrol2 esurrcontrol2 eactualcontrol2
drop _fillin

keep if esurrcase1~=esurrcontrol1 & esurrcase2~=esurrcontrol2

gen str1 strgcase1actual=string(gactualcase1)
gen str1 strecase1surr=string(esurrcase1)
gen str1 strecase1actual=string(eactualcase1)
gen str1 strgcontrol1actual=string(gactualcontrol1)
gen str1 strecontrol1surr=string(esurrcontrol1)
gen str1 strecontrol1actual=string(eactualcontrol1)
gen str1 strgcase2actual=string(gactualcase2)
gen str1 strecase2surr=string(esurrcase2)
gen str1 strecase2actual=string(eactualcase2)
gen str1 strgcontrol2actual=string(gactualcontrol2)
gen str1 strecontrol2surr=string(esurrcontrol2)
gen str1 strecontrol2actual=string(eactualcontrol2)

gen str32 case1="casessurre"+strecase1surr+"actualg"+strgcase1actual+"e"+ strecase1actual
gen str32 control1="controlssurre"+ strecontrol1surr+"actualg"+strgcontrol1actual+"e"+ strecontrol1actual
gen str32 case2="casessurre"+strecase2surr+"actualg"+strgcase2actual+"e"+ strecase2actual
gen str32 control2="controlssurre"+ strecontrol2surr+"actualg"+strgcontrol2actual+"e"+ strecontrol2actual

gen p_actcase_actcontrols=. 

count
local datasetsize=r(N)

local loop1=0
while `loop1'~=`datasetsize'{
local loop1=`loop1'+1
local temp1=case1 in `loop1'
local temp2=control1 in `loop1'
local temp3=case2 in `loop1'
local temp4=control2 in `loop1'
local temp5=``temp1''*``temp2''*``temp3''*``temp4''
replace p_actcase_actcontrols=`temp5' in `loop1'
}

drop  control2 case2 control1 case1 strecontrol2actual strecontrol2surr strgcontrol2actual strecase2actual strecase2surr strgcase2actual strecontrol1actual strecontrol1surr strgcontrol1actual strecase1actual strecase1surr strgcase1actual

gen id=_n
reshape long gactual esurr eactual, i(id) j(caseness) string
gen offset=.
replace offset=`surre0' if esurr==0
replace offset=`surre1' if esurr==1
gen gei=gactual*eactual
gen no_sets= round(p_actcase_actcontrols*1000000000, 1)
gen aff=1 if caseness=="case1"
replace aff=1 if caseness=="case2"
replace aff=0 if aff~=1
gen offset2=ln(offset)
clogit  aff  gactual eactual gei [fweight= no_sets], group(id) or offset(offset2)
lrtest, saving(0)
matrix a=e(b)
local cme_or=exp(a[1,3])
clogit  aff  gactual eactual [fweight= no_sets], group(id) or offset(offset2)
lrtest
local lrtest=r(chi2)
local ncp=`lrtest'*`ssize'/10^9
local cme_power=(1-nchi2(1, `ncp', `chi_sl'))*100
local rec_ncp=npnchi2(1, `chi_sl', (1-`power'/100))
local cme_ss=10^9*(`rec_ncp'/`lrtest')

local aff1=`p_g0e0_d'
local tot1=`p_g0e0_not_d_fm'
local tot1=`tot1'+`aff1'
local aff2=`p_g0e1_d'
local tot2=`p_g0e1_not_d_fm'
local tot2=`tot2'+`aff2'
local aff3=`p_g1e0_d'
local tot3=`p_g1e0_not_d_fm'
local tot3=`tot3'+`aff3'
local aff4=`p_g1e1_d'
local tot4=`p_g1e1_not_d_fm'
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*1000000000, 1)
replace tot=round(tot*1000000000, 1)
blogit aff tot g e gei
lrtest, saving(0)
matrix A=e(b)
local flme_or=exp(el(A, 1, 3))
blogit aff tot g e
lrtest
local lrtest=r(chi2) 
local ncp=`lrtest'*`ssize'/10^9
local flme_power=(1-nchi2(1, `ncp', `chi_sl'))*100
local rec_ncp=npnchi2(1, `chi_sl', (1-`power'/100))
local flme_ss=10^9*(`rec_ncp'/`lrtest')

local aff1=`p_g0e0_d'
local tot1=`p_g0e0_not_d_fmse'
local tot1=`tot1'+`aff1'
local aff2=`p_g0e1_d'
local tot2=`p_g0e1_not_d_fmse'
local tot2=`tot2'+`aff2'
local aff3=`p_g1e0_d'
local tot3=`p_g1e0_not_d_fmse'
local tot3=`tot3'+`aff3'
local aff4=`p_g1e1_d'
local tot4=`p_g1e1_not_d_fmse'
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*1000000000, 1)
replace tot=round(tot*1000000000, 1)
blogit aff tot g e gei, or
matrix A=e(b)
local flmse_or=exp(el(A, 1, 3))
lrtest, saving(0)
blogit aff tot g e, or
lrtest
local lrtest=r(chi2)
local ncp=`lrtest'*`ssize'/10^9
local flmse_power=(1-nchi2(1, `ncp', `chi_sl'))*100
local rec_ncp=npnchi2(1, `chi_sl', (1-`power'/100))
local flmse_ss=10^9*(`rec_ncp'/`lrtest')

local aff1=`p_g0e0_d'
local tot1=`p_g0e0_not_d_fmsge'
local tot1=`tot1'+`aff1'
local aff2=`p_g0e1_d'
local tot2=`p_g0e1_not_d_fmsge'
local tot2=`tot2'+`aff2'
local aff3=`p_g1e0_d'
local tot3=`p_g1e0_not_d_fmsge'
local tot3=`tot3'+`aff3'
local aff4=`p_g1e1_d'
local tot4=`p_g1e1_not_d_fmsge'
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*1000000000, 1)
replace tot=round(tot*1000000000, 1)
blogit aff tot g e gei, or
matrix A=e(b)
local flmsge_or=exp(el(A, 1, 3))
lrtest, saving(0)
blogit aff tot g e, or
lrtest
local lrtest=r(chi2)
local ncp=`lrtest'*`ssize'/10^9
local flmsge_power=(1-nchi2(1, `ncp', `chi_sl'))*100
local rec_ncp=npnchi2(1, `chi_sl', (1-`power'/100))
local flmsge_ss=10^9*(`rec_ncp'/`lrtest')

drop _all
use `using'
replace cmge_power=`cmge_power' in `loop2'
replace cmge_ss=`cmge_ss' in `loop2'
replace cmge_or=`cmge_or' in `loop2'
replace cme_power=`cme_power' in `loop2'
replace cme_ss=`cme_ss' in `loop2'
replace cme_or=`cme_or' in `loop2'
replace flme_power=`flme_power' in `loop2'
replace flme_ss=`flme_ss' in `loop2'
replace flme_or=`flme_or' in `loop2'
replace flmse_power=`flmse_power' in `loop2'
replace flmse_ss=`flmse_ss' in `loop2'
replace flmse_or=`flmse_or' in `loop2'
replace flmsge_power=`flmsge_power' in `loop2'
replace flmsge_ss=`flmsge_ss' in `loop2'
replace flmsge_or=`flmsge_or' in `loop2'
save `using', replace

}




replace cmge_power =round(cmge_power, 0.01)
replace cmge_ss=round(cmge_ss, 1)
replace cmge_or=round(cmge_or, 0.01)
replace cme_power =round(cme_power, 0.01)
replace cme_ss=round(cme_ss, 1)
replace cme_or=round(cme_or, 0.01)
replace flme_power =round(flme_power, 0.01)
replace flme_ss=round(flme_ss, 1)
replace flme_or=round(flme_or, 0.01)
replace flmse_power =round(flmse_power, 0.01)
replace flmse_ss=round(flmse_ss, 1)
replace flmse_or=round(flmse_or, 0.01)
replace flmsge_power =round(flmsge_power, 0.01)
replace flmsge_ss=round(flmsge_ss, 1)
replace flmsge_or=round(flmsge_or, 0.01)

replace cmge_power =. if pop_error==1
replace cmge_ss=. if pop_error==1
replace cmge_or=. if pop_error==1
replace cme_power =. if pop_error==1
replace cme_ss=. if pop_error==1
replace cme_or=. if pop_error==1
replace flme_power=. if pop_error==1|match_error==1
replace flme_ss=. if pop_error==1|match_error==1
replace flme_or=. if pop_error==1|match_error==1
replace flmse_power=. if pop_error==1|match_error==1
replace flmse_ss=. if pop_error==1|match_error==1
replace flmse_or=. if pop_error==1|match_error==1
replace flmsge_power=. if pop_error==1|match_error==1
replace flmsge_ss=. if pop_error==1|match_error==1
replace flmsge_or=. if pop_error==1|match_error==1

keep  pg pe rrg1e0 rrg0e1 rrint ssize power alpha_1 pd spece sense specg sensg orge dm cmge_power cmge_ss cmge_or cme_power cme_ss cme_or flme_power flme_ss flme_or flmse_power flmse_ss flmse_or flmsge_power flmsge_ss flmsge_or pop_error match_error


count
if r(N)~=1 {
noisily {
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
}
}

else {
foreach temp of varlist _all {
local `temp'=`temp'
}

local alpha_1=round(`alpha_1', 0.01)
noisily {

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 in gr "population odds ratio for the association of risk factors (orge):" in ye _col(70) orge
di in gr "degree of matching, flexible matching (dm):" in ye _col(70) dm
di in gr "specificity of the surrogate of the environmental exposure (spece):" in ye _col(70) spece
di in gr "sensitivity of the surrogate of the environmental exposure (sense):" in ye _col(70) sense
di in gr "specificity of the surrogate for the genotype (specg):" in ye _col(70) specg
di in gr "sensitivity of the surrogate for the genotype (sensg):" in ye _col(70) sensg
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 " in ye `ssize' in gr " cases"
di in gr "with a two-sided significance level (alpha=`alpha_1'):
di
di in ye `ssize'*3 in gr " controls (by design)"
di in gr "counter-matching on a surrogate of genotype and exposure:" in ye _col(70) cmge_power 
display in gr "counter-matching on a surrogate of exposure:" in ye _col(70) cme_power 

di
di in ye `ssize' in gr " controls"
display in gr "flexible matching on exposure frequency:" in ye _col(70) flme_power 
display in gr "flexible matching on exposure surrogate frequency:" in ye _col(70) flmse_power 
display in gr "flexible matching on exposure and genotype surrogate frequencies:" in ye _col(70) flmsge_power
di
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 "counter-matching on a surrogate of genotype and exposure:" in ye _col(70) cmge_or 
display in gr "counter-matching on a surrogate of exposure:" in ye _col(70) cme_or 
display in gr "flexible matching on exposure frequency:" in ye _col(70) flme_or 
display in gr "flexible matching on exposure surrogate frequency:" in ye _col(70) flmse_or 
display in gr "flexible matching on exposure and genotype surrogate frequencies:" in ye _col(70) flmsge_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:
di
di in gr "with " in ye 3 in gr " controls per case"
di in gr "counter-matching on a surrogate of genotype and exposure:" in ye _col(70) cmge_ss 
display in gr "counter-matching on a surrogate of exposure:" in ye _col(70) cme_ss 
di 
di in gr "with " in ye 1 in gr " control per case"
display in gr "flexible matching on exposure frequency:" in ye _col(70) flme_ss 
display in gr "flexible matching on exposure surrogate frequency:" in ye _col(70) flmse_ss 
display in gr "flexible matching on exposure and genotype surrogate frequencies:" in ye _col(70) flmsge_ss
di
di in gr "-------------------------------------------------------------------------------
di in gr "NOTES
di in gr "-------------------------------------------------------------------------------

di in gr "The power, required sample sizes and calculated interaction odds ratio for 
di in gr "these designs have been saved into file - " in ye "`using'"
di in gr "Type d for details of the output variables
}
}



sum pop_error
if r(mean)~=0 & r(N)==1 {
noi display
noi display in gr "All results are missing because the parameters specified in " in ye "`using'"
noi display in gr "lead to disease probabilities greater than one;"
noi display in gr "try specifying a lower disease frequency, or smaller effects."
}
else if r(mean)~=0 & r(N)~=1 {
noi display
noi di in gr "One or more or the sets of parameters (observations) that you entered, resulted
noi di in gr "in disease probabilities being greater than one in some or all situations.  For
noi di in gr "these observations the results are posted as missing and variable pop_error 
noi di in gr "(population error) will equal 1 rather than 0
}
else {
drop pop_error
}

sum match_error
if r(mean)~=0 {
noi display
noi display "The degree of matching (flexible matching designs) for one or more observations"
noi display "resulted in exposure/genotype frequencies of greater than 1 in some situations"
noi display "For these observations the results of the variable match_error will equal 1."
noi display "Try a degree of matching between 0 and 1"
}
else {
drop match_error
}


save `using', replace

}


end