R/lsfitci.R

lsfitci <-
function(x,y,nboot=599,alpha=.05,SEED=TRUE,xout=FALSE,outfun=out){
#
#   Compute a confidence interval for the slope parameters of
#   a linear regression equation when using the least squares estimator.
#
#   For p=1 predictor,
#   this function uses an adjusted percentile bootstrap method that
#   gives good results when the error term is heteroscedastic.
#   For p>1, a standard percentile bootstrap method is used
#   with FWE (the probability of at least one type I error)
#   controlled via the Bonferroni inequality.
#
#   The predictor values are assumed to be in the n by p matrix x.
#   The default number of bootstrap samples is nboot=599
#
#   SEED=T causes the seed of the random number generator to be set to 2,
#   otherwise the seed is not set.
#
#   Warning: probability coverage has been studied only when alpha=.05
#
x<-as.matrix(x)
p<-ncol(x)
pp<-p+1
temp<-elimna(cbind(x,y)) # Remove any missing values.
x<-temp[,1:p]
y<-temp[,p+1]
if(xout){
m<-cbind(x,y)
flag<-outfun(x,plotit=FALSE)$keep
m<-m[flag,]
x<-m[,1:p]
y<-m[,pp]
}
x<-as.matrix(x)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples; please wait")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,regboot,x,y,lsfit) # A p+1 by n matrix. The first row
#                     contains the bootstrap intercepts, the second row
#                     contains the bootstrap values for first predictor, etc.
if(p==1){
if(alpha != .05){print("Resetting alpha to .05")
print("With p=1, unknown how to adjust confidence interval")
print("when alpha is not equal to .05.")
}
ilow<-15
ihi<-584
if(length(y) < 250){
ilow<-13
ihi<-586
}
if(length(y) < 180){
ilow<-10
ihi<-589
}
if(length(y) < 80){
ilow<-7
ihi<-592
}
if(length(y) < 40){
ilow<-6
ihi<-593
}
ilow<-round((ilow/599)*nboot)
ihi<-round((ihi/599)*nboot)
}
if(p>1){
ilow<-round(alpha*nboot/2)+1
ihi<-nboot-ilow
}
lsfitci<-matrix(0,ncol(x),2)
for(i in 1:ncol(x)){
ip<-i+1
bsort<-sort(bvec[ip,])
lsfitci[i,1]<-bsort[ilow+1]
lsfitci[i,2]<-bsort[ihi]
}
bsort<-sort(bvec[1,])
interceptci<-c(bsort[15],bsort[584])
crit.level<-NA
pmat<-NA
if(p>1){
crit.level<-alpha/p
pmat<-matrix(NA,nrow=p,ncol=2)
dimnames(pmat) <- list(NULL, c("Slope","p-value"))
for(pv in 1:p){
pmat[pv,1]<-pv
pp<-pv+1
pmat[pv,2]<-(sum(bvec[pp,]<0)+.5*sum(bvec[pp,]==0))/nboot
temp3<-1-pmat[pv,2]
pmat[pv,2]<-2*min(pmat[pv,2],temp3)
}}
list(intercept.ci=interceptci,slope.ci=lsfitci,crit.level=crit.level,
p.values=pmat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.