R/twoKlin.R

twoKlin <-
function(x=NULL,x1=NULL,x2=NULL,tr=.2,alpha=.05,pr=TRUE,opt=1){
#
#  A step-down MCP based on K independent tests.
#  Use Fisher method based on p-values coupled with Hochberg
# 
# Data are assumed to be stored in two R variables, x1 and x2 or in one 
#  R variable, x
#
# If stored in x1 and x2, they are assumed to be matrices with K columns
# or to have list mode, both having length K.
#
# If the data are stored in x, 
# x is assumed to have 2K columns if a matrix or length 2K if it has list mode.
#
# If data are stored in x1 and x2, for each column, compute a p-value.
# That is, perform a test based on the data in column 1 of x1 and x2, 
# followed by a test using the data in column 2 of x1 and x2, etc.
#
# If data are stored in x, the first test is based 
# on the data in columns 1 and K+1,
# the second test is based on columns 2 and K+2, etc. 
#
#  opt=1  Fisher's method
#  opt=2  Chen-Nadarajah method
#  opt=3  Max method
#
if(is.null(x[1])){
if(is.matrix(x1))x=cbind(x1,x2)
if(is.list(x1))x=c(x1,x2)
}
if(is.matrix(x))x=listm(x)
crit=NA
n1=NA
n2=NA
if(is.matrix(x) || is.data.frame(x))K2=ncol(x)
if(is.list(x))K2=length(x)
K=floor(K2/2)
if(2*K!=K2)stop('Total number of groups, K2, should be an even number') 
ic=0
ic2=K
pv=NULL
for(i in 1:K){
ic=ic+1
ic2=ic2+1
testit=yuen(x[[ic]],x[[ic2]],tr=tr,alpha=alpha)
n1[ic]=testit$n1
n2[ic]=testit$n2
pv[ic]=testit$p.value 
}  
pick=NULL
v=order(pv)
ic=0
for(i in K:1){
K2=2*K
flag=TRUE
if(opt==1){
i2=i*2
if(i==K)res=(0-2)*sum(log(pv))  # Fisher test statistic
if(i<K)res=(0-2)*sum(log(pv[-pick]))  # Fisher test statistic
pvF=1-pchisq(res,i2)   #Fisher p-value based on all tests.
}
if(opt==2){
if(i==K)res=sum(qnorm(pv/2)^2)  # C-N test
if(i<K)res=sum(qnorm(pv[-pick]/2)^2)
pvF=1-pchisq(res,i)
}
if(opt==3){
if(i==K)res=max(pv)
if(i<K)res=max(pv[-pick])
pvF=pbeta(res,i,1)
}
if(pvF>alpha)flag=TRUE
if(pvF<=alpha/(K+1-i)){
ic=ic+1
pick=c(pick,v[ic])
flag=FALSE
if(pv[v[ic]]>alpha)flag=TRUE
}
if(flag)break
}
Decision=rep('Not Sig',length(pv))
if(!is.null(pick))Decision[pick]='Reject'
nsig=sum(length(pick))
list(n1=n1,n2=n2,p.values=pv,
Decisions=as.matrix(Decision),num.sig=nsig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.