*! (SJ3-1: st0032)

capture program drop vars2
program define vars2

drop _all

use $using

keep pg1 pg2 inh1 inh2 rrg1 rrg2 rrint pd rrsibm ssize power alpha_1

rename pg1 pg
rename pg2 pe
rename rrg1 rrg1e0
rename rrg2 rrg0e1

gen orge=1
gen rrg0e0=1
gen m_sib=.
gen u_sib=.
gen cc_u=.
gen asp=.
gen coh=.
gen or_m_sib=.
gen or_u_sib=.
gen or_cc_u=.
gen or_asp=.
gen or_coh=.
gen ss_m_sib=.
gen ss_u_sib=.
gen ss_cc_u=.
gen ss_asp=.
gen ss_coh=.
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 af_dom=(-2+(4-4*pg)^0.5)/(-2)
replace af_dom=pg^0.5 if inh1=="r" | inh1=="R"

gen dom_p_g1g1=af_dom^4+4*(af_dom^3)*(1-af_dom)+2*(af_dom^2)*(1-af_dom)^2+4*(0.75^2)*(af_dom^2)*(1-af_dom)^2+af_dom*(1-af_dom)^3
replace dom_p_g1g1=af_dom^4+(af_dom^3)*(1-af_dom)+0.25*(af_dom^2)*(1-af_dom)^2 if inh1=="r" | inh1=="R"

gen dom_p_g1_g1=dom_p_g1g1/pg
gen dom_p_g0_g1=1-dom_p_g1_g1
gen dom_p_g1_g0=(pg-dom_p_g1g1)/(1-pg)
gen dom_p_g0_g0=1-dom_p_g1_g0

replace af_dom=(-2+(4-4*pe)^0.5)/(-2)
replace af_dom=pe^0.5 if inh2=="r" | inh2=="R"

replace dom_p_g1g1=af_dom^4+4*(af_dom^3)*(1-af_dom)+2*(af_dom^2)*(1-af_dom)^2+4*(0.75^2)*(af_dom^2)*(1-af_dom)^2+af_dom*(1-af_dom)^3
replace dom_p_g1g1=af_dom^4+(af_dom^3)*(1-af_dom)+0.25*(af_dom^2)*(1-af_dom)^2 if inh2=="r" | inh2=="R"

gen p_e1_e1=dom_p_g1g1/pe
gen p_e0_e1=1-p_e1_e1
gen p_e1_e0=(pe-dom_p_g1g1)/(1-pe)
gen p_e0_e0=1-p_e1_e0

gen dom_a1=dom_p_g0_g0*p_e0_e0*p_g0e0_d
gen dom_b1=dom_p_g0_g0*p_e1_e0*p_g0e0_d
gen dom_c1=dom_p_g1_g0*p_e0_e0*p_g0e0_d
gen dom_d1=dom_p_g1_g0*p_e1_e0*p_g0e0_d
gen dom_e1=dom_p_g0_g0*p_e0_e1*p_g0e1_d
gen dom_f1=dom_p_g0_g0*p_e1_e1*p_g0e1_d
gen dom_g1=dom_p_g1_g0*p_e0_e1*p_g0e1_d
gen dom_h1=dom_p_g1_g0*p_e1_e1*p_g0e1_d
gen dom_i1=dom_p_g0_g1*p_e0_e0*p_g1e0_d
gen dom_j1=dom_p_g0_g1*p_e1_e0*p_g1e0_d
gen dom_k1=dom_p_g1_g1*p_e0_e0*p_g1e0_d
gen dom_l1=dom_p_g1_g1*p_e1_e0*p_g1e0_d
gen dom_m1=dom_p_g0_g1*p_e0_e1*p_g1e1_d
gen dom_n1=dom_p_g0_g1*p_e1_e1*p_g1e1_d
gen dom_o1=dom_p_g1_g1*p_e0_e1*p_g1e1_d
gen dom_p1=dom_p_g1_g1*p_e1_e1*p_g1e1_d

gen dom_p_g0e0_sib_of_case=dom_a1+dom_e1+dom_i1+dom_m1
gen dom_p_g0e1_sib_of_case=dom_b1+dom_f1+dom_j1+dom_n1
gen dom_p_g1e0_sib_of_case=dom_c1+dom_g1+dom_k1+dom_o1
gen dom_p_g1e1_sib_of_case=dom_d1+dom_h1+dom_l1+dom_p1

gen dom_pd_soc=dom_p_g0e0_sib_of_case*p_d_g0e0*rrsibm+dom_p_g0e1_sib_of_case*p_d_g0e1*rrsibm+dom_p_g1e0_sib_of_case*p_d_g1e0*rrsibm+dom_p_g1e1_sib_of_case*p_d_g1e1*rrsibm

gen dom_p_g0e0_soc_d=dom_p_g0e0_sib_of_case*p_d_g0e0*rrsibm/dom_pd_soc
gen dom_p_g0e1_soc_d=dom_p_g0e1_sib_of_case*p_d_g0e1*rrsibm/dom_pd_soc
gen dom_p_g1e0_soc_d=dom_p_g1e0_sib_of_case*p_d_g1e0*rrsibm/dom_pd_soc
gen dom_p_g1e1_soc_d=dom_p_g1e1_sib_of_case*p_d_g1e1*rrsibm/dom_pd_soc

gen dom_p_g0e0_soc_nd=dom_p_g0e0_sib_of_case*(1-p_d_g0e0*rrsibm)/(1-dom_pd_soc)
gen dom_p_g0e1_soc_nd=dom_p_g0e1_sib_of_case*(1-p_d_g0e1*rrsibm)/(1-dom_pd_soc)
gen dom_p_g1e0_soc_nd=dom_p_g1e0_sib_of_case*(1-p_d_g1e0*rrsibm)/(1-dom_pd_soc)
gen dom_p_g1e1_soc_nd=dom_p_g1e1_sib_of_case*(1-p_d_g1e1*rrsibm)/(1-dom_pd_soc)

gen dom_a=(1-rrsibm*p_d_g0e0)*dom_a1/(1-dom_pd_soc)
gen dom_b=(1-rrsibm*p_d_g0e1)*dom_b1/(1-dom_pd_soc)
gen dom_c=(1-rrsibm*p_d_g1e0)*dom_c1/(1-dom_pd_soc)
gen dom_d=(1-rrsibm*p_d_g1e1)*dom_d1/(1-dom_pd_soc)
gen dom_e=(1-rrsibm*p_d_g0e0)*dom_e1/(1-dom_pd_soc)
gen dom_f=(1-rrsibm*p_d_g0e1)*dom_f1/(1-dom_pd_soc)
gen dom_g=(1-rrsibm*p_d_g1e0)*dom_g1/(1-dom_pd_soc)
gen dom_h=(1-rrsibm*p_d_g1e1)*dom_h1/(1-dom_pd_soc)
gen dom_i=(1-rrsibm*p_d_g0e0)*dom_i1/(1-dom_pd_soc)
gen dom_j=(1-rrsibm*p_d_g0e1)*dom_j1/(1-dom_pd_soc)
gen dom_k=(1-rrsibm*p_d_g1e0)*dom_k1/(1-dom_pd_soc)
gen dom_l=(1-rrsibm*p_d_g1e1)*dom_l1/(1-dom_pd_soc)
gen dom_m=(1-rrsibm*p_d_g0e0)*dom_m1/(1-dom_pd_soc)
gen dom_n=(1-rrsibm*p_d_g0e1)*dom_n1/(1-dom_pd_soc)
gen dom_o=(1-rrsibm*p_d_g1e0)*dom_o1/(1-dom_pd_soc)
gen dom_p=(1-rrsibm*p_d_g1e1)*dom_p1/(1-dom_pd_soc)

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 sib_error=0
replace sib_error=1 if p_d_g0e0*rrsibm>1 | p_d_g0e1*rrsibm>1 | p_d_g1e0*rrsibm>1 | p_d_g1e1*rrsibm>1 | p_d_g0e0*rrsibm<0 | p_d_g0e1*rrsibm<0 | p_d_g1e0*rrsibm<0 | p_d_g1e1*rrsibm<0
keep  pg pe rrg0e0 rrg1e0 rrg0e1 rrint orge pd rrsibm inh1 inh2 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 dom_p_g0e0_soc_d dom_p_g0e1_soc_d dom_p_g1e0_soc_d dom_p_g1e1_soc_d  dom_p_g0e0_soc_nd dom_p_g0e1_soc_nd dom_p_g1e0_soc_nd dom_p_g1e1_soc_nd  dom_a dom_b dom_c dom_d dom_e dom_f dom_g dom_h dom_i dom_j dom_k dom_l dom_n dom_m dom_o dom_p  m_sib u_sib cc_u asp coh or_m_sib or_u_sib or_cc_u or_asp  or_coh  ssize power ss_m_sib  ss_u_sib  ss_cc_u ss_asp  ss_coh  alpha_1 pop_error sib_error
save $using, replace

end