R/ancov2COV.R

ancov2COV <-
function(x1,y1,x2,y2,tr=.2,test=yuen,cr=NULL,pr=TRUE,DETAILS=FALSE,cp.value=FALSE,
plotit=FALSE,xlab='X',ylab='Y',zlab=NULL,span=.75,PV=TRUE,FRAC=.5,MC=FALSE,q=.5,
iter=1000,alpha=.05,TPM=FALSE,tau=.05,est=tmean,fr=1,...){
#
# ANCOVA two covariates, no parametric assumption about the regression surface.
# Use all design points nested deeply within the cloud of data.
# Global test statistic is the average of the p-values based on 
# Yuen's test performed at each of the deepest 
# points where comparisons can be made.
#
# TPM=TRUE: replace average of the p-values with the test statistic
#  studied by 
#  Zaykin, D. V., Zhivotovsky, L. A., Westfall, P. H., & Weir, B.S. (2002).
#   Truncated product method for combining p-values.
#    Genetic Epidemiology 22, 170--185.
#
#  x1 and x2 assumed to be a matrix or data frame with two columns
#
#   if plotit=TRUE then if
#   PV=TRUE create a plot of the p.values as a function of the two covariates
#           using LOESS.
#   if PV=FALSE, plot the difference between the dependent variable as a function of
#         the covariates
#
#   By default, Yuen's test is used, but other tests can be used via the argument
#   test
#
#   pr=TRUE: warning messages will be printed
#
#   DETAILS=TRUE: all p.values are reported for all covariate points used. 
#
#   span: the span used by LOESS

#   fr: span used by rung3hat when estimating the difference between
#       predicted Y for group 1 minus predicted Y for group 2.
#
#   FRAC is the fraction of least deep covariate points that are ignored
#
#   MC=TRUE: use a multicore processor to compute a critical value and global p.value
#
#   com.p.value=TRUE: compute p.value based on the global hypothesis of no differences.
#
#   iter=1000: number of iterations used to compute a critical value or global p.value
#
#      Function returns:
#   test.stat: the test statistic
#   crit.p.val:  the critical value, reject if test.stat<=crit.p.val
#   min.p.val.point:  the values of the covariate that had the smallest p-value
#   min.p.value:  the minimum p-value among all p-values that were computed.
#
com.p.value=cp.value
if(pr)print('Reject if test.stat is less than or equal to crit.value')
if(FRAC<=0 || FRAC >=1)stop('FRAC should be a value between 0 and 1.')
if(ncol(x1)!=2)stop('Should have two covariates')
xy1=elimna(cbind(x1,y1))
x1=xy1[,1:2]
y1=xy1[,3]
xy2=elimna(cbind(x2,y2))
x2=xy2[,1:2]
y2=xy2[,3]
n=min(c(nrow(x1),nrow(x2)))
if(n<50){
if(pr)print('Warning: sample size is less than 50; critical value unknown')
}
if(is.null(cr)){
if(n>=50 & n<=80)cr=as.vector(regYhat(c(50,75),c(.23,.264),xr=n))
if(n>80)cr=.27
}
flag0=is.null(cr)
flag1=FRAC!=.5
flag3=flag0+flag1+com.p.value+TPM
if(flag3>0){
comp.pv=anc2COV.CV(nrow(x1),nrow(x2),iter=iter,MC=MC,TPM=TPM,tau=tau)
MV=sort(comp.pv$M)
ic=round(alpha*iter)
cr=MV[ic]
}
DONE=FALSE
if(identical(test,qcomhd)){
val=ancovampG(x1,y1,x2,y2,DH=TRUE,SEED=TRUE,test=qcomhd,q=q,FRAC=FRAC)
DONE=TRUE
}
if(!DONE){
if(identical(test,yuen))val=ancovampG(x1,y1,x2,y2,DH=TRUE,SEED=TRUE,test=yuen,tr=tr,FRAC=FRAC)
if(!identical(test,yuen))val=ancovampG(x1,y1,x2,y2,DH=TRUE,SEED=TRUE,test=test,FRAC=FRAC,...)
}
est.dif=rung3hat(x1,y1,fr=fr,pts=val$points,est=est)$rmd-rung3hat(x2,y2,pts=val$points,fr=fr,est=est)$rmd
pavg=mean(val$results[,3])
if(TPM){
vals=val$results[,3]
vals=elimna(vals)
pavg=sum(log(vals[vals<=tau]))
}
mpv=which(val$results[,3]==min(val$results[,3]))
points=val$points
results=val$results
rem.res=results[mpv,3]
rem.points=points[mpv,]
points=cbind(points,est.dif)
dimnames(points)=list(NULL,c('COV 1','COV 2','EST.DIF'))
if(plotit){
if(is.null(zlab)){
if(PV)zlab='P-Value'
if(!PV)zlab='Est.Dif'
}
if(PV)lplot(points[,1:2],results[,3],xlab=xlab,ylab=ylab,zlab=zlab,tick='det',span=span)
if(!PV)lplot(points[,1:2],est.dif,xlab=xlab,ylab=ylab,zlab=zlab,tick='det',span=span)
}
nk=nrow(points)
if(!DETAILS){
points=NULL
results=NULL
}
pval=NULL
if(com.p.value || TPM)pval=1-mean(pavg<=comp.pv$M) 
list(num.points.used=nk,test.stat=pavg,crit.value=cr,GLOBAL.p.value=pval,min.p.val.point=rem.points,min.p.value=rem.res,all.points.used=points,all.results=results[,1:3])
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.