permb <- function(xu,yu,xp, yp,r, mu=NULL, method="EH",alternative ="two.sided") { np=5000 method <- pmatch(method, c("EH", "Maritz")) if (is.na(method)) stop("\nbw.method MUST be one of: 'EH', 'Maritz'\n") if (missing(r)) { r <- cor(xp,yp) } else { if ((length(xp) < 1) |(length(yp) < 1)) { stop("\n there no paired observations\n") } } alt <- pmatch(alternative, c("two.sided","less", "greater")) if ( method==1) { # this function computes the mean of the differences for the # unpaired sample # This part compute permutation test for paired data paired <- function(x1,x2) { n=length(x1) y1=numeric(n) y2=numeric(n) for(i in 1:n) { u=runif(1) if ((u-0.5)<=0) { y1[i]=x1[i] y2[i]=x2[i] } else { if((u-0.5)>0) { y1[i]=x2[i] y2[i]=x1[i] } } } axp=mean(y1-y2) return(axp) } unpaired <- function(x1,x2) { m1=length(x1) m2=length(x2) tc=c(x1,x2) n=length(tc) y=sample(tc,replace=F) ax=mean(y[1:m1]) ay=mean(y[(m1+1):n]) # the difference between the sample means cmean=ax-ay return(cmean) } # weight for equal variance permt <- function(xp,yp,xu,yu,r) { n=length(xp) n1=length(xu) n2=length(yu) w=(1/n2+1/n1)/((2-2*r)/n +(1/n2+1/n1)) if(n==0) w=0 if((n1==0)&&(n2==0)) w=1 Tp=paired(xp,yp) Tu=unpaired(xu,yu) t10=w*(Tp) +(1-w)*(Tu) return(t10) } t10=replicate(np,permt(xp,yp,xu,yu,r)) w=(1/n2+1/n1)/((2-2*r)/n +(1/n2+1/n1)) if(n==0) w=0 if((n1==0)&&(n2==0)) w=1 # the value of the observed test statistic Tobs=w*(mean(xp)-mean(yp)) +(1-w)*(mean(xu)-mean(yu)) result=c(t10,Tobs) #The p-value is the chance of seeing a test statistic at least as #large as what was observed (.667). This can be #estimated by the # proportion of permutation statistics that were as #large or larger than the observed one if(alt==1) { pvalue=sum(abs(t10)>=abs(Tobs))/length(result) test="alternative hypothesis: true difference in means is not equal to 0" } if(alt==2) { pvalue=sum(t10<=Tobs)/length(result) test="alternative hypothesis: true difference in means is less than 0" } if(alt==3) { pvalue=sum(t10>=Tobs)/length(result) test="alternative hypothesis: true difference in means is greater than 0" } approach ="Permutation based test on Einsporn and Habtzghi's approach " ans=list(approach, test,"pvalue",pvalue) } if(method==2) { # Maritz's Permutation test permm <- function(xp,yp,xu,yu) { n=length(xp) n1=length(xu) n2=length(yu) # x=c(xu, rep(0,n2),xp) y=c(rep(0,n1),yu,yp) m=n1+n2+n y1=numeric(m) y2=numeric(m) for(i in 1:m) { T2=sample(c(x[i],y[i]),replace=F) y1[i]=T2[1] y2[i]=T2[2] } b=length(y1[y1!=0]) b2=length(y2[y2!=0]) xbar=sum(y1)/b ybar=sum(y2)/b2 diff=xbar-ybar return(diff) } #Maritz tm1=replicate(np,permm(xp,yp,xu,yu)) # this part computes the observed value of Marith n=length(xp) n1=length(xu) n2=length(yu) xl=c(xu, rep(0,n2),xp) yl=c(rep(0,n1),yu,yp) bl=c(rep(1,n1),rep(0,n2),rep(1,n)) #ayp <- numeric(np) xbar1=sum(xl*bl)/sum(bl) ybar1=(sum(xl)+sum(yl)-sum(xl*bl))/(n1+n2+2*n-sum(bl)) Dm0=xbar1-ybar1 resultm=c(tm1,Dm0) if(alt==1) { pvalue=sum(abs(tm1)>=abs(Dm0))/length(resultm) test="alternative hypothesis: true difference in means is not equal to 0" } if(alt==2) {pvalue=sum(tm1<=Dm0)/length(resultm) test="alternative hypothesis: true difference in means is less than 0" } if(alt==3) { pvalue=sum(tm1>=Dm0)/length(resultm) test="alternative hypothesis: true difference in means is greater than 0" } approach="Permutation based test on Maritz's approach " ans=list(approach, test,"pvalue",pvalue) } # the value of the test stat using the permutation test ans } # Example n=20 n1=15 n2=10 r=0.8 xp=rnorm(n) yp=r*xp+(1-r)*(rnorm(n)) xu=rnorm(n1) yu=rnorm(n2) r=0.2 mu=0 permb(xu,yu,xp, yp,r, mu, method="EH", alternative="two.sided")