R/regpord.R

regpord <-
function(x,y,nboot=100,alpha=.05,SEED=TRUE,xout=FALSE,cov.fun=cov.mba,pr=TRUE,
plotit=TRUE,xlab="Standardized Predictors",ylab="Y",est=mean,scat=var,...){
#
# Compare strength of association of two predictors via
# some robust covariance matrix, with predictors standardized.
#
x<-as.matrix(x)
p1<-ncol(x)+1
p<-ncol(x)
xy<-cbind(x,y)
xy<-elimna(xy)
x<-xy[,1:p]
y<-xy[,p1]
if(xout){
m<-cbind(x,y)
flag<-outfun(x,plotit=FALSE)$keep
m<-m[flag,]
x<-m[,1:p]
y<-m[,p1]
}
n=nrow(x)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
if(pr)print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,regpord.sub,x,y,cov.fun)
ptot=(p^2-p)/2
# bvec is a p by nboot matrix.
est=regvarp(x,y,est=est,scat=scat)
regci<-matrix(0,ptot,4)
dimnames(regci)<-list(NULL,c("Pred.","Pred","test.stat","Decision"))
ic=0
crit05=2.06-5.596/sqrt(n)
if(pr){
print("est is the estimated generalized variance")
}
if(p==2){
if(plotit){
z=standm(x,locfun=lloc,est=mean,scat=var)
z1=cbind(z[,1],y)
z2=cbind(z[,2],y)
plot(rbind(z1,z2),type="n",xlab=xlab,ylab=ylab)
points(z1,pch="*")
points(z2,pch="+")
}}
for(j in 1:p){
for(k in 1:p){
if(j<k){
sqse<-mean((bvec[j,]-est[j]-bvec[k,]+est[k])^2)*nboot/(nboot-1)
test=(est[j]-est[k])/sqrt(sqse)
ic=ic+1
regci[ic,1]<-j
regci[ic,2]<-k
regci[ic,3]<-test
regci[ic,4]<-0
if(abs(test)>=crit05)regci[ic,4]<-1
}}}
regci=data.frame(regci)
flag=(regci[,4]==0)
regci[flag,4]="fail to reject"
regci[!flag,4]="reject"
list(crit.value=crit05,est=est,results=regci)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.