sbconf: Bootstrap Confidence Intervals
In scaleboot: Approximately Unbiased P-Values via Multiscale Bootstrap

Description Usage Arguments Details Value Author(s) See Also Examples

A confidence interval for a scalar parameter is obtained by inverting the approximately unbiased p-value. This function is very slow, and it is currently experimental.

sbconf(x, ...)

## Default S3 method:
sbconf(x,sa, probs=c(0.05,0.95), model="poly.2",
       k=2,s=1,sp=-1, cluster=NULL,...)

## S3 method for class 'sbconf'
sbconf(x, probs=x$probs,model=x$model,
       k=x$k,s=x$s,sp=x$sp, nofit=FALSE, ...)

## S3 method for class 'sbconf'
plot(x,model=x$model,k=x$k,s=x$s,sp=x$sp,
     models = attr(x$fits,"models"), log.xy = "",
     xlab="test statistic",ylab=NULL, type.plot = c("p","l","b"),
     yval=c("aic","zvalue","pvalue"), sd=2,add=FALSE, col=1:6,
     pch=NULL,lty=1:5,lwd=par("lwd"), mk.col=col[1],
     mk.lwd=lwd[1], mk.lty=lty[1], ...)

`x`	an object used to select a method. For `sbconf.default`, `x` is a list vector of size `length{sa}` with each element being a vector of bootstrap replicates of a statistic or a list vector of a scalar component.
`...`	further arguments passed to or from other methods.
`sa`	vector of scales in sigma squared (σ^2).
`probs`	a vector of probabilities at which p-values are inverted.
`model`	a character to specify a model for an AU p-value. This should be included in `sboptions("models")`, for which model fitting is made internally.
`k`	a numeric to specify an order of AU p-value.
`s`	σ_0^2
`sp`	σ_p^2
`cluster`	parallel cluster object which may be generated by function `makeCluster`.
`nofit`	logical. No further calls to `sbfit` are made.
`models`	AIC values are plotted for these models.
`log.xy`	character string to be passed to `plot.default`.
`xlab`	a label for the x axis.
`ylab`	a label for the y axis.
`type.plot`	a character to be passed to `plot.default`.
`yval`	determines y-axis. "aic" for AIC values of models, "zvalue" for AU corrected z-values, and "pvalue" for AU corrected p-values.
`sd`	If positive, draws curves +-sd*standard error for z-values and p-values.
`add`	logical. Should not the frame be drawn?
`col`	vector of colors of plots.
`pch`	vector of pch's of plots.
`lty`	vector of lty's of plots.
`lwd`	numeric of lwd of plots.
`mk.col`	color for crosses drawn at `probs`.
`mk.lwd`	lwd for crosses drawn at `probs`.
`mk.lty`	lty for crosses drawn at `probs`.

Let x[[i]] be a vector of bootstrap replicates for a statistic with scale sa[i]. For a threshold value y, the bootstrap probability is bp[i]=sum(x[[i]]<y)/length(x[[i]]). sbconf computes bp for several y values, and finds a value y at which the AU p-value, given by sbfit, equals a probability value specified in probs. In this manner, AU p-values are inverted to obtain bootstrap confidence intervals.

See the examples below for details.

sbconf method returns an object of class "sbconf".

The print method for an object of class "sbconf" prints the confidence intervals.

Hidetoshi Shimodaira

scaleboot.

## An example to calculate confidence intervals
## The test statistic is that for "t4" in data(mam15)

data(mam15) # load mam15.mt
sa <- 10^seq(-2,2,length=13) # parameter for multiscale bootstrap

## Definition of a test statistic of interest.
## "myfun" returns the maximum difference of log-likelihood value
## for a tree named a.
myfun <- function(x,w,a) maxdif(wsumrow(x,w))[[a]]
maxdif <- function(x) {
  i1 <- which.max(x)  # the largest element
  x <- -x + x[i1]
  x[i1] <- -min(x[-i1])  # the second largest value
  x
}
wsumrow <- function(x,w) {
  apply(w*x,2,sum)*nrow(x)/sum(w)
}

## Not run: 
## a quick example with nb=1000 (fairely fast in 2017)
## Compute multiscale bootstrap replicates
nb <- 1000 # nb = 10000 is better but slower
# the following line takes some time (less than 1 minute in 2017)
sim <- scaleboot(mam15.mt,nb,sa,myfun,"t4",count=FALSE,onlyboot=TRUE)

## show 90
## each tail is also interpreted as 95
(conf1 <- sbconf(sim$stat,sim$sa,model="sing.3",k=1)) # with k=1
(conf2 <- sbconf(conf1,model="sing.3",k=2)) # with k=2
(conf3 <- sbconf(conf2,model="sing.3",k=3)) # with k=3

## plot diagnostics for computing the confidence limits
plot(conf3) # AIC values for models v.s. test statistic value
plot(conf3,yval="zval",type="l") # corrected "k.3" z-value

## End(Not run)


## Not run: 
## a longer example with nb=10000 (it was slow in 2010)
## In the following, we used 40 cpu's.
nb <- 10000

library(parallel)
cl <- makeCluster(40)
clusterExport(cl,c("maxdif","wsumrow"))

## Compute multiscale bootstrap replicates
## (It took 80 secs using 40 cpu's)
sim <- scaleboot(mam15.mt,nb,sa,myfun,"t4",count=FALSE,
                 cluster=cl,onlyboot=TRUE)

## Modify option "probs0" to a fine grid with 400 points
## default: 0.001 0.010 0.100 0.900 0.990 0.999
## NOTE: This modification is useful only when cl != NULL,
##   in which case calls to sbfit for the grid points
##   are made in parallel, although iterations seen later
##   are made sequentially.
sboptions("probs0",pnorm(seq(qnorm(0.001),qnorm(0.999),length=400)))

## Calculate bootstrap confidence intervals using "k.1" p-value.
## (It took 70 secs using 40 cpu's)
## First, sbfit is applied to bp's determined by option "probs0"
## Then, additional fitting is made only twice for iteration.
## p[1]=0.05 iter=1 t=4.342723 e=0.0003473446 r=0.0301812
## p[2]=0.95 iter=1 t=42.76558 e=0.002572495 r=0.1896809
conf1 <- sbconf(sim$stat,sim$sa,model="sing.3",k=1,cluster=cl)

## The confidence interval with "k.1" is printed as
##     0.05      0.95 
## 4.342723 42.765582 
conf1 

## Calculate bootstrap confidence intervals
##                        using "k.2" and "k.3" p-values.
## (It took only 10 secs)
## p[1]=0.05 iter=1 t=-2.974480 e=0.003729190 r=0.04725755
## p[2]=0.95 iter=1 t=39.51767 e=0.001030929 r=0.06141937
##      0.05      0.95 
## -2.974480 39.517671 
conf2 <- sbconf(conf1,model="sing.3",k=2)
conf2
## p[1]=0.05 iter=1 t=-3.810157 e=0.01068678 r=0.08793868
## p[2]=0.95 iter=1 t=39.32669 e=0.001711107 r=0.09464663
##      0.05      0.95 
## -3.810157 39.326686 
conf3 <- sbconf(conf2,model="sing.3",k=3)
conf3

## plot diagnostics
plot(conf3) # AIC values for models v.s. test statistic value
plot(conf3,yval="zval",type="l") # corrected "k.3" z-value

stopCluster(cl)

## End(Not run)