*! (SJ3-1: st0032)

capture program drop co_power
program define co_power

quietly {

clear

syntax using/

global using=`"`using'"'

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

keep pg pe rrg0e1 rrg1e0 rrint pd ssize power alpha_1 orge

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
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 rrg0e0=1

gen ss_co=.
gen co_power=.
gen or_co=.

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

save `using', replace


drop _all
local loop=0


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

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

use `using'

local g0e0=p_g0e0_d in `loop'
local g0e1=p_g0e1_d in `loop'
local g1e0=p_g1e0_d in `loop'
local g1e1=p_g1e1_d in `loop'
local power=power in `loop'
local sample_size=ssize in `loop'
local alpha_1=alpha_1 in `loop' 

drop _all

set obs 2
gen g=_n-1
gen no_g=`g0e0'
replace no_g=`g1e0' in 2
gen no_g_and_e=`g0e1'
replace no_g_and_e=`g1e1' in 2
replace no_g=no_g+no_g_and_e

replace no_g=round(1000000*no_g, 1)
replace no_g_and_e=round(1000000*no_g_and_e, 1)

blogit no_g_and_e no_g g, or
matrix a=e(b)
local or_co=exp(a[1,1])
lrtest, saving(0)
blogit no_g_and_e no_g, or
lrtest


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

local temp7=`sample_size'*r(chi2)/1000000
local co= 1-nchi(1, `temp7', `alpha_1')
local co=`co'*100
local ss_co=npnchi2(1, `alpha_1', (1-`power'))
local ss_co=`ss_co'*`sample_size'/`temp7'

displ `co'

drop _all
use `using'
replace ss_co=`ss_co' in `loop'
replace co_power=`co' in `loop'
replace or_co=`or_co' in `loop'
save `using', replace

}

keep pg pe rrg0e1 rrg1e0 rrint pd power ssize alpha_1 alpha_1 orge co_power ss_co or_co pop_error

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 co_power "Power (%) case only study"
label var or_co "Interaction odds ratio, case only study"
label var ss_co "sample size, case only study"

replace co_power=round(co_power, 0.01)
replace or_co=round(or_co, 0.01)
replace ss_co=round(ss_co,1)

replace or_co=. if pop_error==1
replace ss_co=. if pop_error==1
replace co_power=. 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)
local co_power=round(`co_power',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 in gr "population odds ratio for the association of risk factors (orge):" in ye _col(70) orge
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
di in gr "Case-only design:" in ye _col(50) `co_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 "Case-only design:" in ye _col(50) or_co
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 (= number of controls):
di
di in gr "Case-only design:" in ye _col(50) `ss_co'
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 the 
di in gr "case-only 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
}

quietly {
sum orge
}
if r(mean)~=1 {
display
display "A value for the population association of risk factors was specified that"
display "implies non-independence (orge~=1)."
display "Be aware that the false positive rate of the case-only design in these
di "situations is inflated."
}

}
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-only 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
}

quietly {
sum orge
}
if r(mean)~=1 {
display
display "One or more values for the population association of risk factors were specified"
display "such that implies there was non-independence (orge~=1)."
display "Be aware that the false positive rate of the case-only design in these
di "situations is inflated."
}

}

quietly {
sum orge
}
if r(mean)~=1 {
di
display "One or more or the sets of parameters (observations) that you entered specified
di "values for the population association of risk factors that implies
di "non-independence
display "Be aware that the false positive rate of the case only design in these 
di "situations is inflated."
}


end