R/regciMC.R

regciMC <-
function(x,y,regfun=tsreg,nboot=599,alpha=.05,plotit=FALSE,pr=TRUE,
xlab="Predictor 1",ylab="Predictor 2",xout=FALSE,outfun=outpro,SEED=TRUE,...){
#
#   Compute a .95 confidence interval for each of the parameters of
#   a linear regression equation. The default regression method is
#   Theil-Sen estimator.
#
#   When using the least squares estimator, and when n<250, use
#   lsfitci instead.
#
#   Same as the function regci, only a multi-core processor is used.
#
#   The predictor values are assumed to be in the n by p matrix x.
#   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 estimated of
#   the first predictor, etc.
#
#   plotit=TRUE: If there are two predictors, plot 1-alpha confidence region based
#  on the bootstrap samples.
#
library(parallel)
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]
nrem=length(y)
estit=regfun(x,y,xout=xout,...)$coef
if(xout){
if(pr)print("Default for argument outfun is now outpro")
m<-cbind(x,y)
flag<-outfun(x,plotit=FALSE,...)$keep
m<-m[flag,]
x<-m[,1:p]
y<-m[,p1]
}
flagF=identical(regfun,tsreg)
if(flagF){
if(pr){
if(sum(duplicated(y)>0))print("Duplicate values detected; tshdreg might have more power than tsreg")
}}
x=as.matrix(x)
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=TRUE),nrow=nboot)
data=listm(t(data))
bvec<-mclapply(data,regbootMC,x,y,regfun,mc.preschedule=TRUE,xout=FALSE,...)
bvec=matl(bvec)
# bvec is a p+1 by nboot matrix. The first row
#                     contains the bootstrap intercepts, the second row
#                     contains the bootstrap values for first predictor, etc.
p1<-ncol(x)+1
regci<-matrix(0,p1,5)
vlabs="Intercept"
for(j in 2:p1)vlabs[j]=paste("Slope",j-1)
dimnames(regci)<-list(vlabs,c("ci.low","ci.up","Estimate","S.E.","p-value"))
ilow<-round((alpha/2) * nboot)
ihi<-nboot - ilow
ilow<-ilow+1
se<-NA
sig.level<-NA
for(i in 1:p1){
temp=(sum(bvec[i,]<0)+.5*sum(bvec[i,]==0))/nboot
sig.level[i]<-2*(min(temp,1-temp))
bsort<-sort(bvec[i,])
regci[i,1]<-bsort[ilow]
regci[i,2]<-bsort[ihi]
se[i]<-sqrt(var(bvec[i,]))
}
if(p1==3){
if(plotit){
plot(bvec[2,],bvec[3,],xlab=xlab,ylab=ylab)
}}
regci[,3]=estit
regci[,4]=se
regci[,5]=sig.level
list(regci=regci,n=nrem,n.keep=length(y))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.