*! (SJ3-1: st0032)

capture program drop vars1
program define vars1

drop _all

use $using

keep pg pe rrg0e1 rrg1e0 rrint pd sibore rrsibm ssize power alpha_1

gen orge=1
gen rrg0e0=1
gen dom_m_sib=.
gen dom_u_sib=.
gen rec_m_sib=.
gen rec_u_sib=.
gen cc_u=.
gen dom_asp=.
gen dom_coh=.
gen rec_asp=.
gen rec_coh=.
gen or_dom_m_sib=.
gen or_dom_u_sib=.
gen or_rec_m_sib=.
gen or_rec_u_sib=.
gen or_cc_u=.
gen or_dom_asp=.
gen or_dom_coh=.
gen or_rec_asp=.
gen or_rec_coh=.
gen ss_dom_m_sib=.
gen ss_dom_u_sib=.
gen ss_rec_m_sib=.
gen ss_rec_u_sib=.
gen ss_cc_u=.
gen ss_dom_asp=.
gen ss_dom_coh=.
gen ss_rec_asp=.
gen ss_rec_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 a_sore= sibore-1
gen b_sore=-2*sibore*pe-1+2*pe
gen c_sore=sibore*pe^2
gen p_e1e1=pe^2
replace p_e1e1=(-b_sore-(b_sore^2-4*a_sore*c_sore)^0.5)/(2*a_sore) if sibore~=1
gen p_e1_e1=p_e1e1/pe
gen p_e0_e1=1-p_e1_e1
gen p_e1_e0=(pe-p_e1e1)/(1-pe)
gen p_e0_e0=1-p_e1_e0
gen af_dom=(-2+(4-4*pg)^0.5)/(-2)
gen af_rec=pg^0.5
gen rec_p_g1g1=af_rec^4+(af_rec^3)*(1-af_rec)+0.25*(af_rec^2)*(1-af_rec)^2
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
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
gen rec_p_g1_g1=rec_p_g1g1/pg
gen rec_p_g0_g1=1-rec_p_g1_g1
gen rec_p_g1_g0=(pg-rec_p_g1g1)/(1-pg)
gen rec_p_g0_g0=1-rec_p_g1_g0
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 rec_a1=rec_p_g0_g0*p_e0_e0*p_g0e0_d
gen rec_b1=rec_p_g0_g0*p_e1_e0*p_g0e0_d
gen rec_c1=rec_p_g1_g0*p_e0_e0*p_g0e0_d
gen rec_d1=rec_p_g1_g0*p_e1_e0*p_g0e0_d
gen rec_e1=rec_p_g0_g0*p_e0_e1*p_g0e1_d
gen rec_f1=rec_p_g0_g0*p_e1_e1*p_g0e1_d
gen rec_g1=rec_p_g1_g0*p_e0_e1*p_g0e1_d
gen rec_h1=rec_p_g1_g0*p_e1_e1*p_g0e1_d
gen rec_i1=rec_p_g0_g1*p_e0_e0*p_g1e0_d
gen rec_j1=rec_p_g0_g1*p_e1_e0*p_g1e0_d
gen rec_k1=rec_p_g1_g1*p_e0_e0*p_g1e0_d
gen rec_l1=rec_p_g1_g1*p_e1_e0*p_g1e0_d
gen rec_m1=rec_p_g0_g1*p_e0_e1*p_g1e1_d
gen rec_n1=rec_p_g0_g1*p_e1_e1*p_g1e1_d
gen rec_o1=rec_p_g1_g1*p_e0_e1*p_g1e1_d
gen rec_p1=rec_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 rec_p_g0e0_sib_of_case=rec_a1+rec_e1+rec_i1+rec_m1
gen rec_p_g0e1_sib_of_case=rec_b1+rec_f1+rec_j1+rec_n1
gen rec_p_g1e0_sib_of_case=rec_c1+rec_g1+rec_k1+rec_o1
gen rec_p_g1e1_sib_of_case=rec_d1+rec_h1+rec_l1+rec_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 rec_pd_soc=rec_p_g0e0_sib_of_case*p_d_g0e0*rrsibm+rec_p_g0e1_sib_of_case*p_d_g0e1*rrsibm+rec_p_g1e0_sib_of_case*p_d_g1e0*rrsibm+rec_p_g1e1_sib_of_case*p_d_g1e1*rrsibm
gen rec_p_g0e0_soc_d=rec_p_g0e0_sib_of_case*p_d_g0e0*rrsibm/rec_pd_soc
gen rec_p_g0e1_soc_d=rec_p_g0e1_sib_of_case*p_d_g0e1*rrsibm/rec_pd_soc
gen rec_p_g1e0_soc_d=rec_p_g1e0_sib_of_case*p_d_g1e0*rrsibm/rec_pd_soc
gen rec_p_g1e1_soc_d=rec_p_g1e1_sib_of_case*p_d_g1e1*rrsibm/rec_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 rec_p_g0e0_soc_nd=rec_p_g0e0_sib_of_case*(1-p_d_g0e0*rrsibm)/(1-rec_pd_soc)
gen rec_p_g0e1_soc_nd=rec_p_g0e1_sib_of_case*(1-p_d_g0e1*rrsibm)/(1-rec_pd_soc)
gen rec_p_g1e0_soc_nd=rec_p_g1e0_sib_of_case*(1-p_d_g1e0*rrsibm)/(1-rec_pd_soc)
gen rec_p_g1e1_soc_nd=rec_p_g1e1_sib_of_case*(1-p_d_g1e1*rrsibm)/(1-rec_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 rec_a=(1-rrsibm*p_d_g0e0)*rec_a1/(1-rec_pd_soc)
gen rec_b=(1-rrsibm*p_d_g0e1)*rec_b1/(1-rec_pd_soc)
gen rec_c=(1-rrsibm*p_d_g1e0)*rec_c1/(1-rec_pd_soc)
gen rec_d=(1-rrsibm*p_d_g1e1)*rec_d1/(1-rec_pd_soc)
gen rec_e=(1-rrsibm*p_d_g0e0)*rec_e1/(1-rec_pd_soc)
gen rec_f=(1-rrsibm*p_d_g0e1)*rec_f1/(1-rec_pd_soc)
gen rec_g=(1-rrsibm*p_d_g1e0)*rec_g1/(1-rec_pd_soc)
gen rec_h=(1-rrsibm*p_d_g1e1)*rec_h1/(1-rec_pd_soc)
gen rec_i=(1-rrsibm*p_d_g0e0)*rec_i1/(1-rec_pd_soc)
gen rec_j=(1-rrsibm*p_d_g0e1)*rec_j1/(1-rec_pd_soc)
gen rec_k=(1-rrsibm*p_d_g1e0)*rec_k1/(1-rec_pd_soc)
gen rec_l=(1-rrsibm*p_d_g1e1)*rec_l1/(1-rec_pd_soc)
gen rec_m=(1-rrsibm*p_d_g0e0)*rec_m1/(1-rec_pd_soc)
gen rec_n=(1-rrsibm*p_d_g0e1)*rec_n1/(1-rec_pd_soc)
gen rec_o=(1-rrsibm*p_d_g1e0)*rec_o1/(1-rec_pd_soc)
gen rec_p=(1-rrsibm*p_d_g1e1)*rec_p1/(1-rec_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 sibore rrsibm  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 rec_p_g0e0_soc_d rec_p_g0e1_soc_d rec_p_g1e0_soc_d rec_p_g1e1_soc_d dom_p_g0e0_soc_nd dom_p_g0e1_soc_nd dom_p_g1e0_soc_nd dom_p_g1e1_soc_nd rec_p_g0e0_soc_nd rec_p_g0e1_soc_nd rec_p_g1e0_soc_nd rec_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 rec_a rec_b rec_c rec_d rec_e rec_f rec_g rec_h rec_i rec_j rec_k rec_l rec_m rec_n rec_o rec_p dom_m_sib rec_m_sib dom_u_sib rec_u_sib cc_u dom_asp rec_asp dom_coh rec_coh or_dom_m_sib or_rec_m_sib or_dom_u_sib or_rec_u_sib or_cc_u or_dom_asp or_rec_asp or_dom_coh or_rec_coh ssize power ss_dom_m_sib ss_rec_m_sib ss_dom_u_sib ss_rec_u_sib ss_cc_u ss_dom_asp ss_rec_asp ss_dom_coh ss_rec_coh alpha_1 pop_error sib_error
save $using, replace

end