R/anclin.R

anclin <-
function(x1,y1,x2,y2,regfun=tsreg,pts=NULL,ALL=FALSE,npts=25,plotit=TRUE,SCAT=TRUE,
pch1='*',pch2='+',
nboot=100,ADJ=TRUE,xout=FALSE,outfun=out,SEED=TRUE,p.crit=.015,
alpha=.05,crit=NULL,null.value=0,plotPV=FALSE,scale=TRUE,span=.75,xlab='X',xlab1='X1',xlab2='X2',ylab='p-values',ylab2='Y',theta=50,phi=25,MC=FALSE,nreps=1000,pch='*',...){
#
#  ANCOVA:
#  For two independent groups, compute confidence intervals for difference between
#  the typical value of Y, given X,
#   based on some regression estimator
#  By default, 
#  regfun=tsreg meaning that the  Theil--Sen estimator is used.
#
#  The functions anclin and regYci2g are identical. 
#
#  In contrast to the function ancJN, this function can deal with a larger number of 
#  covariate values and it controls the probability of one or more Type I errors using
#  a method that is better, in terms of power, than using Hochberg or Hommel.
# 
#  ADJ=TRUE,  the critical value is adjusted so that the simultaneous
#   probability coverage is 1-alpha.
#   If there is a single covariate,
#  regfun=tsreg or tshdreg, and alpha=.05, an adjustment can be made quickly. Otherwise an
#  adjustment must be computed, which can require relatively high execution time.
#  To reduce execution time, set 
#  MC=TRUE, assuming a multi-core processor is available.
#  If n1<20 and n2<100, assuming that n1<n2,
#  an adjusted critical value must be computed even when using the
#  Theil--Sen estimator.
#   If the R package WRScpp has been installed, execution time can be reduced further
#   using when using Theil--Sen by setting regfun=tsreg_C or tshdreg_C
#
#  nboot: the number of bootstrap samples used to estimate standard errors.
#
# nreps: Number of replications used to compute a critical value. 
#
#  pts: values for the independent variable where confidence intervals are computed
#  pts=NULL means that 100 points evenly spaced between min(x1,x2) and max(x1,x2) are used.
#  If pts is specified,the function and ADJ=TRUE, the function will  compute an
#  adjusted critical value, which again can result in high execution time.
#
xy=elimna(cbind(x1,y1))
x1<-as.matrix(x1)
p=ncol(x1)
if(p>1)stop('Current version allows one covariate only')
p1=p+1
vals=NA
x1<-xy[,1:p]
y1<-xy[,p1]
x1<-as.matrix(x1)
xy=elimna(cbind(x2,y2))
x2<-as.matrix(x2)
p=ncol(x2)
p1=p+1
vals=NA
x2<-xy[,1:p]
y2<-xy[,p1]
x2<-as.matrix(x2)
n1=length(y1)
n2=length(y2)
n=min(c(n1,n2))  
if(xout){
m<-cbind(x1,y1)
flag<-outfun(x1,plotit=FALSE,...)$keep
m<-m[flag,]
n1=nrow(m)
x1<-m[,1:p]
y1<-m[,p1]
x1=as.matrix(x1)
m<-cbind(x2,y2)
flag<-outfun(x2,plotit=FALSE,...)$keep               
m<-m[flag,]                                               
n2=nrow(m)
n=min(c(n1,n2))                             
x2<-m[,1:p]                                                              
y2<-m[,p1]                                                                     
x2=as.matrix(x2)  
}
if(is.null(pts)){
xall=unique(c(x1,x2))
if(ALL)pts=xall
if(!ALL)pts=seq(min(xall),max(xall),length.out=npts)
}
if(ADJ){
if(n<10)stop('Should have a sample size of at least 10')
if(alpha==.05){
#if(identical(regfun,tsreg) || identical(regfun,tsreg_C))alpha=p.crit causes an error if WRScpp not installed 
#if(identical(regfun,tsreg))alpha=p.crit
alpha=p.crit  
crit=qnorm(1-alpha/2)
}
if(!ADJ)p.crit=alpha
if(n<20 & max(c(n1,n2))<100) crit=NULL
if(p>1)crit=NULL
}
if(is.null(crit) & !ADJ)crit=qnorm(1-alpha/2)
if(is.null(crit) & ADJ){
if(SEED)set.seed(2)
padj=regYciCV2G(n1,n2,nboot=nreps,regfun=regfun,MC=MC,SEED=FALSE,ALL=ALL,
null.value=null.value,pts=pts,alpha=alpha,...)$crit.est 
crit=qnorm(1-padj/2)
p.crit=padj
}
sqsd1=regYvar(x1,y1,regfun=regfun,pts=pts,nboot=nboot,SEED=SEED)
sqsd2=regYvar(x2,y2,regfun=regfun,pts=pts,nboot=nboot,SEED=SEED)
sd=sqrt(sqsd1+sqsd2)
est1=regYhat(x1,y1,regfun=regfun,xr=pts,...)
est2=regYhat(x2,y2,regfun=regfun,xr=pts,...)
pv=2*(1-pnorm(abs(est1-est2-null.value)/sd))
est=cbind(pts,est1-est2,est1-est2-crit*sd,est1-est2+crit*sd,pv)
dimnames(est)=list(NULL,c('X','Est.Dif','Lower.ci','Upper.ci','p.value'))
if(plotit){
plotPV=FALSE
plot(c(x1,x2),c(y1,y2),type="n",xlab=xlab,ylab=ylab2)
reg1=regfun(x1,y1,...)$coef
reg2=regfun(x2,y2,...)$coef
if(SCAT){
points(x1,y1,pch=pch1)
points(x2,y2,pch=pch2)
}
abline(reg1)
abline(reg2,lty=2)
}
if(plotPV){
if(ncol(x1)>2)stop('Can plot only with one or two independent variables')
if(ncol(x1)==1)plot(pts,pv,xlab=xlab,ylab=ylab,pch=pch)
if(ncol(x2)==2)lplot(pts,pv,xlab=xlab1,ylab=xlab2,zlab=ylab,span=span,ticktype='detail',scale=scale,theta=theta,phi=phi)
}
list(output=est,p.crit=p.crit,crit.value=crit,num.sig=sum(est[,5]<=p.crit))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.