# R/fuzzyBHexact.R In fuzzyFDR: Exact calculation of fuzzy decision rules for multiple testing

```"fuzzyBHexact" <-
function(pvals,pprev,alpha=0.05,tol=10e-6,q.myuni=T,dp=20){

pvals <- round(pvals,d=dp)
pprev <- round(pprev,d=dp)

## get the intervals D_j
if(q.myuni) p.all <- sort(myUnique(c(pvals,pprev),tol=tol))
else p.all <- sort(unique(c(pvals,pprev)))

#print(p.all)

n.int <- length(p.all)-1
n.tests <- length(pvals)
p.minus <- p.all[-(n.int+1)]
p.plus <- p.all[-1]
print(paste("total no. intervals = ",n.int))

## get z.min and z.max for all i (for intervals D_j)
z.min <- rep(-1,n.tests)
z.max <- rep(-1,n.tests)
numalloc.poss <- 1
for(i in 1:n.tests){
if(q.myuni){
z.min[i] <- myMatch(pprev[i],p.minus,tol=tol)
if(pvals[i]==0) z.max[i] <- 1
else z.max[i] <- myMatch(pvals[i],p.plus,tol=tol)
}
else{
z.min[i] <- match(pprev[i],p.minus)
if(pvals[i]==0) z.max[i] <- 1
else z.max[i] <- match(pvals[i],p.plus)
}
numalloc.poss <- numalloc.poss * (z.max[i]-z.min[i]+1)
}
print(paste("total no. possible alloc. = ",numalloc.poss))
#print(pvals)
#print(pprev)
#print(z.min)
#print(z.max)

## get global r.minus,r.plus (over all alloc)
r.minus <- rep(-1,n.int)
r.plus <- rep(-1,n.int)
for(j in 1:n.int){
r.minus[j] <- sum(z.max<j)+1
r.plus[j] <- sum(z.min<=j)
#### this bit was wrong !!!
####dum <- z.max>=j & z.min<=j
####r.minus[j] <- min((1:n.tests)[dum])
####r.plus[j] <- max((1:n.tests)[dum])
}

## get sf,sc
a.minus <- r.minus*alpha/n.tests
a.plus <- r.plus*alpha/n.tests
dummy <- as.numeric(p.minus <= a.plus)
global.sf <- sum(sign(rev(cumsum(rev(dummy)))))
dummy <- as.numeric(p.plus <= a.minus)
global.sc <- sum(sign(rev(cumsum(rev(dummy)))))
print(paste("global sf = ",global.sf))
print(paste("global sc = ",global.sc))

## print for intervals D_j
leng <- r.plus-r.minus+1
#ytable <- cbind(p.minus,p.plus,r.minus,r.plus,leng,a.minus,a.plus)
ytable <- cbind(round(p.minus,d=4),round(p.plus,d=4),r.minus,r.plus,
leng,round(a.minus,d=4),round(a.plus,d=4))
ytable <- data.frame(ytable)
names(ytable) <- c("p.minus","p.plus","r.minus","r.plus","leng","a.minus","a.plus")
print(ytable)

## redefine intervals rej,fuzzy,acc
new.n.int <- global.sf-global.sc + as.numeric(global.sc>0) + as.numeric(global.sf<n.int)
print(paste("reduced no. intervals = ",new.n.int))
new.p.minus <- rep(-1,new.n.int)
new.p.plus <- rep(-1,new.n.int)
if(global.sc>0){
new.p.minus[1] <- p.minus[1]
new.p.plus[1] <- p.plus[global.sc]
}
if(global.sf>global.sc){
for(j in 1:(global.sf-global.sc)){
new.p.minus[j+as.numeric(global.sc>0)] <- p.minus[global.sc+j]
new.p.plus[j+as.numeric(global.sc>0)] <- p.plus[global.sc+j]
}}
if(global.sf<n.int){
new.p.minus[new.n.int] <- p.minus[global.sf+1]
new.p.plus[new.n.int] <- p.plus[n.int]
}

## get z.min and z.max for all i (for new intervals)
new.z.min <- rep(-1,n.tests)
new.z.max <- rep(-1,n.tests)
new.probsample <- array(-1,dim=c(n.tests,new.n.int))
new.numalloc.poss <- 1
for(i in 1:n.tests){
if(z.min[i]<=global.sc) new.z.min[i] <- 1
else if(z.min[i]<=global.sf & global.sc>0) new.z.min[i] <- z.min[i]-global.sc+1
else if(z.min[i]<=global.sf & global.sc==0) new.z.min[i] <- z.min[i]
else new.z.min[i] <- new.n.int
if(z.max[i]<=global.sc) new.z.max[i] <- 1
else if(z.max[i]<=global.sf & global.sc>0) new.z.max[i] <- z.max[i]-global.sc+1
else if(z.max[i]<=global.sf & global.sc==0) new.z.max[i] <- z.max[i]
else new.z.max[i] <- new.n.int
for(j in new.z.min[i]:new.z.max[i]){
new.probsample[i,j] <- max(0, (min(new.p.plus[j],pvals[i])-max(new.p.minus[j],pprev[i])) )/(pvals[i]-pprev[i])
}
#print(paste("new.probsample = ",new.probsample[i,]))
new.numalloc.poss <- new.numalloc.poss * (new.z.max[i]-new.z.min[i]+1)
}
print(paste("reduced no. alloc. = ",new.numalloc.poss))

## get global r.minus,r.plus (over all alloc) for new intervals
new.r.minus <- rep(-1,new.n.int)
new.r.plus <- rep(-1,new.n.int)
for(j in 1:new.n.int){
new.r.minus[j] <- sum(new.z.max<j)+1
new.r.plus[j] <- sum(new.z.min<=j)
#### this bit was wrong !!!
####dum <- new.z.max>=j & new.z.min<=j
####new.r.minus[j] <- min((1:n.tests)[dum])
####new.r.plus[j] <- max((1:n.tests)[dum])
}

## print for new intervals
#ztable <- cbind(new.p.minus,new.p.plus,new.r.minus,new.r.plus)
ztable <- cbind(round(new.p.minus,d=4),round(new.p.plus,d=4),new.r.minus,new.r.plus)
ztable <- data.frame(ztable)
names(ztable) <- c("new.p.minus","new.p.plus","new.r.minus","new.r.plus")
print(ztable)

## get the allocations matrices (for prob.rej.ties)
if(global.sc>0 & global.sf<n.int) len.max <- max(new.r.plus[-c(1,new.n.int)]-new.r.minus[-c(1,new.n.int)]+1)
else if(global.sc>0 & global.sf==n.int) len.max <- max(new.r.plus[-1]-new.r.minus[-1]+1)
else if(global.sc==0 & global.sf<n.int) len.max <- max(new.r.plus[-new.n.int]-new.r.minus[-new.n.int]+1)
else len.max <- max(new.r.plus-new.r.minus+1)
#print(paste("len.max = ",len.max))
if(len.max>2) getAlloc(len.max-1)

############################ enumerate allocations

lowerprod <- rep(1,n.tests)
lowerprod[1] <- 1
for(i in 2:n.tests){
delt <- new.z.max[i-1]-new.z.min[i-1]+1
lowerprod[i] <- lowerprod[i-1] * delt
}
#print("new total alloc = ")
#print(lowerprod[n.tests]*(new.z.max[n.tests]-new.z.min[n.tests]+1))

z <- array(-1,dim=c(new.numalloc.poss,n.tests))
for(i in 1:n.tests){
delt <- new.z.max[i]-new.z.min[i]+1
upperprod <- new.numalloc.poss/(lowerprod[i]*delt)
##zprime <- rep(-1,upperprod)
zprime <- rep(-1,upperprod*delt)
for(j in 1:delt){
zprime[ (j-1)*upperprod + (1:upperprod) ] <- new.z.min[i]+j-1
}
z[,i] <- rep(zprime,lowerprod[i])
}
#print(z)

############################## loop over allocations

## initialise/declare vectors
len <- rep(-1,new.n.int)
r.minus <- rep(-1,new.n.int)
r.plus <- rep(-1,new.n.int)
tau <- rep(0,n.tests)
totalpalloc <- 0

print("starting loop over allocations")

for(d in 1:new.numalloc.poss){

#print(d)
#print(z[d,])
#for(i in 1:n.tests){
# print(new.probsample[i,z[d,i]])
#}

## prob alloc
logpalloc <- 0
for(i in 1:n.tests){
logpalloc <- logpalloc + log(new.probsample[i,z[d,i]])
}
palloc <- exp(logpalloc)

#print(palloc)

## normalising constant
totalpalloc <- totalpalloc + palloc

## ranks, length ties for each interval
len[1] <- sum(z[d,]==1)
r.minus[1] <- 1
if(new.n.int>1){
for(j in 2:new.n.int){
len[j] <- sum(z[d,]==j)
r.minus[j] <- r.minus[j-1] + len[j-1]
}
}
r.plus <- c(r.minus[-1]-1,n.tests)

#print(len)
#print(r.minus)
#print(r.plus)
#print(new.p.minus)
#print(new.p.plus)

##### CHANGED CALC. OF SF HERE
## get sf, sc (note def. of sc different from global one)
dummy <- as.numeric(new.p.minus < r.plus*alpha/n.tests)
if(global.sf<n.int)
sf <- sum(sign(rev(cumsum(rev(dummy[-new.n.int])))))
else sf <- sum(sign(rev(cumsum(rev(dummy)))))
dummy <- as.numeric(new.p.plus <= r.plus*alpha/n.tests)
sc <- sum(sign(rev(cumsum(rev(dummy)))))
#print(paste("sf = ",sf))
#print(paste("sc = ",sc))

## get pi and eta for each interval

probrej <- rep(-1,new.n.int)
tzero <- rep(-1,new.n.int+1)
tzero[new.n.int+1] <- 1

if(sf<new.n.int){
for(j in new.n.int:(sf+1)){
probrej[j] <- 0
tzero[j] <- 1
}
}
if(sf>sc){
for(j in sf:(sc+1)){
if(len[j]==0){
tzero[j] <- tzero[j+1]
}
else{
dummy <- probRejTies(new.p.minus[j],new.p.plus[j],
r.minus=r.minus[j],len=len[j],
ntests=n.tests,alpha=alpha)
probrej[j] <- 1 - tzero[j+1]*(1 - dummy\$piprob)
tzero[j] <- dummy\$tzero * tzero[j+1]
}
}
}
if(sc>0){
for(j in sc:1){
probrej[j] <- 1
tzero[j] <- 0
}
}

## get tau for each i
for(i in 1:n.tests){
tau[i] <- tau[i] + palloc*probrej[z[d,i]]
}
}

#print("total palloc = ")
#print(totalpalloc)

## normalise tau with totalprob
## (since not summing over all possible alloc)
#tau <- tau/totalpalloc

print("")
print("Exact Method")
print(paste("alpha = ",alpha))
#xtable <- cbind(pvals,pprev,z.min,z.max,new.z.min,new.z.max,tau)
#print(xtable)

ind <- order(pvals,pprev)
xtable <- cbind(round(pvals,d=4),round(pprev,d=4),z.min,z.max,new.z.min,new.z.max,round(tau,d=4))
xtable <- data.frame(xtable)
names(xtable) <- c("pvals","pprev","z.min","z.max","new.z.min","new.z.max","tau")
print(xtable[ind,])

return(xtable)

}
```

## Try the fuzzyFDR package in your browser

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

fuzzyFDR documentation built on May 2, 2019, 5:14 a.m.