R/ancovampG.R

ancovampG <-
function(x1,y1,x2,y2,fr1=1,fr2=1,est=tmean,alpha=.05,pts=NULL,SEED=TRUE,test=yuen,DH=FALSE,FRAC=.5,cov.fun=skip.cov,pr=FALSE,q=.5,plotit=FALSE,pv=FALSE,theta=50,xlab=' ',ylab=' ',SCAT=FALSE,zlab=' ',...){
#
#  ANCOVA:
#
#  This function generalizes the R function ancovamp 
#  so that any hypothesis testing method
#  can be used to compare groups at specified design points.
#
# No parametric assumption is made about the form of
# the regression surface--a running interval smoother is used.
# Design points are chosen based on depth of points in x1 if pts=NULL
#  Assume data are in x1 y1 x2 and y2, can have more than one covariate
#
#  test: argument test determines the method that will be used to compare groups.
#        Example: test=medpb2 would compare medians using a percentile bootstrap
#
#  pts: a matrix of design points at which groups are compared
#
#  DH=T, groups compared at the deepest (1-FRAC) design points.
#  if DH=T, there are two covariates and plot=TRUE, plot a smooth with dependent variable=p.values if pv=TRUE
#  or the estimated difference in the measures of location if pv=FALSE
#  If SCAT=TRUE, instead create a scatterplot of the points used in pts, the covariate values
#  and mark the significant ones with *
#
#  theta can be use to rotate the plot.
#
#  SEED=TRUE sets the seed for the random number generator 
#       so that same result is always returned when 
#        using a bootstrap method or when using cov.mve or cov.mcd
#
#   cov.fun: returns covariance matrix in $cov (e.g. 
#   skipcov does not return it in $cov, but skip.cov does. So cov.mve could be used)
#
#   Returns:
#   designs points where comparisons were made.
#   n's used, p-values
#   crit.p.value: critical p-value based on Hochberg's method for controlling FWE
#   sig=1 if a signficant result based on Hochberg; 0 otherwise
#
x1=as.matrix(x1)
p=ncol(x1)
p1=p+1
m1=elimna(cbind(x1,y1))
x1=m1[,1:p]
y1=m1[,p1]
x2=as.matrix(x2)
p=ncol(x2)
p1=p+1
m2=elimna(cbind(x2,y2))
x2=m2[,1:p]
y2=m2[,p1]
#
#
#
if(is.null(pts[1])){
x1<-as.matrix(x1)
pts<-ancdes(x1,DH=DH,FRAC=FRAC)
}
pts<-as.matrix(pts)
n1<-1
n2<-1
vecn<-1
mval1<-cov.fun(x1)
mval2<-cov.fun(x2)
for(i in 1:nrow(pts)){
n1[i]<-length(y1[near3d(x1,pts[i,],fr1,mval1)])
n2[i]<-length(y2[near3d(x2,pts[i,],fr2,mval2)])
}
flag<-rep(TRUE,nrow(pts))
for(i in 1:nrow(pts))if(n1[i]<10 || n2[i]<10)flag[i]<-F
flag=as.logical(flag)
pts<-pts[flag,]
if(sum(flag)==1)pts<-t(as.matrix(pts))
dd=NULL
if(sum(flag)==0){
print('No comparable design points found, might increase span.')
pts=NULL
mat=NULL
dd=NULL
}
if(sum(flag)>0){
mat<-matrix(NA,nrow(pts),6)
mat[,5]=0
dimnames(mat)<-list(NULL,c('n1','n2','p.value','crit.p.value','Sig','est.dif'))
output=list()
for (i in 1:nrow(pts)){
g1<-y1[near3d(x1,pts[i,],fr1,mval1)]
g2<-y2[near3d(x2,pts[i,],fr2,mval2)]
g1<-g1[!is.na(g1)]
g2<-g2[!is.na(g2)]
if(identical(test,qcomhd))temp=qcomhd(g1,g2,q=q,plotit=FALSE)
if(!identical(test,qcomhd))temp=test(g1,g2,...)
if(is.null(temp$p.value))print('Apparently argument test is a function that does not return a p-value')
mat[i,3]=temp$p.value
output[[i]]=temp
mat[i,1]<-length(g1)
mat[i,2]<-length(g2)
mat[i,6]=est(g1)-est(g2)
if(length(g1)<=5)print(paste('Warning, there are',length(g1),' points corresponding to the design point X=',pts[i,]))
if(length(g2)<=5)print(paste('Warning, there are',length(g2),' points corresponding to the design point X=',pts[i,]))
}
npt=nrow(pts)
dvec=alpha/c(1:npt)
temp2<-order(0-mat[,3])
sigvec<-(mat[temp2,3]>=dvec)
dd=0
if(sum(sigvec)<npt)dd<-npt-sum(sigvec) #number that are sig.
mat[temp2,4]=dvec
flag=mat[,3]<=mat[,4]
if(sum(flag)>0)mat[flag,5]=1
}
if(SCAT)plotit=FALSE
if(DH){
if(ncol(x1)==2){
if(plotit){
if(pv)lplot(pts,mat[,3],xlab=xlab,ylab=ylab,ticktype='det',zlab=zlab,theta=theta)
if(!pv)lplot(pts,mat[,6],xlab=xlab,ylab=ylab,ticktype='det',zlab=zlab,theta=theta)
}
if(SCAT){
chk=mat[,5]==1
plot(pts[!chk,1],pts[!chk,2],xlab=xlab,ylab=ylab)
points(pts[chk,1],pts[chk,2],pch='*')
}
}}
list(points=pts,results=mat,num.sig=dd)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.