我正在为完全相同的问题寻找可行的解决方案。我发现的最好的是Foulkes Andrea 在他的书Applied Statistical Genetics with R(2009)中介绍的Null Unrestricted Bootstrap。与他专门考虑回归的所有其他文章和书籍相反。除了其他方法外,他还建议使用 Null Unrestricted Bootstrap,该方法适用于无法轻松计算残差的情况(如在我的情况下,我对许多独立回归(基本上是简单的相关性)进行建模,每个回归都具有相同的响应变量和不同的片段)。我发现这个方法也被称为maxT方法。
> attach(fms)
> Actn3Bin <- > data.frame(actn3_r577x!="TT",actn3_rs540874!="AA",actn3_rs1815739!="TT",actn3_1671064!="GG")
> Mod <- summary(lm(NDRM.CH~.,data=Actn3Bin))
> 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
+ }
+ }
一旦我们有了所有的TestStatBoot
矩阵(在行中我们有引导复制,在列中我们有引导T⃗ ∗^统计数据)我们找到了Tcrit.我们准确地观察α=0.05更显着的百分比T⃗ ∗^统计数据(更显着意味着绝对值大于Tcrit.)。
我们报告i-th 模型组件显着,如果它T⃗ i^>Tcrit.
最后一步可以用这段代码完成
p.value<-0.05 # The target alpha threshold
digits<-1000000
library(gtools) # for binsearch
pValueFun<-function(cj)
{
mean(apply(abs(TestStatBoot)>cj/digits,1,sum)>=1,na.rm=T)
}
ans<-binsearch(pValueFun,c(0.5*digits,100*digits),target=p.value)
p.level<-(1-pnorm(q=ans$where[[1]]/digits))*2 #two-sided.