R/difQpciMC.R

difQpciMC <-
function(x,y=NULL,q=seq(5,40,5)/100,xlab='Quantile',ylab='Group 1 minus Group 2',plotit=TRUE,
SEED=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. A 1-alpha confidence interval for the sum is indicated by a +
#  If the distributions are symmetric
#  the plot should be approximately a horizontal line. If in addition the median
#  of the difference scores is zero, the horizontal line will intersect the y-axis at zero.
#
#  Similar to difQplot, only plots fewer quantiles by default and returns p-values for
#  each quantile indicated by the argument q.
#
#  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=as.matrix(x)
if(is.null(y))dif=x
if(ncol(x)>2)stop("Should be at most two groups")
if(ncol(x)==2)dif=x[,1]-x[,2]
if(!is.null(y))dif=x-y
dif=elimna(dif)
dif=as.matrix(dif)
nv=length(dif)
output=matrix(NA,ncol=8,nrow=length(q))
dimnames(output)=list(NULL,c('quantile','Est_q','Est_1.minus.q','SUM','ci.low','ci.up','p_crit','p-value'))
for(i in 1:length(q)){
test=DqdifMC(dif,q=q[i],plotit=FALSE,nboot=nboot,SEED=SEED)
output[i,1]=q[i]
output[i,2]=test$est.q
output[i,3]=test$est.1.minus.q
output[i,8]=test$p.value
output[i,5]=test$conf.interval[1]
output[i,6]=test$conf.interval[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
}
output[,4]=output[,2]+output[,3]
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=nv,output=output)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.