R/kfwefunctions.R In someKfwer: Controlling the Generalized Familywise Error Rate

Documented in getAlphakfweGRkfweLRkfweOrd

```######################################### Ordinal Procedure
kfweOrd<-function(p,k=1,alpha=.01,ord=NULL,alpha.prime=alpha,J=qnbinom(alpha,k,alpha.prime),disp=TRUE,GD=FALSE){

if(!is.null(ord)){# sort by ord
o <- order(ord,decreasing=T)
}
else{o <- 1:length(p)}
ps <- p[o]

if(GD) alpha1 <- k*alpha.prime/(J+k)
else alpha1 <- alpha.prime

u<-cumsum(ps>alpha1)

if(sum(u<=J)>0){
h <- rep(0,length(p))
h[1:max(which(u<=J))] <- 1
h[ps>alpha1] <- 0

if(sum(h)<k)  h[(h==0)&&(ps<=alpha1)][1:min(sum((h==0)&&(ps<=alpha1)),k-1-sum(h))]=1
h[o] <- h
}
else{h <- rep(0,length(p))}

if(disp==T) 	cat(paste("Ordered k-FWER procedure\n ",length(p)," tests, k=", k, ", alpha=",alpha,", individual alpha threshold=",round(alpha1,digits=7),"\n ",J," jumps allowed","\n ",sum(h)," rejections\n\n",sep=""))
return(h==1)}

################kFWE by Guo e Romano - (Step Down and Single Step) procedure
kfweGR<-function(p,k=1,alpha=.01,disp=TRUE,SD=TRUE,const=10,alpha.prime=getAlpha(k=k,s=length(p),alpha=alpha,const=const)) {

if(is.null(alpha.prime)) alpha.prime=getAlpha(k=k,s=length(p),alpha=alpha,const=50)

rej <- rep(0,length(p))
rej[p<=alpha.prime] <- 1
n.rej <- sum(p<=alpha.prime)

sd=((n.rej>=k)&(SD))
while (sd){
alpha.prime=getAlpha(k=k,s=length(p)-n.rej+k-1,alpha=alpha,const=const)
rej[p<=alpha.prime] <- 1
sd=n.rej<sum(rej)
n.rej=sum(p<=alpha.prime)
}

if(disp) 	cat(paste("Guo and Romano k-FWER ",switch(SD,"Step Down ",""),"procedure\n ",length(p)," tests, k=", k, ", alpha=",alpha, "\n ",round(alpha.prime,digits=7)," individual alpha threshold\n ",n.rej," rejections\n\n",sep=""))
return(rej==1)
}

############kFWE by Lehamann and Romano
kfweLR <- function(p,k=1,alpha=0.01,disp=TRUE) {
s <- length(p)
sdconst <- rep(1,s)

sdconst[1:min(k,s)] <- k*alpha/s
if(s>k) sdconst[(k+1):s] <- k*alpha/(s+k-((k+1):s))

ps <- sort(p)
u <- ps<sdconst
res <- 0
if(any(u)) {
w <- min(which(!u))-1
res <- ps[w]}
p[which(p>res)] <- 1
p[p<=alpha] <- 0
h=(!p)
if(disp) 	cat(paste("Lehmann e Romano k-FWER Step Down procedure\n ",length(p)," tests, k=", k, ", alpha=",alpha, "\n ",sum(h)," rejections\n\n",sep=""))
return(h==1)
}

############### other functions

getAlpha <- function(s,k=1,alpha=.01,const=10){
start<-1E-8
stop<-1

delta<-1
while(delta>1E-8){
alphas<-start+((0:const)/(const))*(stop-start)
temp<-round(pbinom(k-1,s,alphas,lower.tail=F),digits=7)-alpha
temp2<-max(which(temp<=0))
alpha.prime<-alphas[temp2]
start<-alpha.prime
stop<-alphas[temp2+1]
delta<-abs(temp[temp2])
}

return(alpha.prime)
}
```

Try the someKfwer package in your browser

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

someKfwer documentation built on May 29, 2017, 1:39 p.m.