R/ancovamp.R

ancovamp <-
function(x1,y1,x2,y2,fr1=1,fr2=1,tr=.2,alpha=.05,pts=NA,SEED=T){
#
# Compare two independent  groups using the ancova method.
# No parametric assumption is made about the form of
# the regression lines--a running interval smoother is used.
# Design points are chosen based on depth of points in x1 if pts=NA
#  Assume data are in x1 y1 x2 and y2
#  x1 and x2 should be matrices with two or more columns.
#
if(SEED)set.seed(2) # now cov.mve always returns same result
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.na(pts[1])){
x1<-as.matrix(x1)
pts<-ancdes(x1)
}
pts<-as.matrix(pts)
if(nrow(pts)>=29){
print("WARNING: More than 28 design points")
print("Only first 28 are used.")
pts<-pts[1:28,]
}
n1<-1
n2<-1
vecn<-1
mval1<-cov.mve(x1)
mval2<-cov.mve(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(T,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))
if(sum(flag)==0)stop("No comparable design points found, might increase span.")
mat<-matrix(NA,nrow(pts),8)
dimnames(mat)<-list(NULL,c("n1","n2","DIF","TEST","se","ci.low","ci.hi","p.value"))
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)]
test<-yuen(g1,g2,tr=tr)
mat[i,1]<-length(g1)
mat[i,2]<-length(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,]))
mat[i,3]<-test$dif
mat[i,4]<-test$teststat
mat[i,5]<-test$se
mat[i,8]<-test$p.value
if(nrow(pts)>=2)critv<-smmcrit(test$df,nrow(pts))
if(nrow(pts)==1)critv<-qt(.975,test$df)
cilow<-test$dif-critv*test$se
cihi<-test$dif+critv*test$se
mat[i,6]<-cilow
mat[i,7]<-cihi
}
list(points=pts,output=mat,crit=critv)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.