Using R to realize Local FDR control procedure (Sun&Cai 2007)

秋窗风雨夕1q84

已于 2024-09-14 15:47:52 修改

阅读量131

点赞数 3

文章标签： r语言开发语言

于 2024-09-14 15:12:10 首次发布

本文链接：https://blog.csdn.net/syqtheone/article/details/142259140

版权

############ estimate from TONY CAI's review
############ Lfdr procedure ##########
## 0.8*N(0,1) + 0.15*N(-2.5,1)+ 0.05*N(3,1)
set.seed(123)
n=5000
alpha=0.05
rnorm(80)
rnorm(15,-2.5,1)
rnorm(5,3,1)
X=c(rnorm(.8*n),rnorm(.15*n,-2.5,1),rnorm(.05*n,3,1))
names(X)=1:n
hist(X,breaks=50)

res<-locfdr::locfdr(X,nulltype = 0,plot=1) #plot=0
mean.emp = res$fp0["mlest","delta"]
sd.emp = res$fp0["mlest","sigma"]
fdr=res$fdr
names(fdr)=1:length(fdr)
Q=cumsum(sort(fdr))/1:length(fdr)
names(Q)<-1:length(fdr)
Q_set=Q[Q<=0.1]
j=max(as.numeric(names(Q_set)))
j
sprintf("#rejections: %d",j)

R_p_set=sort(fdr)[1:j]
R=names(sort(fdr)[1:j])
N11=setdiff(R,as.character(1:(.8*n)))
sprintf("#true rejections: %d",length(N11))
sprintf("FDP: %f",(j-length(N11))/j)
upper_cutoff=min(X[R][X[R]>0])
lower_cutoff=max(X[R][X[R]<0])
sprintf("lower cutoff: %f; upper cutoff: %f",lower_cutoff,upper_cutoff)

在这里插入图片描述

在这里插入图片描述
p0估计得不好

Efron的empirical distribution estimation control/correction 以及包locfdr 和 sun&cai的local fdr control在一起有些混淆了。
发现可以利用efron的包locfdr来做sun&cai 2007的local fdr control procedure, 只不过须要再附上一些代码步骤。
还没太理解efron包locfdr的用途，似乎是估计出empirical null? 但是p0的确估计得不好。（真实的p0应该是0.8）
Cai et al. (2007) 似乎也有估计pi0的方法，利用来特征函数？须要再看论文，理解比较。目前efron包locfdr里估计p0的方法应该是MLE 或 central matching estimates都可以估计。
local fdr procedure，不再利用p value，直接用f(x)(统计量密度) 也是一个点。
tony cai的review写的真的很好，可以很清晰地看到，fdr control中思想的阶段性跃迁变化。

参考文献

Large-Scale Global and Simultaneous Inference: Estimation and Testing in Very High Dimensions
T. Tony Cai1 and Wenguang Sun2. 2017. p17
R包locfdr的说明chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://cran.r-project.org/web/packages/locfdr/vignettes/locfdr-example.pdf
Cai, T., Jin, J., and Low, M. (2007), “Estimation and Confidence Sets for Sparse Normal Mixtures,” The Annals of Statistics, in press.