R/qwmwhd.R

qwmwhd <-
function(x,y,q=seq(5,40,5)/100,xlab="Quantile",ylab="Sum of q and 1-q Quantiles",plotit=TRUE,alpha=.05,nboot=1000){
#
#  Plot that provides perspective on the degree a distribution is symmetric about zero.
#  This function plots the sum of q and 1-q quantiles of the distribution of D=X-Y, X and Y independent.
#  A 1-alpha confidence interval for the sum is indicated by a +
#  If the distribution is symmetric
#  the plot should be approximately a horizontal line.
#
#  FWE is controlled via Hochberg's method, which was used to determine critical
#  p-values based on the argument
#  alpha.
#
#  Can alter the quantiles compared via the argument
#  q
#  q must be less than .5
#
x=elimna(x)
y=elimna(y)
n1=length(x)
n2=length(y)
output=matrix(NA,ncol=8,nrow=length(q))
dimnames(output)=list(NULL,c("quantile","Est.1","Est.2","SUM","ci.low","ci.up","p_crit","p-value"))
for(i in 1:length(q)){
test=cbmhd(x,y,q=q[i],plotit=FALSE,nboot=nboot)
output[i,1]=q[i]
output[i,2]=test$Est1
output[i,3]=test$Est2
output[i,4]=test$sum
output[i,8]=test$p.value
output[i,5]=test$ci[1]
output[i,6]=test$ci[2]
}
temp=order(output[,8],decreasing=TRUE)
zvec=alpha/c(1:length(q))
output[temp,7]=zvec
output <- data.frame(output)
output$signif=rep("YES",nrow(output))
for(i in 1:nrow(output)){
if(output[temp[i],8]>output[temp[i],7])output$signif[temp[i]]="NO"
if(output[temp[i],8]<=output[temp[i],7])break
}
if(plotit){
plot(rep(q,3),c(output[,4],output[,5],output[,6]),type="n",xlab=xlab,ylab=ylab)
points(q,output[,6],pch="+")
points(q,output[,5],pch="+")
points(q,output[,4],pch="*")
}
list(n=c(n1,n2),output=output)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.