sbconf {scaleboot}

Bootstrap Confidence Intervals

Description

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.

Usage

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], ...)

Arguments

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 (\sigma^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

\sigma_0^2

sp

\sigma_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.

Details

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.

Value

sbconf method returns an object of class "sbconf".

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

Author(s)

Hidetoshi Shimodaira

See Also

scaleboot.

Examples

## 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)