# R/binband.R In WRS2: A Collection of Robust Statistical Methods

#### Documented in binband

binband <- function(x, y, KMS = FALSE, alpha = 0.05, ADJ.P = FALSE, ...){
#
#  Comparing two independent variables in terms of their probability function.
#  For each value that occurs, say x, test P(X=x)=P(Y=x)
#  So this method is useful when dealing with highly discrete data.
#
#  If KMS=T, use Kulinskaya, Morgenthaler and Staudte (2010)
#   method for comparing binomials
# Kulinskaya, E., Morgenthaler, S. and Staudte, R. (2010).
# Variance Stabilizing the Difference of two Binomial
#  Proportions. {\em American Statistician, 64},
#  350--356 DOI:10.1198/tast.2010.09096

#  Otherwise use Storer and Kim.
#
#   ADJ.P=T means that critical p-value is adjusted to control FWE when the sample
#   size is small (<50).
#
#
#  Hochberg's method is used to determine critical p-values so that FWE=alpha
#
cl <- match.call()
x=elimna(x)
y=elimna(y)
vals=sort(unique(c(x,y)))
ncon=length(vals)
n1=length(x)
n2=length(y)
p.values=NA
cv=1
if(!KMS){
output=matrix(NA,ncol=6,nrow=length(vals))
dimnames(output)=list(NULL,c("Value","p1.est","p2.est","p1-p2","p.value","p.crit"))
}
if(KMS){
output=matrix(NA,ncol=8,nrow=length(vals))
dimnames(output)=list(NULL,c("Value","p1.est","p2.est","p1-p2","ci.low","ci.up","p.value",
"p.crit"))
}
for(i in 1:length(vals)){
x1=sum(x==vals[i])
y1=sum(y==vals[i])
if(!KMS){
output[i,5]=twobinom(x1,n1,y1,n2)$p.value output[i,2]=x1/n1 output[i,3]=y1/n2 output[i,1]=vals[i] output[i,4]=output[i,2]-output[i,3] } if(KMS){ temp=bi2KMSv2(x1,n1,y1,n2) output[i,1]=vals[i] output[i,5]=temp$ci[1]
output[i,6]=temp$ci[2] output[i,2]=x1/n1 output[i,3]=y1/n2 output[i,4]=output[i,2]-output[i,3] output[i,7]=temp$p.value
}}
# Determine adjusted  critical p-value using Hochberg method
ncon=length(vals)
dvec=alpha/c(1:ncon)
mn=max(c(n1,n2))
cv=1
if(ncon!=2){
if(mn>50){
cv=2-(mn-50)/50
if(cv<1)cv=1
}
if(mn<=50)cv=2
}
if(KMS){
flag=(output[,7]<=2*alpha)
output[flag,8]=output[flag,8]/cv
}
if(!KMS){
cv=1
flag=(output[,5]<=2*alpha)
if(min(c(n1,n2))<20 && n1!=n2 && ncon>=5)cv=2
output[flag,5]=output[flag,5]/cv
}}
if(KMS){
temp2=order(0-output[,7])
output[temp2,8]=dvec
}
if(!KMS){
temp2=order(0-output[,5])
output[temp2,6]=dvec
}
output

outtable <- data.frame(output)
colnames(outtable)[4] <- "p1-p2"

result <- list(partable = outtable, alpha = alpha, call = cl)
class(result) <- "robtab"
result
}


## Try the WRS2 package in your browser

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

WRS2 documentation built on March 19, 2024, 3:08 a.m.