BH和IHW之间的一个区别是，BH对所有假设设置了一个不变的拒绝阈值，而不依赖于协变量。另一方面，IHW的拒绝阈值取决于协变量。在许多情况下，这意味着BH排斥阈值将是非恒定的tdr函数，而IHW将近似于恒定的tdr阈值。在这里，我们想要使用一个简单的模拟来探索这个问题。

库("ggplot2")库("dplyr")库("tidyr")库("gridExtra") #库("LSD")库("IHW")

我们建立了辅助函数来生成模拟，应用BH和IHW并在拒绝阈值处评估tdr。

# try具有非单调关系的函数sim <- function(m =20000){x< - 1-runif(m) pi0 <- 0.75+0.3*(x-0.5)^2#0.8 #ifelse(x<=。5 x 0.5 + 0.6 * 0.8) eff_size < - 2 - 2.2 * abs (x - 0.5) ^ 1.3 # 3.2 # ifelse (x < =。5, 1.5, 1.5-(x-1/2)*4) H <- rbinom(m, 1,1 -pi0) Z <- rnorm(m) Z[H==1] <- rnorm(sum(H)， eff_size[H==1]) pvalue <- 1-pnorm(Z) t <- pvalue t_bh <- get_bh_threshold(pvalue,0.1) alt_ft_bh <- dnorm(-eff_size+ qnorm(1-t_bh))/dnorm(qnorm(1-t_bh)) tdr_bh <- alt_ft_bh*(1-pi0)/((1-pi0)*alt_ft_bh + pi0) t_ddhw <- thresholds(ihw(pvalue, x, alpha= .1, nbins=20, nsplits_internal=10, nfolds =1)，levels_only=FALSE) alt_ft_ddhw <- dnorm(-eff_size+ qnorm(1-t_ddhw))/dnorm(qnorm(1-t_ddhw)) tdr_ddhw <- alt_ft_ddhw*(1-pi0)/((1-pi0)*alt_ft_ddhw + pi0) alt_ft <- dnorm(qnorm(1-pi0) -eff_size)/dnorm(qnorm(1-pi0)*alt_ft + pi0) tdr <- alt_ft*(1-pi0)/((1-pi0)*alt_ft + pi0) return(data.frame(x=x,pi0=pi0, t=t, eff_size=eff_size,H=H, Z=Z, pvalue=pvalue, alt_ft=alt_ft, tdr=tdr, tdr_bh= alt_ft_bh, t_bh=t_bh, t_ddhw= tdr_ddhw, tdr_ddhw=tdr_ddhw, alt_ft_ddhw)}

运行以上模拟:

Set.seed (1) sim_df <- sim(m=80000)

##只使用1次折叠!只有当你想要学习权重时才使用这个方法，但千万不要用于测试!

首先根据p值绘制拒绝阈值。请注意，BH返回一条水平线，而IHW没有，因为分配的权重不均匀。

sim_df_t <- select(sim_df, x, t_bh, t_ddhw) %>% gather(method, t， -x) %>% mutate(method=ifelse(method=="t_bh"， "BH"，"IHW")) ggplot(sim_df_t, aes(x=x, y=t,color=method)) + geom_step() + xlab("covariate") + ylab("p-value threshold"))

另一方面，BH在tdr方面显示出高度非恒定的排斥阈值。这是次优的。IHW更好地逼近恒定的tdr阈值。

sim_df_tdr <- select(sim_df, x, tdr_bh, tdr_ddhw) %>% gather(method, tdr， -x) %>% mutate(method=ifelse(method=="tdr_bh"， "BH"，"IHW")) ggplot(sim_df_tdr, aes(x=x, y=tdr,color=method)) + geom_step() + xlab("covariate") + ylab("tdr threshold"))

#sim_df$group <- groups_by_filter(sim_df$x, 10) #plot(sim_df$x, sim_df$tdr_bh) #with(filter(sim_df，-log10(pvalue) > 0.001)， heatscatter(-log10(pvalue)，x, add.contour=TRUE)) #with(filter(sim_df,tdr>0.2)， heatscatter(tdr,x, add.contour=TRUE)) #ggplot(filter(sim_df)， aes(x=-log10(pvalue)， y=x)) + geom_density2d(aes(color = ..level..)) #ggplot(filter(sim_df, tdr>0.4)， aes(x= ..level..)) + geom_density2d(aes(x= ..level..))