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])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.