R/ancova.R

ancova <-
function(x1,y1,x2,y2,fr1=1,fr2=1,tr=.2,alpha=.05,plotit=TRUE,pts=NA,sm=FALSE,
pr=TRUE,xout=FALSE,outfun=out,LP=TRUE,SCAT=TRUE,xlab='X',ylab='Y',pch1='*',pch2='+',
skip.crit=FALSE,crit.val=1.09,...){
#
# 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.
#
#  Assume data are in x1 y1 x2 and y2
#
#  sm=T will create smooths using bootstrap bagging.
#  pts can be used to specify the design points where the regression lines
#  are to be compared.
#
#   skip.crit: this argument is used to avoid computational issues associated with ancdet.pv
#
if(ncol(as.matrix(x1))>1)stop('One covariate only is allowed with this function')
if(length(x1)!=length(y1))stop("x1 and y1 have different lengths")
if(length(x2)!=length(y2))stop("x2 and y2 have different lengths")
xy=elimna(cbind(x1,y1))
x1=xy[,1]
y1=xy[,2]
xy=elimna(cbind(x2,y2))
x2=xy[,1]
y2=xy[,2]
if(pr){
print("NOTE: Confidence intervals are adjusted to control the probability")
print("of at least one Type I error.")
print("But p-values are not")
}
if(xout){
flag<-outfun(x1,...)$keep
x1<-x1[flag]
y1<-y1[flag]
flag<-outfun(x2,...)$keep
x2<-x2[flag]
y2<-y2[flag]
}
if(is.na(pts[1])){
npt<-5
isub<-c(1:5)  # Initialize isub
test<-c(1:5)
xorder<-order(x1)
y1<-y1[xorder]
x1<-x1[xorder]
xorder<-order(x2)
y2<-y2[xorder]
x2<-x2[xorder]
n1<-1
n2<-1
vecn<-1
for(i in 1:length(x1))n1[i]<-length(y1[near(x1,x1[i],fr1)])
for(i in 1:length(x1))n2[i]<-length(y2[near(x2,x1[i],fr2)])
for(i in 1:length(x1))vecn[i]<-min(n1[i],n2[i])
sub<-c(1:length(x1))
isub[1]<-min(sub[vecn>=12])
isub[5]<-max(sub[vecn>=12])
isub[3]<-floor((isub[1]+isub[5])/2)
isub[2]<-floor((isub[1]+isub[3])/2)
isub[4]<-floor((isub[3]+isub[5])/2)
mat<-matrix(NA,5,10)
dimnames(mat)<-list(NULL,c("X","n1","n2","DIF","TEST","se","ci.low","ci.hi","p.value","crit.val"))
for (i in 1:5){
g1<-y1[near(x1,x1[isub[i]],fr1)]
g2<-y2[near(x2,x1[isub[i]],fr2)]
g1<-g1[!is.na(g1)]
g2<-g2[!is.na(g2)]
test<-yuen(g1,g2,tr=tr)
mat[i,1]<-x1[isub[i]]
mat[i,2]<-length(g1)
mat[i,3]<-length(g2)
mat[i,4]<-test$dif
mat[i,5]<-test$teststat
mat[i,6]<-test$se
if(skip.crit)critv=crit.val
if(!skip.crit){
critv<-NA
if(alpha==.05)critv<-smmcrit(test$df,5)
if(alpha==.01)critv<-smmcrit01(test$df,5)
if(is.na(critv))critv<-smmval(test$df,5,alpha=alpha)
mat[i,10]<-critv
}
cilow<-test$dif-critv*test$se
cihi<-test$dif+critv*test$se
mat[i,7]<-cilow
mat[i,8]<-cihi
mat[i,9]<-test$p.value
}}
if(!is.na(pts[1])){
if(!skip.crit){
if(length(pts)>=29)stop("At most 28 points can be compared")
}
n1<-1
n2<-1
vecn<-1
for(i in 1:length(pts)){
n1[i]<-length(y1[near(x1,pts[i],fr1)])
n2[i]<-length(y2[near(x2,pts[i],fr2)])
}
mat<-matrix(NA,length(pts),10)
dimnames(mat)<-list(NULL,c("X","n1","n2","DIF","TEST","se","ci.low","ci.hi",
"p.value","crit.val"))
for (i in 1:length(pts)){
g1<-y1[near(x1,pts[i],fr1)]
g2<-y2[near(x2,pts[i],fr2)]
g1<-g1[!is.na(g1)]
g2<-g2[!is.na(g2)]
test<-yuen(g1,g2,tr=tr)
mat[i,1]<-pts[i]
mat[i,2]<-length(g1)
mat[i,3]<-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,4]<-test$dif
mat[i,5]<-test$teststat
mat[i,6]<-test$se
if(skip.crit)critv=crit.val
if(!skip.crit){
if(length(pts)>=2)critv<-smmcrit(test$df,length(pts))
if(length(pts)==1)critv<-qt(.975,test$df)
}
cilow<-test$dif-critv*test$se
cihi<-test$dif+critv*test$se
mat[i,7]<-cilow
mat[i,8]<-cihi
mat[i,9]<-test$p.value
mat[i,10]<-critv
}}
if(plotit){
runmean2g(x1,y1,x2,y2,fr=fr1,est=mean,tr=tr,sm=sm,xout=FALSE,LP=LP,
SCAT=SCAT,xlab=xlab,ylab=ylab,pch1=pch1,pch2=pch2,...)
}
list(output=mat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.