*! (SJ3-1: st0032)

capture program drop second_primary
program define second_primary

version 7
set more off

quietly {

syntax using/

use `using'
count

*no. of ranges for which power is being sought
local scenar=r(N)

capture gen orge=1
keep pg pe rrg0e1 rrg1e0 rrint orge pd ssize power alpha_1

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|orge<0
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
}
}


gen rrg0e0=1
gen cc_u=.
gen pri_2_cas=.
gen pri_2_con=.

gen or_cc_u=.
gen or_pri_2_cas=.
gen or_pri_2_con=.

gen ss_cc_u=.
gen ss_pri_2_cas=.
gen ss_pri_2_con=.

set type double

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 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 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 p_d_2nd_pri=p_g0e0_d*p_d_g0e0+p_g0e1_d*p_d_g0e1+p_g1e0_d*p_d_g1e0+p_g1e1_d*p_d_g1e1

gen p_g0e0_2nd_pri_d= p_g0e0_d*  p_d_g0e0/ p_d_2nd_pri
gen p_g0e1_2nd_pri_d= p_g0e1_d*  p_d_g0e1/ p_d_2nd_pri
gen p_g1e0_2nd_pri_d= p_g1e0_d*  p_d_g1e0/ p_d_2nd_pri
gen p_g1e1_2nd_pri_d= p_g1e1_d*  p_d_g1e1/ p_d_2nd_pri

gen p_g0e0_2nd_pri_not_d= p_g0e0_d*(1- p_d_g0e0)/(1- p_d_2nd_pri)
gen p_g0e1_2nd_pri_not_d= p_g0e1_d*(1- p_d_g0e1)/(1- p_d_2nd_pri)
gen p_g1e0_2nd_pri_not_d= p_g1e0_d*(1- p_d_g1e0)/(1- p_d_2nd_pri)
gen p_g1e1_2nd_pri_not_d= p_g1e1_d*(1- p_d_g1e1)/(1- p_d_2nd_pri)

gen p_d_g0e02=p_d_2nd_pr/(p_g0e0_d* rrg0e0+ p_g0e1_d* rrg0e1+ p_g1e0_d* rrg1e0+ rrg1e0* rrg0e1* rrint* p_g1e1_d)
gen p_d_g0e12= p_d_g0e02* rrg0e1
gen p_d_g1e02= p_d_g0e02* rrg1e0
gen p_d_g1e12= p_d_g0e02* 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 
gen snd_pri_error=0
replace snd_pri_error=1 if  p_d_2nd_pri<0 | p_d_2nd_pri>1 |p_d_g0e02>1 | p_d_g0e12>1 | p_d_g1e02>1 | p_d_g1e12>1 | p_d_g0e02<0 | p_d_g0e12<0 | p_d_g1e02<0 | p_d_g1e12<0 

drop  p_d_2nd_pri p_d_g1e1 p_d_g1e0 p_d_g0e1 p_d_g0e0 p_g1e1 p_g1e0 p_g0e1 p_g0e0 p_d_g1e12 p_d_g1e02 p_d_g0e12 p_d_g0e02 

save `using', replace

local loop3=0

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

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

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'

local aff1= p_g0e0_2nd_pri_d in `loop3'
local tot1= p_g0e0_2nd_pri_not_d in `loop3'
local tot1=`tot1'+`aff1'

local aff2= p_g0e1_2nd_pri_d in `loop3'
local tot2= p_g0e1_2nd_pri_not_d in `loop3'
local tot2=`tot2'+`aff2'

local aff3= p_g1e0_2nd_pri_d in `loop3'
local tot3= p_g1e0_2nd_pri_not_d in `loop3'
local tot3=`tot3'+`aff3'

local aff4= p_g1e1_2nd_pri_d in `loop3'
local tot4= p_g1e1_2nd_pri_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 pri_2_cas= 1-nchi(1, `temp7', `alpha_1')

local ss_pri_2_cas=npnchi2(1, `alpha_1', (1-`power'))
local ss_pri_2_cas=`ss_pri_2_cas'*(`ssize')/`temp7'

matrix A=e(b)
local or_pri_2_cas=exp(el(A, 1, 3))

drop _all

use `using'

local aff1= p_g0e0_2nd_pri_d in `loop3'
local tot1= p_g0e0_not_d in `loop3'
local tot1=`tot1'+`aff1'

local aff2= p_g0e1_2nd_pri_d in `loop3'
local tot2= p_g0e1_not_d in `loop3'
local tot2=`tot2'+`aff2'

local aff3= p_g1e0_2nd_pri_d in `loop3'
local tot3= p_g1e0_not_d in `loop3'
local tot3=`tot3'+`aff3'

local aff4= p_g1e1_2nd_pri_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 pri_2_con= 1-nchi(1, `temp7', `alpha_1')

local ss_pri_2_con=npnchi2(1, `alpha_1', (1-`power'))
local ss_pri_2_con=`ss_pri_2_con'*(`ssize')/`temp7'

matrix A=e(b)
local or_pri_2_con=exp(el(A, 1, 3))

drop _all

use `using'

replace cc_u=`cc_u' in `loop3'
replace pri_2_cas=`pri_2_cas' in `loop3'
replace pri_2_con=`pri_2_con' in `loop3'

replace or_cc_u=`or_cc_u' in `loop3'
replace or_pri_2_cas=`or_pri_2_cas' in `loop3'
replace or_pri_2_con=`or_pri_2_con' in `loop3'

replace ss_cc_u=`ss_cc_u' in `loop3'
replace ss_pri_2_cas=`ss_pri_2_cas' in `loop3'
replace ss_pri_2_con=`ss_pri_2_con' in `loop3'

save `using', replace

}

drop  p_g0e0_d p_g0e1_d p_g1e0_d p_g1e1_d p_g0e0_not_d p_g0e1_not_d p_g1e0_not_d p_g1e1_not_d p_g0e0_2nd_pri_d p_g0e1_2nd_pri_d p_g1e0_2nd_pri_d p_g1e1_2nd_pri_d p_g0e0_2nd_pri_not_d p_g0e1_2nd_pri_not_d p_g1e0_2nd_pri_not_d p_g1e1_2nd_pri_not_d
save `using', replace

label var pg "P(G)"
label var pe "P(E)"
label var rrg0e0 "RR(G0E0)"
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 cc_u "Power (%) unmatched case-control study"
label var pri_2_cas "Power (%) 2 primaries cf 1 primary"
label var pri_2_con "Power (%) 2 primaries cf popn control"
label var or_cc_u "Interaction odds ratio, population based"
label var or_pri_2_cas "Interaction odds ratio, 2 primaries cf 1 primary"
label var or_pri_2_con "Interaction odds ratio, 2 primaries cf popn control"
label var ss_cc_u "sample size, population based case control study"
label var ss_pri_2_cas "sample size, 2 primaries cf 1 primary"
label var ss_pri_2_con "sample size, 2 primaries cf popn control"

replace cc_u=round(cc_u*100, 0.01)
replace pri_2_cas=round(pri_2_cas*100, 0.01)
replace pri_2_con=round(pri_2_con*100, 0.01)

replace or_cc_u=round(or_cc_u, 0.01)
replace or_pri_2_cas=round(or_pri_2_cas, 0.01)
replace or_pri_2_con=round(or_pri_2_con, 0.01)

replace ss_cc_u=round(ss_cc_u, 1)
replace ss_pri_2_cas=round(ss_pri_2_cas, 1)
replace ss_pri_2_con=round(ss_pri_2_con, 1)

drop rrg0e0

replace cc_u=. if pop_error==1
replace pri_2_cas=. if pop_error==1 | snd_pri_error==1
replace pri_2_con=. if pop_error==1 | snd_pri_error==1
replace or_cc_u=. if pop_error==1
replace or_pri_2_cas=. if pop_error==1 | snd_pri_error==1
replace or_pri_2_con=. if pop_error==1 | snd_pri_error==1
replace ss_cc_u=. if pop_error==1
replace ss_pri_2_cas=. if pop_error==1 | snd_pri_error==1
replace ss_pri_2_con=. if pop_error==1 | snd_pri_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 design with " in ye `ssize' in gr " cases and " in ye `ssize' in gr " controls" 
display in gr "population based case-control study, power (%):" in ye _col(70)  cc_u
display in gr "comparing people with 2 primary cancers with cases, power (%):" in ye _col(70)  pri_2_cas
display in gr "comparing people with 2 primary cancers with controls, power (%):" in ye _col(70)  pri_2_con

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 "population based case-control study:"  in ye _col(70) or_cc_u
display in gr "comparing people with 2 primary cancers with cases:" in ye _col(70)  or_pri_2_cas
display in gr "comparing people with 2 primary cancers with controls:" in ye _col(70)  or_pri_2_con
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)
display in gr  "population based case-control study:"  in ye _col(70)  ss_cc_u
display  in gr "comparing people with 2 primary cancers with cases:"  in ye _col(70)  ss_pri_2_cas
display  in gr "comparing people with 2 primary cancers with controls:"  in ye _col(70)  ss_pri_2_con

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 "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 & snd_pri_error==1 {
display
display "No results are given for the 2nd primary designs because the parameters in"
display in ye "`using'" in gr "under the assumed multiplicative disease model lead to the probability of
display "a second primary cancer being greater than one."
display "try specifying a lower disease frequency or a smaller effects"

drop pop_error snd_pri_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 pop_error
}
if r(mean)~=0 {
display
di in gr "One or more or the sets of parameters (observations) 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 if pop_error==0 & snd_pri_error==1 {
display
di in gr "One or more or the sets of parameters (observations) entered, resulted
di in gr "in disease probabilities being greater than one for the second primary design
di in gr "under the assumed multiplicative model in some or all situations.  For
di in gr "these observations the results are posted as missing and variable snd_pri_error 
di in gr "(error in the second primary design) will equal 1 rather than 0
quietly {
drop pop_error
save `using', replace
}
}
else {
drop pop_error snd_pri_error

quietly {
save `using', replace
}




}

end