R/calcBootEst.R

Defines functions calcBCaA calcBCaZ0Boot calcBCaZ0 calcBCaP0 calcBCaQuant calcBCaP calcBCaTbl calcPercP calcPercTbl

Documented in calcBCaA calcBCaP0 calcBCaZ0 calcBCaZ0Boot

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