**************************************************************** *! Version 9.0.2, 7 October 2006 *! Author: James Cui, Monash University *! Perform Hardy-Weinberg Equilibrium Test in case-control study *! Original publication: Sep 2000 STB 57: 17-19 *! Update 07Oct06: table position use _col **************************************************************** capture program drop genhwcci program genhwcci, rclass version 9.0 gettoken AA1 0 : 0, parse(" ,") gettoken Aa1 0 : 0, parse(" ,") gettoken aa1 0 : 0, parse(" ,") gettoken AA2 0 : 0, parse(" ,") gettoken Aa2 0 : 0, parse(" ,") gettoken aa2 0 : 0, parse(" ,") confirm integer number `AA1' confirm integer number `Aa1' confirm integer number `aa1' confirm integer number `AA2' confirm integer number `Aa2' confirm integer number `aa2' if (`AA1'<0 | `Aa1'<0 | `aa1'<0 | `AA2'<0 | `Aa2'<0 | `aa2'<0) { di in red "negative numbers invalid" exit 498 } syntax [, Label(string) Binvar ] tempname obs1 obs2 obs AA Aa aa nA1 na1 nA2 na2 nA na scalar `obs1' = `AA1' + `Aa1' + `aa1' scalar `obs2' = `AA2' + `Aa2' + `aa2' scalar `obs' = `obs1' + `obs2' scalar `AA' = `AA1' + `AA2' scalar `Aa' = `Aa1' + `Aa2' scalar `aa' = `aa1' + `aa2' scalar `nA1' = 2*(`AA1') + `Aa1' scalar `na1' = 2*(`aa1') + `Aa1' scalar `nA2' = 2*(`AA2') + `Aa2' scalar `na2' = 2*(`aa2') + `Aa2' scalar `nA' = 2*(`AA') + `Aa' scalar `na' = 2*(`aa') + `Aa' local nob = `nA' + 1 if _N < `nA' { preserve qui set obs `nob' } if "`label'" ~= "" { parse "`label'", parse (" ") local naAA = "`1'" local naAa = "`2'" local naaa = "`3'" } else { local naAA = "AA" local naAa = "Aa" local naaa = "aa" } local al1 = substr("`naAA'",1,1) local al2 = substr("`naaa'",1,1) *--------------------------------------------------------------- * 1. CASE-CONTROL TABLE *--------------------------------------------------------------- di in gr _col(9) "Genotype |" _col(28) "Case" _col(41) "Control |" _col(59) "Total" di in gr _col(6) "------------+-------------------------------+-------------" di in gr %16s /* */ "`naAA'" " |" in ye _col(24) %8.0f `AA1' _col(40) %8.0f `AA2' /* */ _col(50) in gr "|" in ye _col(56) %8.0f `AA' di in gr %16s /* */ "`naAa'" " |" in ye _col(24) %8.0f `Aa1' _col(40) %8.0f `Aa2' /* */ _col(50) in gr "|" in ye _col(56) %8.0f `Aa' di in gr %16s /* */ "`naaa'" " |" in ye _col(24) %8.0f `aa1' _col(40) %8.0f `aa2' /* */ _col(50) in gr "|" in ye _col(56) %8.0f `aa' di in gr _col(6) "------------+-------------------------------+-------------" di in gr _col(12) "total" " |" /* */ in ye _col(24) %8.0f `obs1' _col(40) %8.0f `obs2' /* */ _col(50) in gr "|" in ye _col(56) %8.0f `obs' *--------------------------------------------------------------- * 2. CHI-SQUARE TEST *--------------------------------------------------------------- tempname pAA1 pAa1 paa1 pA1 pa1 pA2 pa2 scalar `pAA1' = `AA1'/`obs1' scalar `pAa1' = `Aa1'/`obs1' scalar `paa1' = `aa1'/`obs1' scalar `pA1' = `nA1'/(2*`obs1') scalar `pa1' = `na1'/(2*`obs1') scalar `pA2' = `nA2'/(2*`obs2') scalar `pa2' = `na2'/(2*`obs2') tempname eAA1 eAa1 eaa1 eAA2 eAa2 eaa2 xAA1 xAa1 xaa1 tempname xAA2 xAa2 xaa2 chi1 chi2 eobs1 eobs2 scalar `eAA1' = `obs1'*(`pA1'^2) scalar `eAa1' = `obs1'*2*(`pA1'*`pa1') scalar `eaa1' = `obs1'*(`pa1'^2) scalar `eAA2' = `obs2'*(`pA2'^2) scalar `eAa2' = `obs2'*2*(`pA2'*`pa2') scalar `eaa2' = `obs2'*(`pa2'^2) scalar `xAA1' = ((`AA1'-`eAA1')^2)/`eAA1' scalar `xAa1' = ((`Aa1'-`eAa1')^2)/`eAa1' scalar `xaa1' = ((`aa1'-`eaa1')^2)/`eaa1' scalar `xAA2' = ((`AA2'-`eAA2')^2)/`eAA2' scalar `xAa2' = ((`Aa2'-`eAa2')^2)/`eAa2' scalar `xaa2' = ((`aa2'-`eaa2')^2)/`eaa2' scalar `chi1' = `xAA1'+`xAa1'+`xaa1' scalar `chi2' = `xAA2'+`xAa2'+`xaa2' scalar `eobs1' = `eAA1'+`eAa1'+`eaa1' scalar `eobs2' = `eAA2'+`eAa2'+`eaa2' tempname llo1 lle1 llo2 lle2 lr1 lr2 scalar `llo1'= (`AA1'*ln(`AA1'/`obs1')) /* */ + (`Aa1'*ln(`Aa1'/`obs1')) /* */ + (`aa1'*ln(`aa1'/`obs1')) scalar `lle1'= (`AA1'*ln(`eAA1'/`eobs1')) /* */ + (`Aa1'*ln(`eAa1'/`eobs1')) /* */ + (`aa1'*ln(`eaa1'/`eobs1')) scalar `llo2'= (`AA2'*ln(`AA2'/`obs2')) /* */ + (`Aa2'*ln(`Aa2'/`obs2')) /* */ + (`aa2'*ln(`aa2'/`obs2')) scalar `lle2'= (`AA2'*ln(`eAA2'/`eobs2')) /* */ + (`Aa2'*ln(`eAa2'/`eobs2')) /* */ + (`aa2'*ln(`eaa2'/`eobs2')) scalar `lr1'=2*(`llo1'-`lle1') scalar `lr2'=2*(`llo2'-`lle2') tempvar nAA1 nAa1 naa1 pobs1 pcal1 pval1 sump1 qui gen int `nAA1'=_n-1 qui gen int `nAa1'=`nA1'-2*(`nAA1') qui replace `nAA1'=. if `nAa1'<0 qui gen double `pcal1'= /* */ lnfact(`obs1') - lnfact(2*`obs1') + lnfact(`nA1') + /* */ lnfact(2*`obs1'-`nA1') - lnfact(`nAa1')-lnfact((`nA1' - `nAa1')/2)- /* */ lnfact(`obs1'-(`nA1'+`nAa1')/2) qui replace `pcal1'=exp(`pcal1'+(`nAa1')*ln(2)) qui gen double `pobs1'=`pcal1' if `nAa1'==`Aa1' qui sort `pobs1' `pcal1' qui replace `pobs1'=`pobs1'[_n-1] if `pobs1'==. & `pcal1'!=. qui gen double `pval1'=`pcal1' if `pcal1'<=`pobs1' qui gen double `sump1'=sum(`pval1') tempvar nAA2 nAa2 naa2 pobs2 pcal2 pval2 sump2 qui gen int `nAA2'=_n-1 qui gen int `nAa2'=`nA2'-2*(`nAA2') qui replace `nAA2'=. if `nAa2'<0 qui gen double `pcal2'= /* */ lnfact(`obs2') - lnfact(2*`obs2') + lnfact(`nA2') + /* */ lnfact(2*`obs2'-`nA2') - lnfact(`nAa2')-lnfact((`nA2' - `nAa2')/2)- /* */ lnfact(`obs2'-(`nA2'+`nAa2')/2) qui replace `pcal2'=exp(`pcal2'+(`nAa2')*ln(2)) qui gen double `pobs2'=`pcal2' if `nAa2'==`Aa2' qui sort `pobs2' `pcal2' qui replace `pobs2'=`pobs2'[_n-1] if `pobs2'==. & `pcal2'!=. qui gen double `pval2'=`pcal2' if `pcal2'<=`pobs2' qui gen double `sump2'=sum(`pval2') tempname p0 q0 L0 L1 x2 scalar `p0'= `nA'/(2*`obs') scalar `q0'= 1.0-`p0' scalar `L0'= `AA'*ln((`p0')^2)+`Aa'*ln(2*(`p0')*(`q0'))+`aa'*ln((`q0')^2) scalar `L1'=`AA1'*ln(`pAA1')+`Aa1'*ln(`pAa1')+`aa1'*ln(`paa1') /* */ +`AA2'*ln((`pA2')^2)+`Aa2'*ln(2*`pA2'*`pa2')+ /* */ `aa2'*ln((`pa2')^2) qui gen `x2' = -2*(`L0'-`L1') tempname D1 D2 VD1 VD2 ED1 ED2 scalar `D1'=(`AA1'/`obs1')-(`pA1'^2) scalar `D2'=(`AA2'/`obs2')-(`pA2'^2) scalar `ED1'=`D1'-(1/(2*`obs1'))*(`pA1'*(1-`pA1')+`D1') scalar `ED2'=`D2'-(1/(2*`obs2'))*(`pA2'*(1-`pA2')+`D2') scalar `VD1'=(1/`obs1')* ( (`pA1'^2) * ((1-`pA1')^2) + ((1-2*`pA1')^2)*`D1' - `D1'^2) scalar `VD2'=(1/`obs2')* ( (`pA2'^2) * ((1-`pA2')^2) + ((1-2*`pA2')^2)*`D2' - `D2'^2) tempname se1 se2 se3 se4 if `"`binvar'"'=="" { scalar `se1'= sqrt((1/(2*`obs1'))*(`pA1'+(`AA1'/`obs1')-2*(`pA1'^2))) scalar `se2'= sqrt((1/(2*`obs1'))*(`pa1'+(`aa1'/`obs1')-2*(`pa1'^2))) scalar `se3'= sqrt((1/(2*`obs2'))*(`pA2'+(`AA2'/`obs2')-2*(`pA2'^2))) scalar `se4'= sqrt((1/(2*`obs2'))*(`pa2'+(`aa2'/`obs2')-2*(`pa2'^2))) } else { scalar `se1'=sqrt(`pA1' * `pa1' / (2 * `obs1')) scalar `se2'=`se1' scalar `se3'=sqrt(`pA2' * `pa2' / (2 * `obs2')) scalar `se4'=`se3' local vartype " (binomial)" } *--------------------------------------------------------------- * 3. CASE ONLY TABLE *--------------------------------------------------------------- di "" di in gr _col(13) "Case" /* di in gr " Allele frequencies " */ di in gr /* */ " Allele | Case Frequency Std. Err." di in gr _col(6) "------------+--------------------------------------" di in gr /* */ %16s "`al1'" " |" in ye _col(20) %8.0f `nA1' /* */ _col(34) %8.4f `pA1' _col(49) %8.4f `se1' "`vartype'" di in gr /* */ %16s "`al2'" " |" in ye _col(20) %8.0f `na1'/* */ _col(34) %8.4f `pa1' _col(49) %8.4f `se2' "`vartype'" di in gr _col(6) "------------+--------------------------------------" di in gr _col(12) "total" " |" /* */ in ye _col(20) %8.0f `na1'+`nA1' _col(34) %8.4f `pA1'+`pa1' di in gr _col(6) "------------+--------------------------------------" di in gr _col(6) "Estimated disequilibrium coefficient (D) = " /* */ in ye %8.4f `D1' /* */ _n in gr _col(42) " SE = " in ye %8.4f sqrt(`VD1') di in gr " Hardy-Weinberg Equilibrium Test: " di in gr /* */ " Pearson chi2 (" in ye "1" in gr ") = " /* */ in ye %8.3f `chi1' /* */ in gr " Pr= " in ye %5.4f chiprob(1,`chi1') di in gr /* */ " likelihood-ratio chi2 (" in ye "1" in gr ") = " /* */ in ye %8.3f `lr1' /* */ in gr " Pr= " in ye %5.4f chiprob(1,`lr1') di in gr /* */ " Exact significance prob = " /* */ in ye %5.4f `sump1'[_N] *--------------------------------------------------------------- * 4. RETURN LIST IN CASES *--------------------------------------------------------------- return scalar N_1 = `obs1' return scalar chi2_1 = `chi1' return scalar pchi2_1 = chiprob(1,`chi1') return scalar lr_1 = `lr1' return scalar plr_1 = chiprob(1,`lr1') return scalar pexact_1 = `sump1'[_N] *--------------------------------------------------------------- * 5. CONTROL ONLY TABLE *--------------------------------------------------------------- di "" di in gr _col(10) "Control" /* di in gr " Allele frequencies " */ di in gr /* */ " Allele | Control Frequency Std. Err." di in gr _col(6) "------------+--------------------------------------" di in gr /* */ %16s "`al1'" " |" in ye _col(20) %8.0f `nA2' /* */ _col(34) %8.4f `pA2' _col(49) %8.4f `se3' "`vartype'" di in gr /* */ %16s "`al2'" " |" in ye _col(20) %8.0f `na2'/* */ _col(34) %8.4f `pa2' _col(49) %8.4f `se4' "`vartype'" di in gr _col(6) "------------+--------------------------------------" di in gr _col(12) "total" " |" /* */ in ye _col(20) %8.0f `na2'+`nA2' _col(34) %8.4f `pA2'+`pa2' di in gr _col(6) "------------+--------------------------------------" di in gr _col(6) "Estimated disequilibrium coefficient (D) = " /* */ in ye %8.4f `D2' /* */ _n in gr _col(42) " SE = " in ye %8.4f sqrt(`VD2') di in gr " Hardy-Weinberg Equilibrium Test: " di in gr /* */ " Pearson chi2 (" in ye "1" in gr ") = " /* */ in ye %8.3f `chi2' /* */ in gr " Pr= " in ye %5.4f chiprob(1,`chi2') di in gr /* */ " likelihood-ratio chi2 (" in ye "1" in gr ") = " /* */ in ye %8.3f `lr2' /* */ in gr " Pr= " in ye %5.4f chiprob(1,`lr2') di in gr /* */ " Exact significance prob = " /* */ in ye %5.4f `sump2'[_N] *--------------------------------------------------------------- * 6. RETURN LIST IN CONTROLS *--------------------------------------------------------------- return scalar N_0 = `obs2' return scalar chi2_0 = `chi2' return scalar pchi2_0 = chiprob(1,`chi2') return scalar lr_0 = `lr2' return scalar plr_0 = chiprob(1,`lr2') return scalar pexact_0 = `sump2'[_N] *--------------------------------------------------------------- * 7. GIVEN CONTROLS UNDER HWE, TEST CASES UNDER HWE *--------------------------------------------------------------- di _n _col(6) in gr "Test H0: cases under HWE: (given controls under HWE)" di in gr /* */ " likelihood-ratio chi2 (" in ye "2" in gr ") = " /* */ in ye %8.3f `x2' /* */ in gr " Pr= " in ye %5.4f chiprob(2,`x2') *--------------------------------------------------------------- * 8. UPDATED RETURN VALUES *--------------------------------------------------------------- return scalar lr_c = `x2' return scalar plr_c = chiprob(2,`x2') end