#Bonferroni adjustment
virco <- read.csv(file="Virco_data.csv", header=T, sep=",")
attach(virco)
PrMut <- virco[,23:121]!="-" & virco[,23:121]!="."
NObs <- dim(virco)[1]
PrMutSub <-data.frame(PrMut[ , apply(PrMut,2,sum) > NObs*.05])
Trait <- IDV.Fold - NFV.Fold
TtestP <- function(Geno){return(t.test(Trait[Geno==1],
Trait[Geno==0], na.rm=T)$"p.value")
}
Pvec <- apply(PrMutSub, 2, TtestP)
sort(Pvec)
## P30 P76 P88 P55 P48
## 3.732500e-12 9.782323e-10 1.432468e-06 2.286695e-06 5.749467e-06
## P89 P11 P82 P60 P85
## 8.924013e-05 4.171618e-04 9.500604e-04 1.115441e-03 1.219064e-03
## P54 P43 P61 P46 P67
## 1.489381e-03 2.025621e-03 2.556156e-03 4.198935e-03 7.765537e-03
## P69 P84 P47 P35 P32
## 1.113762e-02 1.557464e-02 1.574864e-02 2.392427e-02 2.508445e-02
## P33 P14 P16 P72 P13
## 2.722251e-02 3.441981e-02 5.570492e-02 5.748494e-02 6.375590e-02
## P15 P34 P53 P64 P90
## 1.089171e-01 1.167541e-01 1.556130e-01 2.540249e-01 2.618606e-01
## P36 P63 P37 P77 P24
## 2.896151e-01 2.945370e-01 3.257741e-01 3.356589e-01 3.441678e-01
## P93 P71 P10 P74 P73
## 3.619516e-01 3.761893e-01 4.268153e-01 4.480744e-01 4.906612e-01
## P20 P19 P58 P62 P41
## 5.311825e-01 5.342250e-01 5.440101e-01 6.677043e-01 6.998280e-01
## P12 P57
## 8.050362e-01 9.938846e-01
names(PrMutSub)[Pvec < 0.05]
## [1] "P11" "P14" "P30" "P32" "P33" "P35" "P43" "P46" "P47" "P48" "P54"
## [12] "P55" "P60" "P61" "P67" "P69" "P76" "P82" "P84" "P85" "P88" "P89"
PvecAdj <- p.adjust(Pvec, method="bonferroni")
sort(PvecAdj)
## P30 P76 P88 P55 P48
## 1.754275e-10 4.597692e-08 6.732600e-05 1.074747e-04 2.702250e-04
## P89 P11 P82 P60 P85
## 4.194286e-03 1.960660e-02 4.465284e-02 5.242573e-02 5.729603e-02
## P54 P43 P61 P46 P67
## 7.000090e-02 9.520419e-02 1.201393e-01 1.973500e-01 3.649803e-01
## P69 P84 P47 P10 P12
## 5.234681e-01 7.320083e-01 7.401862e-01 1.000000e+00 1.000000e+00
## P13 P14 P15 P16 P19
## 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## P20 P24 P32 P33 P34
## 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## P35 P36 P37 P41 P53
## 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## P57 P58 P62 P63 P64
## 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## P71 P72 P73 P74 P77
## 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## P90 P93
## 1.000000e+00 1.000000e+00
names(PrMutSub)[PvecAdj < 0.05]
## [1] "P11" "P30" "P48" "P55" "P76" "P82" "P88" "P89"
#Tukey's single-step method
fmsURL <- "http://www.stat-gen.org/book.e1/data/FMS_data.txt"
fms <- read.delim(file=fmsURL, header=T, sep="\t")
attach(fms)
Trait <- NDRM.CH
summary(lm(Trait~resistin_c180g))
##
## Call:
## lm(formula = Trait ~ resistin_c180g)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.054 -22.754 -6.054 15.346 193.946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.054 2.004 27.973 <2e-16 ***
## resistin_c180gCG -5.918 2.864 -2.067 0.0392 *
## resistin_c180gGG -4.553 4.356 -1.045 0.2964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.05 on 603 degrees of freedom
## (791 observations deleted due to missingness)
## Multiple R-squared: 0.007296, Adjusted R-squared: 0.004003
## F-statistic: 2.216 on 2 and 603 DF, p-value: 0.11
TukeyHSD(aov(Trait~resistin_c180g))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Trait ~ resistin_c180g)
##
## $resistin_c180g
## diff lwr upr p adj
## CG-CC -5.917630 -12.645660 0.8103998 0.0977410
## GG-CC -4.553042 -14.788156 5.6820721 0.5486531
## GG-CG 1.364588 -8.916062 11.6452381 0.9478070
#Benjamini and Hochberg (B-H) adjustment
Pvec <- as.vector(Pvec)
m <- length(Pvec)
BHp <- sort(Pvec,decreasing=T)*m/seq(m,1)
sort(cummin(BHp))
## [1] 1.754275e-10 2.298846e-08 2.244200e-05 2.686866e-05 5.404499e-05
## [6] 6.990477e-04 2.800943e-03 5.581605e-03 5.729603e-03 5.729603e-03
## [11] 6.363718e-03 7.933683e-03 9.241488e-03 1.409643e-02 2.433202e-02
## [16] 3.271676e-02 4.112145e-02 4.112145e-02 5.894845e-02 5.894845e-02
## [21] 6.092657e-02 7.353323e-02 1.125747e-01 1.125747e-01 1.198611e-01
## [26] 1.968885e-01 2.032385e-01 2.612075e-01 4.102483e-01 4.102483e-01
## [31] 4.326013e-01 4.326013e-01 4.621682e-01 4.621682e-01 4.621682e-01
## [36] 4.725479e-01 4.778621e-01 5.279031e-01 5.399871e-01 5.765269e-01
## [41] 5.946157e-01 5.946157e-01 5.946157e-01 7.132296e-01 7.309314e-01
## [46] 8.225370e-01 9.938846e-01
BHp[order(Pvec,decreasing=T)] <- cummin(BHp)
names(PrMutSub)[BHp < 0.05]
## [1] "P11" "P30" "P43" "P46" "P47" "P48" "P54" "P55" "P60" "P61" "P67"
## [12] "P69" "P76" "P82" "P84" "P85" "P88" "P89"
sort(p.adjust(Pvec, method="BH"))
## [1] 1.754275e-10 2.298846e-08 2.244200e-05 2.686866e-05 5.404499e-05
## [6] 6.990477e-04 2.800943e-03 5.581605e-03 5.729603e-03 5.729603e-03
## [11] 6.363718e-03 7.933683e-03 9.241488e-03 1.409643e-02 2.433202e-02
## [16] 3.271676e-02 4.112145e-02 4.112145e-02 5.894845e-02 5.894845e-02
## [21] 6.092657e-02 7.353323e-02 1.125747e-01 1.125747e-01 1.198611e-01
## [26] 1.968885e-01 2.032385e-01 2.612075e-01 4.102483e-01 4.102483e-01
## [31] 4.326013e-01 4.326013e-01 4.621682e-01 4.621682e-01 4.621682e-01
## [36] 4.725479e-01 4.778621e-01 5.279031e-01 5.399871e-01 5.765269e-01
## [41] 5.946157e-01 5.946157e-01 5.946157e-01 7.132296e-01 7.309314e-01
## [46] 8.225370e-01 9.938846e-01
#Benjamini and Yekutieli adjustment
BYp <- p.adjust(Pvec, method="BY")
sort(BYp)
## [1] 7.785410e-10 1.020220e-07 9.959678e-05 1.192422e-04 2.398497e-04
## [6] 3.102348e-03 1.243048e-02 2.477096e-02 2.542777e-02 2.542777e-02
## [11] 2.824195e-02 3.520940e-02 4.101339e-02 6.255943e-02 1.079846e-01
## [16] 1.451958e-01 1.824955e-01 1.824955e-01 2.616111e-01 2.616111e-01
## [21] 2.703899e-01 3.263378e-01 4.996023e-01 4.996023e-01 5.319392e-01
## [26] 8.737842e-01 9.019653e-01 1.000000e+00 1.000000e+00 1.000000e+00
## [31] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [36] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [41] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [46] 1.000000e+00 1.000000e+00
names(PrMutSub)[BYp < 0.05]
## [1] "P11" "P30" "P43" "P48" "P54" "P55" "P60" "P61" "P76" "P82" "P85"
## [12] "P88" "P89"
#Calculation of the q-value
library(qvalue)
sort(qvalue(Pvec,lambda=0)$qvalues)
## [1] 1.754275e-10 2.298846e-08 2.244200e-05 2.686866e-05 5.404499e-05
## [6] 6.990477e-04 2.800943e-03 5.581605e-03 5.729603e-03 5.729603e-03
## [11] 6.363718e-03 7.933683e-03 9.241488e-03 1.409643e-02 2.433202e-02
## [16] 3.271676e-02 4.112145e-02 4.112145e-02 5.894845e-02 5.894845e-02
## [21] 6.092657e-02 7.353323e-02 1.125747e-01 1.125747e-01 1.198611e-01
## [26] 1.968885e-01 2.032385e-01 2.612075e-01 4.102483e-01 4.102483e-01
## [31] 4.326013e-01 4.326013e-01 4.621682e-01 4.621682e-01 4.621682e-01
## [36] 4.725479e-01 4.778621e-01 5.279031e-01 5.399871e-01 5.765269e-01
## [41] 5.946157e-01 5.946157e-01 5.946157e-01 7.132296e-01 7.309314e-01
## [46] 8.225370e-01 9.938846e-01
sort(qvalue(Pvec,pi0.method="bootstrap")$qvalues)
## [1] 3.317778e-11 4.347699e-09 4.244350e-06 5.081544e-06 1.022127e-05
## [6] 1.322076e-04 5.297292e-04 1.055623e-03 1.083613e-03 1.083613e-03
## [11] 1.203540e-03 1.500460e-03 1.747799e-03 2.665991e-03 4.601800e-03
## [16] 6.187566e-03 7.777107e-03 7.777107e-03 1.114864e-02 1.114864e-02
## [21] 1.152276e-02 1.390699e-02 2.129072e-02 2.129072e-02 2.266877e-02
## [26] 3.723660e-02 3.843755e-02 4.940095e-02 7.758834e-02 7.758834e-02
## [31] 8.181585e-02 8.181585e-02 8.740769e-02 8.740769e-02 8.740769e-02
## [36] 8.937077e-02 9.037581e-02 9.983983e-02 1.021252e-01 1.090358e-01
## [41] 1.124569e-01 1.124569e-01 1.124569e-01 1.348898e-01 1.382376e-01
## [46] 1.555626e-01 1.879687e-01
qvalue(Pvec,pi0.method="bootstrap")$pi0
## [1] 0.1891253
#Free step-down resampling adjustment
attach(fms)
## 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
Actn3Bin <- data.frame(actn3_r577x!="TT",actn3_rs540874!="AA",
actn3_rs1815739!="TT",actn3_1671064!="GG")
Mod <- summary(lm(NDRM.CH~.,data=Actn3Bin))
Mod
##
## Call:
## lm(formula = NDRM.CH ~ ., data = Actn3Bin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.181 -22.614 -7.414 15.486 198.786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 54.700 3.212 17.028 <2e-16 ***
## actn3_r577x.....TT.TRUE -12.891 4.596 -2.805 0.0052 **
## actn3_rs540874.....AA.TRUE 10.899 11.804 0.923 0.3562
## actn3_rs1815739.....TT.TRUE 27.673 17.876 1.548 0.1222
## actn3_1671064.....GG.TRUE -29.166 17.516 -1.665 0.0964 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.93 on 591 degrees of freedom
## (801 observations deleted due to missingness)
## Multiple R-squared: 0.01945, Adjusted R-squared: 0.01281
## F-statistic: 2.93 on 4 and 591 DF, p-value: 0.02037
TestStatObs <- Mod$coefficients[-1,3]
Tobs <- as.vector(sort(abs(TestStatObs)))
MissDat <- apply(is.na(Actn3Bin),1,any) | is.na(NDRM.CH)
Actn3BinC <- Actn3Bin[!MissDat,]
Ord <- order(abs(TestStatObs))
M <- 1000
NSnps <- 4
Nobs <- sum(!MissDat)
TestStatResamp <- matrix(nrow=M, ncol=NSnps)
for (i in 1:M){
Ynew <- sample(Mod$residuals, size=Nobs, replace=T)
ModResamp <- summary(lm(Ynew~., data=Actn3BinC))
TestStatResamp[i,] <- abs(ModResamp$coefficients[-1,3])[Ord]
}
Qmat <- t(apply(TestStatResamp, 1, cummax))
Padj <- apply(t(matrix(rep(Tobs,M), NSnps)) < Qmat, 2, mean)
Padj
## [1] 0.349 0.224 0.221 0.036
#Null unrestricted bootstrap approach
CoefObs <- as.vector(Mod$coefficients[-1,1])
B <-1000
TestStatBoot <- matrix(nrow=B,ncol=NSnps)
for (i in 1:B){
SampID <- sample(1:Nobs,size=Nobs, replace=T)
Ynew <- NDRM.CH[!MissDat][SampID]
Xnew <- Actn3BinC[SampID,]
CoefBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,1]
SEBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,2]
if (length(CoefBoot)==length(CoefObs)){
TestStatBoot[i,] <- (CoefBoot-CoefObs)/SEBoot
}
}
for (cj in seq(2.7,2.8,.01)){
print(cj)
print(mean(apply(abs(TestStatBoot)>cj,1,sum)>=1,na.rm=T))
}
## [1] 2.7
## [1] 0.06079027
## [1] 2.71
## [1] 0.0597771
## [1] 2.72
## [1] 0.0597771
## [1] 2.73
## [1] 0.0597771
## [1] 2.74
## [1] 0.0597771
## [1] 2.75
## [1] 0.0597771
## [1] 2.76
## [1] 0.05876393
## [1] 2.77
## [1] 0.05876393
## [1] 2.78
## [1] 0.05673759
## [1] 2.79
## [1] 0.0526849
## [1] 2.8
## [1] 0.05167173