R/Dqcomhd.R

Dqcomhd <-
function(x,y,q=c(1:9)/10,nboot=1000,plotit=TRUE,SEED=TRUE,xlab="Group 1",ylab="Est.1-Est.2",na.rm=TRUE,alpha=.05){
#
# Compare the quantiles of the marginal distributions associated with  two dependent groups
# via hd estimator. Tied values are allowed.
# When comparing lower or upper quartiles, both power and the probability of Type I error
# compare well to other methods have been derived.
#
#  x: data for group 1
#  y: data for group 2
#  q: the quantiles to be compared
#  nboot: Number of bootstrap samples
#
#
if(SEED)set.seed(2)
if(na.rm){
xy=elimna(cbind(x,y))
x=xy[,1]
y=xy[,2]
}
pv=NULL
output=matrix(NA,nrow=length(q),ncol=10)
dimnames(output)<-list(NULL,c("q","n1","n2","est.1","est.2","est.1_minus_est.2","ci.low","ci.up","p_crit","p-value"))
for(i in 1:length(q)){
output[i,1]=q[i]
output[i,2]=length(elimna(x))
output[i,3]=length(elimna(y))
output[i,4]=hd(x,q=q[i])
output[i,5]=hd(y,q=q[i])
output[i,6]=output[i,4]-output[i,5]
if(na.rm){
temp=bootdpci(x,y,est=hd,q=q[i],dif=FALSE,plotit=FALSE,pr=FALSE,nboot=nboot,alpha=alpha,SEED=FALSE)
output[i,7]=temp$output[1,5]
output[i,8]=temp$output[1,6]
output[i,10]=temp$output[1,3]
}
if(!na.rm){
temp=rmmismcp(x,y,est=hd,q=q[i],plotit=FALSE,pr=FALSE,nboot=nboot,alpha=alpha,SEED=FALSE)
output[i,7]=temp$output[1,6]
output[i,8]=temp$output[1,7]
output[i,10]=temp$output[1,4]
}
}
if(plotit){
xax=rep(output[,4],3)
yax=c(output[,6],output[,7],output[,8])
plot(xax,yax,xlab=xlab,ylab=ylab,type="n")
points(output[,4],output[,6],pch="*")
lines(output[,4],output[,6])
points(output[,4],output[,7],pch="+")
points(output[,4],output[,8],pch="+")
}
temp=order(output[,10],decreasing=TRUE)
zvec=alpha/c(1:length(q))
output[temp,9]=zvec
output <- data.frame(output)
output$signif=rep("YES",nrow(output))
for(i in 1:nrow(output)){
if(output[temp[i],10]>output[temp[i],9])output$signif[temp[i]]="NO"
if(output[temp[i],10]<=output[temp[i],9])break
}
output
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.