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