#EM approach to haplotype frequency estimation

library(haplo.stats)
fmsURL <- "http://www.stat-gen.org/book.e1/data/FMS_data.txt"
fms <- read.delim(file=fmsURL, header=T, sep="\t")
attach(fms)

Geno <- cbind(substr(actn3_r577x,1,1), substr(actn3_r577x,2,2),
substr(actn3_rs540874,1,1), substr(actn3_rs540874,2,2),
substr(actn3_rs1815739,1,1), substr(actn3_rs1815739,2,2),
substr(actn3_1671064,1,1), substr(actn3_1671064,2,2))
SNPnames <- c("actn3_r577x", "actn3_rs540874", "actn3_rs1815739", "actn3_1671064")

Geno.C <- Geno[Race=="Caucasian" & !is.na(Race),]
HaploEM <- haplo.em(Geno.C, locus.label=SNPnames, 
                    control=haplo.em.control(min.posterior=1e-4))
HaploEM
## =========================================================================== 
##                                 Haplotypes                                  
## =========================================================================== 
##    actn3_r577x actn3_rs540874 actn3_rs1815739 actn3_1671064 hap.freq
## 1            C              A               C             G  0.00261
## 2            C              A               T             A  0.00934
## 3            C              A               T             G  0.01354
## 4            C              G               C             A  0.47294
## 5            C              G               C             G  0.01059
## 6            T              A               C             A  0.00065
## 7            T              A               T             G  0.39891
## 8            T              G               C             A  0.08557
## 9            T              G               T             A  0.00065
## 10           T              G               T             G  0.00520
## =========================================================================== 
##                                   Details                                   
## =========================================================================== 
## lnlike =  -1285.406 
## lr stat for no LD =  2780.769 , df =  5 , p-val =  0
Geno.AA <- Geno[Race=="African Am" & !is.na(Race),]
HaploEM2 <- haplo.em(Geno.AA, locus.label=SNPnames, 
                     control=haplo.em.control(min.posterior=1e-4))
HaploEM2
## =========================================================================== 
##                                 Haplotypes                                  
## =========================================================================== 
##   actn3_r577x actn3_rs540874 actn3_rs1815739 actn3_1671064 hap.freq
## 1           C              A               C             A  0.01100
## 2           C              A               C             G  0.08131
## 3           C              A               T             G  0.03764
## 4           C              G               C             A  0.57763
## 5           C              G               C             G  0.01104
## 6           T              A               C             A  0.00072
## 7           T              A               T             G  0.17167
## 8           T              G               C             A  0.10833
## 9           T              G               C             G  0.00068
## =========================================================================== 
##                                   Details                                   
## =========================================================================== 
## lnlike =  -84.97943 
## lr stat for no LD =  119.5876 , df =  4 , p-val =  0
#Calculating posterior haplotype probabilities
HaploEM$nreps[1:5]
## indx.subj
## 1 2 3 4 5 
## 1 2 2 2 1
HaploEM$indx.subj[1:8]
## [1] 1 2 2 3 3 4 4 5
HaploEM$hap1code[1:8]
## [1] 4 3 7 8 7 8 7 4
HaploEM$hap2code[1:8]
## [1] 4 8 4 3 4 3 4 4
HaploEM$post[1:8]
## [1] 1.000000000 0.006102808 0.993897192 0.006102808 0.993897192 0.006102808
## [7] 0.993897192 1.000000000
HapProb <- HaploEM$hap.prob
HapProb
##  [1] 0.0026138447 0.0093400121 0.0135382726 0.4729357032 0.0105890282
##  [6] 0.0006518550 0.3989126970 0.0855667219 0.0006548104 0.0051970549
p1 <- 2*prod(HapProb[c(3,8)])
p2 <- 2*prod(HapProb[c(4,7)])
p1 / (p1+p2)
## [1] 0.006102807
p2 / (p1+p2)
## [1] 0.9938972
hgdpURL <- "http://www.stat-gen.org/book.e1/data/HGDP_AKT1.txt"
hgdp <- read.delim(file=hgdpURL, header=T, sep="\t")

virco <- read.csv(file="Virco_data.csv", header=T, sep=",")

#Testing hypotheses about haplotype frequencies within the EM framework
HapDesign <- function(HaploEM){
Nobs <- length(unique(HaploEM$indx.subj)) #number of observations
Nhap <- length(HaploEM$hap.prob) #number of haplotypes
XmatHap <- matrix(data=0,nrow=Nobs,ncol=Nhap)
for (i in 1:Nobs){
IDSeq <- seq(1:sum(HaploEM$nreps))[HaploEM$indx.subj==i]
for (j in 1:length(IDSeq)){
XmatHap[i,HaploEM$hap1code[IDSeq][j]] <-XmatHap[i,HaploEM$hap1code[IDSeq][j]] +
HaploEM$post[IDSeq][j]
XmatHap[i,HaploEM$hap2code[IDSeq][j]] <-
XmatHap[i,HaploEM$hap2code[IDSeq][j]] +
HaploEM$post[IDSeq][j]
}
}
return(XmatHap)
}


HapFreqSE <- function(HaploEM){
HapMat <- HapDesign(HaploEM)
Nobs <- length(unique(HaploEM$indx.subj)) #number of observations
Nhap <- length(HaploEM$hap.prob) #number of haplotypes
S.Full<-matrix(data=0, nrow=Nobs, ncol=Nhap-1)
for(i in 1:Nobs){
for(k in 1:(Nhap-1)){
S.Full[i,k]<-HapMat[i,k]/HaploEM$hap.prob[k]-
HapMat[i,Nhap]/HaploEM$hap.prob[Nhap]
}
}
Score<-t(S.Full)%*%S.Full
invScore<-solve(Score)
HapSE<-c(sqrt(diag(invScore)),
sqrt(t(rep(1,Nhap-1))%*%invScore%*%rep(1,Nhap-1)))
return(HapSE)
}

FreqDiff <- HaploEM2$hap.prob[4] - HaploEM$hap.prob[4]
s1 <- HapFreqSE(HaploEM)[4]
s2 <- HapFreqSE(HaploEM2)[4]
SE <- sqrt(s1^2 + s2^2)
CI <- c(FreqDiff - 1.96*SE, FreqDiff + 1.96*SE)
CI
## [1] -0.003390307  0.212779484
#Application of haplotype trend regression
attach(fms)
## The following objects are masked from fms (pos = 3):
## 
##     acdc_rs1501299, ace_id, actn3_1671064, actn3_r577x,
##     actn3_rs1815739, actn3_rs540874, adrb2_1042713, adrb2_1042714,
##     adrb2_rs1042718, adrb3_4994, Age, agrp_5030980, akt1_a15756t,
##     akt1_a22889g, akt1_a7699g, akt1_c10744t_c12886t, akt1_c15676t,
##     akt1_c5854t_c7996t, akt1_c6024t_c8166t, akt1_c6148t_c8290t,
##     akt1_c832g_c3359g, akt1_c9756a_c11898t, akt1_g14803t,
##     akt1_g15129a, akt1_g1780a_g363a, akt1_g20703a, akt1_g22187a,
##     akt1_g23477a, akt1_g2347t_g205t, akt1_g2375a_g233a,
##     akt1_g288c, akt1_g4362c, akt1_t10598a_t12740a,
##     akt1_t10726c_t12868c, akt1_t22932c, akt1_t8407g, akt2_2304186,
##     akt2_7254617, akt2_969531, akt2_rs892118, ampd1_c34t,
##     ankrd6_a545t, ankrd6_a550t, ankrd6_g197805a, ankrd6_i128l,
##     ankrd6_k710x, ankrd6_m485l, ankrd6_p636l, ankrd6_q122e,
##     ankrd6_t233m, apoe_c472t, ardb1_1801253, b2b, bcl6_1056932,
##     bcl6_17797517, bcl6_3774298, bcl6_4686467, bmp2_235768,
##     bmp2_rs15705, c8orf68_rs6983944, Calc.post.BMI,
##     carp.exon3_snp3, carp_3utr, carp_a8470g, carp_c105t,
##     carp_exon3_snp1, carp_exon3_snp2, carp_exon5_snp1,
##     carp_exon5_snp2, carp_exon6_ag, carp_exon8_ct, cast_rs754615,
##     cast_rs7724759, cav2_q130e, Center, CHOL, CHOL_HDL_C,
##     CK_day_0, CK_day_2, CK_day_4, cntf_g6a, cox5b_11904110,
##     cox5b_17022045, CRP, ctsf_572846, D_Post_B_VOL, D_Post_BM_VOL,
##     D_Post_F_VOL, D_Post_M_VOL, D_Post_WA_VOL, D_Post_WM_VOL,
##     D_Pre_B_VOL, D_Pre_BM_VOL, D_Pre_F_VOL, D_Pre_M_VOL,
##     D_Pre_WA_VOL, D_Pre_WM_VOL, D23.CH, D23_DIFF, DBP,
##     ddit_rs1053227, Dom.Arm, DRM.CH, DRM_DIFF, esr1_rs1042717,
##     esr1_rs1801132, esr1_rs2077647, esr1_rs2228480,
##     esr1_rs2234693, esr1_rs9340799, esr2_1256112, esr2_3020450,
##     fatso, fbox32_rs2891770, fbox32_rs3739287, fbox32_rs4871385,
##     fbox32_rs6690663, fbox32_rs6983944, FINS, FLGU, fst_1423560,
##     fst_722910, fstl1_rs9631455, galntl1_rs8005903, gapd_7307229,
##     gapd_7971637, gdf8_1225t, gdf8_e164k, gdf8_k153r, gdf8_p198a,
##     gdf8_rs1805085, gdf8_rs1805086, gdf8_t55a, Gender,
##     gnb3_rs5443, gs_s287nga, HDL_C, HOMA, hsd11b1_rs3753519,
##     hsd11b1_rs4844880, hsd11b1_rs846907, hsd11b1_rs846909,
##     hsd11b1_rs846910, hsd11b1_rs846911, hsd11b1_rs860185, id,
##     igf1_2288377, igf1_pro, igf1_t1245c, igf2_2230949,
##     igf2_3213216, igf2_3213220, igf2_3213221, igf2_rs680,
##     igfbp3_2132570, igfbp3_6670, il15_1057972, il15_1589241,
##     il15_rs2296135, il15ra_2228059, il15ra_3136618, il6_c174g,
##     il6_g572c, insig2_rs7566605, irs1_g972r, kchj11_rs5219, LDL_C,
##     lep_2167270, lepr_1137100, lepr_1137101, lepr_1805096,
##     lepr_8179183, lrp5_snp2, lrp5_snp3, lrp5_snp4, Mean_BP,
##     Met_syn, mgst3_4147542, mgst3_7549530, mlck_c37885a,
##     mlck_c49t, mlck_g91689t, mulk_3757470, mylk_c37885a,
##     mylk_g91689t, myod1_rs2249104, ND_Post_B_VOL, ND_Post_BM_VOL,
##     ND_Post_F_VOL, ND_Post_M_VOL, ND_Post_WA_VOL, ND_Post_WM_VOL,
##     ND_Pre_B_VOL, ND_Pre_BM_VOL, ND_Pre_F_VOL, ND_Pre_M_VOL,
##     ND_Pre_WA_VOL, ND_Pre_WM_VOL, ND23.CH, ND23_DIFF, NDRM.CH,
##     NDRM_DIFF, nnmt_g5082t, nos3_rs1799983, nr3c1_rs10482614,
##     nr3c1_rs10482616, nr3c1_rs4582314, nr3c1_rs4634384,
##     nr3c1_rs6195, opg_2073618, opg_3102735, p2rx2_6560891,
##     p2rx2_7966242, p2ry2_2511241, p2ry2_rs1783596, pai1_4g5g,
##     pcr13_snp1, pcr15_snp1, pcr15_snp2, pcr15_snp3, pcr15_snp4,
##     pcr5_snp1, pcr5_snp2, pcr8_snp1, pgc1_g76039a, pgm1_2284998,
##     pik3_rs3173908, pik3cg_4727666, Post.Height, Post.weight,
##     Post_D_avg, Post_DRM_Max, Post_LBi_Avg, Post_LTri_Avg,
##     Post_ND_avg, Post_NDRM_Max, Post_RBi_Avg, Post_RTri_Avg,
##     Post1_D_AVG, Post1_D1, Post1_D2, Post1_D3, Post1_ND_AVG,
##     Post1_ND1, Post1_ND2, Post1_ND3, Post2_D_AVG, Post2_D1,
##     Post2_D2, Post2_D3, Post2_ND_AVG, Post2_ND1, Post2_ND2,
##     Post2_ND3, ppar_gp12a, ppara_135539, ppara_1800206,
##     ppara_4252778, pparg_c1a_rs8192678, prdx6_4382766, pre.BMI,
##     Pre.height, Pre.HR, Pre.weight, Pre_DRM_Max, Pre_LBi_Avg,
##     Pre_LTri_Avg, Pre_NDRM_Max, Pre_RBi_Avg, Pre_RTri_Avg,
##     pygm_620006, Race, rankl_17458177, rankl_17536280,
##     rankl_17639305, rankl_4531631, resistin_a537c, resistin_c180g,
##     resistin_c30t, resistin_c398t, resistin_c980g, resistin_g540a,
##     rnf28_c16242t, rs11630261, rs302964, rs849409, runx_12526462,
##     SBP, slc35f1_rs10484290, slc35f1_rs12110370, Status,
##     syngr2_c886t, tcfl72_12255372, tcfl72_7903146,
##     tcfl72_rs12255372, tcfl72_rs7903146, tdp52l1_9321028, Term,
##     TG, tnfa_g308a, tnni1_2799672, tnni1_3738287, tpd52l1_3778458,
##     tpd52l1_3799736, tpd52l1_4896782, tpd52l1_514096, ucp2_ct,
##     V1_D_AVG, V1_D1, V1_D2, V1_D3, V1_ND_AVG, V1_ND1, V1_ND2,
##     V1_ND3, V123_D_AVG, V123_ND_AVG, V2_D_AVG, V2_D1, V2_D2,
##     V2_D3, V2_ND_AVG, V2_ND1, V2_ND2, V2_ND3, V23_D_AVG,
##     V23_ND_AVG, V3_D_AVG, V3_D1, V3_D2, V3_D3, V3_ND_AVG, V3_ND1,
##     V3_ND2, V3_ND3, vdr_bsm1, vdr_fok1, vdr_rs731236, vdr_taq1,
##     vim_6602186, visfatin_10487820, visfatin_10953502,
##     visfatin_2041681, visfatin_3801272, visfatin_6947766,
##     visfatin_929604, VLDL_TG
HapMat <- HapDesign(HaploEM)
Trait <- NDRM.CH[Race=="Caucasian" & !is.na(Race)]
mod1 <- (lm(Trait~HapMat))
mod2 <- (lm(Trait~1))
anova(mod2,mod1)
## Analysis of Variance Table
## 
## Model 1: Trait ~ 1
## Model 2: Trait ~ HapMat
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    776 881666                           
## 2    766 869272 10     12394 1.0921 0.3653
#EM for estimation and testing of haplotype{trait association
attach(fms)
## The following objects are masked from fms (pos = 3):
## 
##     acdc_rs1501299, ace_id, actn3_1671064, actn3_r577x,
##     actn3_rs1815739, actn3_rs540874, adrb2_1042713, adrb2_1042714,
##     adrb2_rs1042718, adrb3_4994, Age, agrp_5030980, akt1_a15756t,
##     akt1_a22889g, akt1_a7699g, akt1_c10744t_c12886t, akt1_c15676t,
##     akt1_c5854t_c7996t, akt1_c6024t_c8166t, akt1_c6148t_c8290t,
##     akt1_c832g_c3359g, akt1_c9756a_c11898t, akt1_g14803t,
##     akt1_g15129a, akt1_g1780a_g363a, akt1_g20703a, akt1_g22187a,
##     akt1_g23477a, akt1_g2347t_g205t, akt1_g2375a_g233a,
##     akt1_g288c, akt1_g4362c, akt1_t10598a_t12740a,
##     akt1_t10726c_t12868c, akt1_t22932c, akt1_t8407g, akt2_2304186,
##     akt2_7254617, akt2_969531, akt2_rs892118, ampd1_c34t,
##     ankrd6_a545t, ankrd6_a550t, ankrd6_g197805a, ankrd6_i128l,
##     ankrd6_k710x, ankrd6_m485l, ankrd6_p636l, ankrd6_q122e,
##     ankrd6_t233m, apoe_c472t, ardb1_1801253, b2b, bcl6_1056932,
##     bcl6_17797517, bcl6_3774298, bcl6_4686467, bmp2_235768,
##     bmp2_rs15705, c8orf68_rs6983944, Calc.post.BMI,
##     carp.exon3_snp3, carp_3utr, carp_a8470g, carp_c105t,
##     carp_exon3_snp1, carp_exon3_snp2, carp_exon5_snp1,
##     carp_exon5_snp2, carp_exon6_ag, carp_exon8_ct, cast_rs754615,
##     cast_rs7724759, cav2_q130e, Center, CHOL, CHOL_HDL_C,
##     CK_day_0, CK_day_2, CK_day_4, cntf_g6a, cox5b_11904110,
##     cox5b_17022045, CRP, ctsf_572846, D_Post_B_VOL, D_Post_BM_VOL,
##     D_Post_F_VOL, D_Post_M_VOL, D_Post_WA_VOL, D_Post_WM_VOL,
##     D_Pre_B_VOL, D_Pre_BM_VOL, D_Pre_F_VOL, D_Pre_M_VOL,
##     D_Pre_WA_VOL, D_Pre_WM_VOL, D23.CH, D23_DIFF, DBP,
##     ddit_rs1053227, Dom.Arm, DRM.CH, DRM_DIFF, esr1_rs1042717,
##     esr1_rs1801132, esr1_rs2077647, esr1_rs2228480,
##     esr1_rs2234693, esr1_rs9340799, esr2_1256112, esr2_3020450,
##     fatso, fbox32_rs2891770, fbox32_rs3739287, fbox32_rs4871385,
##     fbox32_rs6690663, fbox32_rs6983944, FINS, FLGU, fst_1423560,
##     fst_722910, fstl1_rs9631455, galntl1_rs8005903, gapd_7307229,
##     gapd_7971637, gdf8_1225t, gdf8_e164k, gdf8_k153r, gdf8_p198a,
##     gdf8_rs1805085, gdf8_rs1805086, gdf8_t55a, Gender,
##     gnb3_rs5443, gs_s287nga, HDL_C, HOMA, hsd11b1_rs3753519,
##     hsd11b1_rs4844880, hsd11b1_rs846907, hsd11b1_rs846909,
##     hsd11b1_rs846910, hsd11b1_rs846911, hsd11b1_rs860185, id,
##     igf1_2288377, igf1_pro, igf1_t1245c, igf2_2230949,
##     igf2_3213216, igf2_3213220, igf2_3213221, igf2_rs680,
##     igfbp3_2132570, igfbp3_6670, il15_1057972, il15_1589241,
##     il15_rs2296135, il15ra_2228059, il15ra_3136618, il6_c174g,
##     il6_g572c, insig2_rs7566605, irs1_g972r, kchj11_rs5219, LDL_C,
##     lep_2167270, lepr_1137100, lepr_1137101, lepr_1805096,
##     lepr_8179183, lrp5_snp2, lrp5_snp3, lrp5_snp4, Mean_BP,
##     Met_syn, mgst3_4147542, mgst3_7549530, mlck_c37885a,
##     mlck_c49t, mlck_g91689t, mulk_3757470, mylk_c37885a,
##     mylk_g91689t, myod1_rs2249104, ND_Post_B_VOL, ND_Post_BM_VOL,
##     ND_Post_F_VOL, ND_Post_M_VOL, ND_Post_WA_VOL, ND_Post_WM_VOL,
##     ND_Pre_B_VOL, ND_Pre_BM_VOL, ND_Pre_F_VOL, ND_Pre_M_VOL,
##     ND_Pre_WA_VOL, ND_Pre_WM_VOL, ND23.CH, ND23_DIFF, NDRM.CH,
##     NDRM_DIFF, nnmt_g5082t, nos3_rs1799983, nr3c1_rs10482614,
##     nr3c1_rs10482616, nr3c1_rs4582314, nr3c1_rs4634384,
##     nr3c1_rs6195, opg_2073618, opg_3102735, p2rx2_6560891,
##     p2rx2_7966242, p2ry2_2511241, p2ry2_rs1783596, pai1_4g5g,
##     pcr13_snp1, pcr15_snp1, pcr15_snp2, pcr15_snp3, pcr15_snp4,
##     pcr5_snp1, pcr5_snp2, pcr8_snp1, pgc1_g76039a, pgm1_2284998,
##     pik3_rs3173908, pik3cg_4727666, Post.Height, Post.weight,
##     Post_D_avg, Post_DRM_Max, Post_LBi_Avg, Post_LTri_Avg,
##     Post_ND_avg, Post_NDRM_Max, Post_RBi_Avg, Post_RTri_Avg,
##     Post1_D_AVG, Post1_D1, Post1_D2, Post1_D3, Post1_ND_AVG,
##     Post1_ND1, Post1_ND2, Post1_ND3, Post2_D_AVG, Post2_D1,
##     Post2_D2, Post2_D3, Post2_ND_AVG, Post2_ND1, Post2_ND2,
##     Post2_ND3, ppar_gp12a, ppara_135539, ppara_1800206,
##     ppara_4252778, pparg_c1a_rs8192678, prdx6_4382766, pre.BMI,
##     Pre.height, Pre.HR, Pre.weight, Pre_DRM_Max, Pre_LBi_Avg,
##     Pre_LTri_Avg, Pre_NDRM_Max, Pre_RBi_Avg, Pre_RTri_Avg,
##     pygm_620006, Race, rankl_17458177, rankl_17536280,
##     rankl_17639305, rankl_4531631, resistin_a537c, resistin_c180g,
##     resistin_c30t, resistin_c398t, resistin_c980g, resistin_g540a,
##     rnf28_c16242t, rs11630261, rs302964, rs849409, runx_12526462,
##     SBP, slc35f1_rs10484290, slc35f1_rs12110370, Status,
##     syngr2_c886t, tcfl72_12255372, tcfl72_7903146,
##     tcfl72_rs12255372, tcfl72_rs7903146, tdp52l1_9321028, Term,
##     TG, tnfa_g308a, tnni1_2799672, tnni1_3738287, tpd52l1_3778458,
##     tpd52l1_3799736, tpd52l1_4896782, tpd52l1_514096, ucp2_ct,
##     V1_D_AVG, V1_D1, V1_D2, V1_D3, V1_ND_AVG, V1_ND1, V1_ND2,
##     V1_ND3, V123_D_AVG, V123_ND_AVG, V2_D_AVG, V2_D1, V2_D2,
##     V2_D3, V2_ND_AVG, V2_ND1, V2_ND2, V2_ND3, V23_D_AVG,
##     V23_ND_AVG, V3_D_AVG, V3_D1, V3_D2, V3_D3, V3_ND_AVG, V3_ND1,
##     V3_ND2, V3_ND3, vdr_bsm1, vdr_fok1, vdr_rs731236, vdr_taq1,
##     vim_6602186, visfatin_10487820, visfatin_10953502,
##     visfatin_2041681, visfatin_3801272, visfatin_6947766,
##     visfatin_929604, VLDL_TG
## The following objects are masked from fms (pos = 4):
## 
##     acdc_rs1501299, ace_id, actn3_1671064, actn3_r577x,
##     actn3_rs1815739, actn3_rs540874, adrb2_1042713, adrb2_1042714,
##     adrb2_rs1042718, adrb3_4994, Age, agrp_5030980, akt1_a15756t,
##     akt1_a22889g, akt1_a7699g, akt1_c10744t_c12886t, akt1_c15676t,
##     akt1_c5854t_c7996t, akt1_c6024t_c8166t, akt1_c6148t_c8290t,
##     akt1_c832g_c3359g, akt1_c9756a_c11898t, akt1_g14803t,
##     akt1_g15129a, akt1_g1780a_g363a, akt1_g20703a, akt1_g22187a,
##     akt1_g23477a, akt1_g2347t_g205t, akt1_g2375a_g233a,
##     akt1_g288c, akt1_g4362c, akt1_t10598a_t12740a,
##     akt1_t10726c_t12868c, akt1_t22932c, akt1_t8407g, akt2_2304186,
##     akt2_7254617, akt2_969531, akt2_rs892118, ampd1_c34t,
##     ankrd6_a545t, ankrd6_a550t, ankrd6_g197805a, ankrd6_i128l,
##     ankrd6_k710x, ankrd6_m485l, ankrd6_p636l, ankrd6_q122e,
##     ankrd6_t233m, apoe_c472t, ardb1_1801253, b2b, bcl6_1056932,
##     bcl6_17797517, bcl6_3774298, bcl6_4686467, bmp2_235768,
##     bmp2_rs15705, c8orf68_rs6983944, Calc.post.BMI,
##     carp.exon3_snp3, carp_3utr, carp_a8470g, carp_c105t,
##     carp_exon3_snp1, carp_exon3_snp2, carp_exon5_snp1,
##     carp_exon5_snp2, carp_exon6_ag, carp_exon8_ct, cast_rs754615,
##     cast_rs7724759, cav2_q130e, Center, CHOL, CHOL_HDL_C,
##     CK_day_0, CK_day_2, CK_day_4, cntf_g6a, cox5b_11904110,
##     cox5b_17022045, CRP, ctsf_572846, D_Post_B_VOL, D_Post_BM_VOL,
##     D_Post_F_VOL, D_Post_M_VOL, D_Post_WA_VOL, D_Post_WM_VOL,
##     D_Pre_B_VOL, D_Pre_BM_VOL, D_Pre_F_VOL, D_Pre_M_VOL,
##     D_Pre_WA_VOL, D_Pre_WM_VOL, D23.CH, D23_DIFF, DBP,
##     ddit_rs1053227, Dom.Arm, DRM.CH, DRM_DIFF, esr1_rs1042717,
##     esr1_rs1801132, esr1_rs2077647, esr1_rs2228480,
##     esr1_rs2234693, esr1_rs9340799, esr2_1256112, esr2_3020450,
##     fatso, fbox32_rs2891770, fbox32_rs3739287, fbox32_rs4871385,
##     fbox32_rs6690663, fbox32_rs6983944, FINS, FLGU, fst_1423560,
##     fst_722910, fstl1_rs9631455, galntl1_rs8005903, gapd_7307229,
##     gapd_7971637, gdf8_1225t, gdf8_e164k, gdf8_k153r, gdf8_p198a,
##     gdf8_rs1805085, gdf8_rs1805086, gdf8_t55a, Gender,
##     gnb3_rs5443, gs_s287nga, HDL_C, HOMA, hsd11b1_rs3753519,
##     hsd11b1_rs4844880, hsd11b1_rs846907, hsd11b1_rs846909,
##     hsd11b1_rs846910, hsd11b1_rs846911, hsd11b1_rs860185, id,
##     igf1_2288377, igf1_pro, igf1_t1245c, igf2_2230949,
##     igf2_3213216, igf2_3213220, igf2_3213221, igf2_rs680,
##     igfbp3_2132570, igfbp3_6670, il15_1057972, il15_1589241,
##     il15_rs2296135, il15ra_2228059, il15ra_3136618, il6_c174g,
##     il6_g572c, insig2_rs7566605, irs1_g972r, kchj11_rs5219, LDL_C,
##     lep_2167270, lepr_1137100, lepr_1137101, lepr_1805096,
##     lepr_8179183, lrp5_snp2, lrp5_snp3, lrp5_snp4, Mean_BP,
##     Met_syn, mgst3_4147542, mgst3_7549530, mlck_c37885a,
##     mlck_c49t, mlck_g91689t, mulk_3757470, mylk_c37885a,
##     mylk_g91689t, myod1_rs2249104, ND_Post_B_VOL, ND_Post_BM_VOL,
##     ND_Post_F_VOL, ND_Post_M_VOL, ND_Post_WA_VOL, ND_Post_WM_VOL,
##     ND_Pre_B_VOL, ND_Pre_BM_VOL, ND_Pre_F_VOL, ND_Pre_M_VOL,
##     ND_Pre_WA_VOL, ND_Pre_WM_VOL, ND23.CH, ND23_DIFF, NDRM.CH,
##     NDRM_DIFF, nnmt_g5082t, nos3_rs1799983, nr3c1_rs10482614,
##     nr3c1_rs10482616, nr3c1_rs4582314, nr3c1_rs4634384,
##     nr3c1_rs6195, opg_2073618, opg_3102735, p2rx2_6560891,
##     p2rx2_7966242, p2ry2_2511241, p2ry2_rs1783596, pai1_4g5g,
##     pcr13_snp1, pcr15_snp1, pcr15_snp2, pcr15_snp3, pcr15_snp4,
##     pcr5_snp1, pcr5_snp2, pcr8_snp1, pgc1_g76039a, pgm1_2284998,
##     pik3_rs3173908, pik3cg_4727666, Post.Height, Post.weight,
##     Post_D_avg, Post_DRM_Max, Post_LBi_Avg, Post_LTri_Avg,
##     Post_ND_avg, Post_NDRM_Max, Post_RBi_Avg, Post_RTri_Avg,
##     Post1_D_AVG, Post1_D1, Post1_D2, Post1_D3, Post1_ND_AVG,
##     Post1_ND1, Post1_ND2, Post1_ND3, Post2_D_AVG, Post2_D1,
##     Post2_D2, Post2_D3, Post2_ND_AVG, Post2_ND1, Post2_ND2,
##     Post2_ND3, ppar_gp12a, ppara_135539, ppara_1800206,
##     ppara_4252778, pparg_c1a_rs8192678, prdx6_4382766, pre.BMI,
##     Pre.height, Pre.HR, Pre.weight, Pre_DRM_Max, Pre_LBi_Avg,
##     Pre_LTri_Avg, Pre_NDRM_Max, Pre_RBi_Avg, Pre_RTri_Avg,
##     pygm_620006, Race, rankl_17458177, rankl_17536280,
##     rankl_17639305, rankl_4531631, resistin_a537c, resistin_c180g,
##     resistin_c30t, resistin_c398t, resistin_c980g, resistin_g540a,
##     rnf28_c16242t, rs11630261, rs302964, rs849409, runx_12526462,
##     SBP, slc35f1_rs10484290, slc35f1_rs12110370, Status,
##     syngr2_c886t, tcfl72_12255372, tcfl72_7903146,
##     tcfl72_rs12255372, tcfl72_rs7903146, tdp52l1_9321028, Term,
##     TG, tnfa_g308a, tnni1_2799672, tnni1_3738287, tpd52l1_3778458,
##     tpd52l1_3799736, tpd52l1_4896782, tpd52l1_514096, ucp2_ct,
##     V1_D_AVG, V1_D1, V1_D2, V1_D3, V1_ND_AVG, V1_ND1, V1_ND2,
##     V1_ND3, V123_D_AVG, V123_ND_AVG, V2_D_AVG, V2_D1, V2_D2,
##     V2_D3, V2_ND_AVG, V2_ND1, V2_ND2, V2_ND3, V23_D_AVG,
##     V23_ND_AVG, V3_D_AVG, V3_D1, V3_D2, V3_D3, V3_ND_AVG, V3_ND1,
##     V3_ND2, V3_ND3, vdr_bsm1, vdr_fok1, vdr_rs731236, vdr_taq1,
##     vim_6602186, visfatin_10487820, visfatin_10953502,
##     visfatin_2041681, visfatin_3801272, visfatin_6947766,
##     visfatin_929604, VLDL_TG
Geno.C <- setupGeno(Geno.C)
Trait <- NDRM.CH[Race=="Caucasian" & !is.na(Race)]
Dat <- data.frame(Geno.C=Geno.C, Trait=Trait)
library(haplo.stats)
haplo.glm(Trait~Geno.C,data=Dat, 
          allele.lev=attributes(Geno.C)$unique.alleles)
## 
## Call:  haplo.glm(formula = Trait ~ Geno.C, data = Dat, allele.lev = attributes(Geno.C)$unique.alleles)
## 
## Coefficients:
## (Intercept)     Geno.C.3     Geno.C.5     Geno.C.8     Geno.C.9  
##    50.67787      8.49595     -0.44085      2.01114      8.42214  
## Geno.C.rare  
##     3.98509  
## 
## Haplotypes:
##             loc.1 loc.2 loc.3 loc.4 hap.freq
## Geno.C.3        C     A     T     G 0.012490
## Geno.C.5        C     G     C     G 0.010776
## Geno.C.8        T     A     T     G 0.402424
## Geno.C.9        T     G     C     A 0.083942
## Geno.C.rare     *     *     *     * 0.018135
## haplo.base      C     G     C     A 0.472233
## 
## Degrees of Freedom:  761 Total (i.e. Null);  756 Residual
## 
## Subjects removed by NAs in y or x, or all NA in geno
##   yxmiss genomiss 
##       14       15 
## 
##      Null Deviance:  864660 
##  Residual Deviance:  853550 
##                AIC:  7526.6
haplo.glm(Trait~Geno.C,data=Dat,
          allele.lev=attributes(Geno.C)$unique.alleles, control=haplo.glm.control(haplo.effect="dominant"))
## 
## Call:  haplo.glm(formula = Trait ~ Geno.C, data = Dat, control = haplo.glm.control(haplo.effect = "dominant"), 
##     allele.lev = attributes(Geno.C)$unique.alleles)
## 
## Coefficients:
## (Intercept)     Geno.C.3     Geno.C.5     Geno.C.8     Geno.C.9  
##     48.9376       4.3548       1.3853       4.5781      13.9041  
## Geno.C.rare  
##      3.4349  
## 
## Haplotypes:
##             loc.1 loc.2 loc.3 loc.4 hap.freq
## Geno.C.3        C     A     T     G 0.012541
## Geno.C.5        C     G     C     G 0.010773
## Geno.C.8        T     A     T     G 0.402398
## Geno.C.9        T     G     C     A 0.083754
## Geno.C.rare     *     *     *     * 0.018130
## haplo.base      C     G     C     A 0.472405
## 
## Degrees of Freedom:  761 Total (i.e. Null);  756 Residual
## 
## Subjects removed by NAs in y or x, or all NA in geno
##   yxmiss genomiss 
##       14       15 
## 
##      Null Deviance:  864660 
##  Residual Deviance:  846380 
##                AIC:  7520.2