R/Bband.R

Bband <-
function(x,alpha=.05,plotit=TRUE,sm=TRUE,SEED=TRUE,nboot=500,grp=c(1:4),
xlab="X (First Level)",ylab="Delta",crit=NA,print.all=FALSE,plot.op=FALSE){
#
# Apply the shift function when analyzing main effect in a
# 2 by 2 design.
#
# For variables x1, x2, x3 and x4,
# In effect, this function applies a shift function to the distributions
# d1=(x1+x3)/2 and d2=(x2+x4)/2.
# That is, focus on main effects of Factor B.
#
# grp indicates the groups to be compared. By default grp=c(1,2,3,4)
# meaning that the first level of factor A consists of groups 1 and 2
# and the 2nd level of factor A consists of groups 3 and 4.
# (So level 1 of factor B consists of groups 1 and 3
#
# print.all=F,
# returns number sig, meaning number of confidence intervals that do not
# contain zero,
# the critical value used as well as the KS test statistics.
# print.all=T reports all confidence intervals, the number of which can
# be large.
#
if(!is.list(x) && !is.matrix(x))stop("store data in list mode or a matrix")
if(SEED)set.seed(2)
if(is.matrix(x))x<-listm(x)
for(j in 1:length(x))x[[j]]=elimna(x[[j]])/2
if(length(x)<4)stop("There must be at least 4 groups")
if(length(grp)!=4)stop("The argument grp must have 4 values")
x<-x[grp]
grp=c(1,3,2,4)
x<-x[grp]  # Arrange groups for main effects on factor B
n<-c(length(x[[1]]),length(x[[2]]),length(x[[3]]),length(x[[4]]))
# Approximate the critical value
#
vals<-NA
y<-list()
if(is.na(crit)){
print("Approximating critical value. Please wait.")
for(i in 1:nboot){
for(j in 1:4)
y[[j]]<-rnorm(n[j])
temp<-ks.test(outer(y[[1]],y[[2]],FUN="+"),outer(y[[3]],y[[4]],FUN="+"))
vals[i]<-temp[1]$statistic
}
vals<-sort(vals)
ic<-(1-alpha)*nboot
crit<-vals[ic]
}
if(plot.op){
plotit<-F
g2plot(v1,v2)
}
output<-sband(outer(x[[1]],x[[2]],FUN="+"),outer(x[[3]],x[[4]],FUN="+"),
plotit=plotit,crit=crit,flag=FALSE,sm=sm,xlab=xlab,ylab=ylab)
if(!print.all){
numsig<-output$numsig
ks.test.stat<-ks.test(outer(x[[1]],x[[2]],FUN="+"),
outer(x[[3]],x[[4]],FUN="+"))$statistic
output<-matrix(c(numsig,crit,ks.test.stat),ncol=1)
dimnames(output)<-list(c("number sig","critical value","KS test statistics"),
NULL)
}
output
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.