R/calcBootEst.R In MiguelRodo/ggboot: Bootstrap univariate and multivariate plots

Documented in calcBCaAcalcBCaP0calcBCaZ0calcBCaZ0Boot

```# calculate confidence intervals and p-values using percentile method
calcPercTbl = function( bootTbl, alpha, nullValue, altSide ){
bootTbl %>%
group_by( split ) %>%
summarise( p = calcPercP( bootStat, nullValue, altSide ),
lb =  partial( quantile, p = alpha/2, na.rm = TRUE)( bootStat ),
ub =  partial( quantile, p = 1 - alpha/2, na.rm = TRUE)( bootStat ) )
}

calcPercP = function( bootVec, nullValue, altSide ){
p = sum( bootVec <= nullValue ) / length( bootVec )
if( altSide == "high" ) return( p ) # if the proportion less than the null is the pValue
if( altSide == "low" ) return( 1 - p )
if( altSide == "both" ) return( min( p, 1-p ) )
}

# calculate intervals and p-values using bca method
calcBCaTbl = function( bootTbl, origStatTbl, alpha, nullValue, calcStat, altSide ){

bootTbl %>%
group_by( split ) %>%
summarise( p = calcBCaP( bootStat, nullValue,
origStatTbl\$origStat[ which( origStatTbl\$split == split[1] ) ],
calcStat, altSide ),
lb = calcBCaQuant( bootStat, alpha/2,
origStatTbl\$origStat[ which( origStatTbl\$split == split[1] ) ],
calcStat ),
ub = calcBCaQuant( bootStat, 1 - alpha/2,
origStatTbl\$origStat[ which( origStatTbl\$split == split[1] ) ],
calcStat ) )
}

# calculate the pval using bca
calcBCaP = function( bootVec, nullValue, origStat, calcStat, altSide ){

### PRELIMINARIES
calcGHat = function(x) sum( ( bootVec <= x ) * 1 ) / length( bootVec ) # ghat function

# extreme cases
if(  calcGHat( nullValue ) == 0 ){
if( altSide %in% c( "high", "both" ) ) return( 0 )
if( altSide == "low" ) return( 1 )

}
if( calcGHat( nullValue ) == 1 ){
if( altSide == "high" ) return( 1 )
if( altSide %in% c( "both", "low" ) ) return( 0 )
}

R = qnorm( calcGHat( nullValue ) ) # quantile such that the difference from the median value is less than it
z0 = calcBCaZ0Boot( bootVec, origStat )
a = calcBCaA( bootVec, calcStat )
num = (R - z0) * (1 - a * z0 ) - z0
den = 1 + a * ( R - z0 )
p = pnorm( num/den )

if( altSide == "high" ) return( p ) # if the proportion less than the null is the pValue
if( altSide == "low" ) return( 1 - p )
if( altSide == "both" ) return( min( p, 1-p ) )

}

# calculate bca quantile
calcBCaQuant = function( bootVec, alpha, origStat, calcStat ){

# preliminaries
z0 = calcBCaZ0Boot( bootVec, origStat )
a = calcBCaA( bootVec, calcStat )
za = qnorm( alpha )

# second term
num = z0 + za
den = 1 - a * ( z0 + za )
term2 = num/den

# output the quantile
setNames( quantile( bootVec, pnorm( z0 + term2 ) ), NULL )
}

#' Calculate p0.
#'
#' \code{calcBCaP0} calculates p0 for BCa bootstrap
#' algorithm.
#'
#' @param bootVec numeric vector. Vector of bootstrap estimates from the estimator
#' of interest.
#' @param origStat numeric. Estimate from estimator applied to original data.
#' @return p0
#' @examples
#' calcBCaP0( rnorm(15), 0.1 )
calcBCaP0 = function( bootVec, origStat ) sum( ( bootVec <= origStat ) * 1 ) / length( bootVec )

#' calculate z0
#'
#' \code{calcBCaZ0} calculates z0 for BCa bootstrap algorithm.
#'
#' @param p0 numeric. Proportion of bootstrap estimates less than the
#' estimate from the estimator applied to the whole sample.
#' @return z0
#' @examples  calcBCaZ0( 0.2 )
calcBCaZ0 = function( p0 ) qnorm( p0, lower.tail = TRUE )

#' Calculate z0 from bootstrap data.
#'
#' Wrapper function around calc_bca_p0 and calc_bca_z0 that
#' calculates z0 from raw data. This is different to
#' calc_bca_z0, which calculates z0 based on a given p0.
#'
#' @param bootVec numeric vector. Vector of bootstrap estimates.
#' @param origStat numeric. Estimate from estimate applied to the
#' whole sample.
#' @return z0
#' @examples calcBCaZ0Boot( rnorm( 20 ), 0.05 )
calcBCaZ0Boot = function( bootVec, origStat ) calcBCaZ0( calcBCaP0( bootVec, origStat ) )

#' Calculate acceleration factor, \eqn{a}, for BCa bootstrap algorithm.
#'
#' \code{calcBCaA} calculates the acceleration factor, \eqn{a}, for the
#' BCa algorithm, when the estimator is for a single statistic.
#'
#' @param bootVec numeric vector. Vector of bootstrap estimates.
#' @param calcStat function. The estimator of interest as a function. Returns
#' the estimate.
#' @return \eqn{a}
#' @examples
#' x = rnorm(20)
#' calcBCaA( rnorm( 20 ), mean )

calcBCaA = function( bootVec, calcStat ){
jackVec = llply( seq_along( bootVec ), function( i ) calcStat( bootVec[ -i ]  ) ) %>% unlist()
if( calcLU( jackVec ) < length( jackVec ) /2 ) return( 0 ) # return 0 if there are few unique entries
num = sum( ( jackVec - mean(jackVec) )^3 )
den = sum( ( jackVec - mean(jackVec) )^2 )^1.5
num / den / 6
}
```
MiguelRodo/ggboot documentation built on May 1, 2018, 12:22 a.m.