R/reg2ciMC.R

reg2ciMC <-
function(x,y,x1,y1,regfun=tsreg,nboot=599,alpha=.05,plotit=TRUE,SEED=TRUE,
xout=FALSE,outfun=outpro,pr=TRUE,...){
#
#   Compute a .95 confidence interval for the difference between the
#   the intercepts and slopes corresponding to two independent groups.
#   The default regression method is Theil-Sen.
#
#   Same as reg2ci, only takes advantage of a multi-core processor
#
#   The predictor values for the first group are
#   assumed to be in the n by p matrix x.
#   The predictors for the second group are in x1
#
#   The default number of bootstrap samples is nboot=599
#
#   regfun can be any R function that returns the coefficients in
#   the vector regfun$coef, the first element of which contains the
#   estimated intercept, the second element contains the estimate of
#   the first predictor, etc.
#
library(parallel)
x<-as.matrix(x)
xx<-cbind(x,y)
xx<-elimna(xx)
x<-xx[,1:ncol(x)]
x<-as.matrix(x)
y<-xx[,ncol(x)+1]
x1<-as.matrix(x1)
xx1<-cbind(x1,y1)
xx1<-elimna(xx1)
x1<-xx1[,1:ncol(x1)]
x1<-as.matrix(x1)
y1<-xx1[,ncol(x1)+1]
x=as.matrix(x)
x1=as.matrix(x1)
if(xout){
flag1=outfun(x,plotit=FALSE,...)$keep
flag2=outfun(x1,plotit=FALSE,...)$keep
x=x[flag1,]
y=y[flag1]
x1=x1[flag2,]
y1=y1[flag2]
}
n=length(y)
n[2]=length(y1)
x<-as.matrix(x)
x1<-as.matrix(x1)
est1=regfun(x,y)$coef
est2=regfun(x1,y1)$coef
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
#
data<-matrix(sample(length(y),size=length(y)*nboot,replace=T),nrow=nboot)
data=listm(t(data))
bvec<-mclapply(data,regbootMC,x,y,regfun,mc.preschedule=TRUE,xout=FALSE)
bvec=matl(bvec)
data<-matrix(sample(length(y1),size=length(y1)*nboot,replace=T),nrow=nboot)
data=listm(t(data))
bvec1<-mclapply(data,regbootMC,x1,y1,regfun,mc.preschedule=TRUE,xout=FALSE)
bvec1=matl(bvec1)
bvec<-bvec-bvec1
p1<-ncol(x)+1
regci<-matrix(0,p1,6)
dimnames(regci)<-list(NULL,
c("Parameter","ci.lower","ci.upper","p.value","Group 1","Group 2"))
ilow<-round((alpha/2)*nboot)+1
ihi<-nboot-(ilow-1)
for(i in 1:p1){
temp<-sum(bvec[i,]<0)/nboot+sum(bvec[i,]==0)/(2*nboot)
regci[i,4]<-2*min(temp,1-temp)
bsort<-sort(bvec[i,])
regci[i,2]<-bsort[ilow]
regci[i,3]<-bsort[ihi]
regci[,1]<-c(0:ncol(x))
}
regci[,5]=est1
regci[,6]=est2
if(ncol(x)==1 && plotit){
plot(c(x,x1),c(y,y1),type="n",xlab="X",ylab="Y")
points(x,y)
points(x1,y1,pch="+")
abline(regfun(x,y)$coef)
abline(regfun(x1,y1)$coef,lty=2)
}
list(n=n,output=regci)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.