# R/agjack.pi0.R In pi0: Estimating the proportion of true null hypotheses for FDR

#### Defines functions `agjack.pi0`

`agjack.pi0` <-
function(subtobj,mean.n=c('mean','harmean','geomean') ,pointpair=FALSE, trunc=TRUE,tol=1e-5){
y=subtobj[,'f1']
n1=subtobj[,'n1']
n2=subtobj[,'n2']
min.n=pmin(n1,n2)
max.n=pmax(n1,n2)
mean.n=match.arg(mean.n)
mean.n=get(mean.n)
n=apply(cbind(n1,n2),1,mean.n)

if(!pointpair){
y.mean=aggregate(y,list(min.n,max.n),mean)
lev=apply(y.mean[,1:2],2,function(xx)as.numeric(xx))
y.mean=y.mean[,3]
n.unique=apply(lev,1,mean.n)
nns=length(n.unique)
rslt=0
k=0
for(i in 1:(nns-1))
for(j in (i+1):nns)
if(abs(sqrt((n.unique[j]+2)/(n.unique[i]+2))-1)>tol) {
k=k+1
rslt=rslt+if(trunc){
min(1,max(0,gjack(y.mean[i],y.mean[j],sqrt((n.unique[j]+2)/(n.unique[i]+2))),na.rm=TRUE),na.rm=TRUE)
}else{gjack(y.mean[i],y.mean[j],sqrt((n.unique[j]+2)/(n.unique[i]+2)))}
}
return(rslt/k)
}else{
nns=length(n)
rslt=0
k=0
for(i in 1:nns)
for(j in 1:nns)
if(abs(sqrt((n[j]+2)/(n[i]+2))-1)>tol) {
k=k+1
rslt=rslt+if(trunc){
min(1,max(0,gjack(y[i],y[j],sqrt((n[j]+2)/(n[i]+2))),na.rm=TRUE),na.rm=TRUE)
}else{gjack(y[i],y[j],sqrt((n[j]+2)/(n[i]+2)))}
}
return(rslt/k)
}
}

## Try the pi0 package in your browser

Any scripts or data that you put into this service are public.

pi0 documentation built on May 2, 2019, 4:47 p.m.