R/calcPW.R

Defines functions calcPW

Documented in calcPW

#' Compute pair-wise comparisons estimates and bootstrap statistics
#'
#' Calculates the mean or median sample estimates of all
#' pair-wise differences within a factor. Can be set
#' compute the maximum Bayes factors (MBFs) and indicate
#' which of the comparisons are statistically
#' significant after controlling the false discovery rate, if the number of bootstrap repetition (\code{B}) is set to greater than 0.
#'
#' @param data dataframe. Contains numeric variable named 'resp' for responses and character variable named 'group' containing factor levels.  c
#' @param B integer. Number of bootstrap samples per subgroup. If equal to zero, then the MBF is not calculated.
#' @inheritParams ggbootUV
#' @return A tibble of the comparison groups and estimates of differences in \code{calcStat}, as well as the confidence intervals and MBFs if bootstrapping was performed.
#' @export

calcPW = function( data, calcStat, alpha, B, fdr, method = "bca" ){


  ### PRELIMINARIES

  calc_boot_stat = function( x, w ) calcStat(x[w])

  groupVec = data$group
  respVec = data$resp

  combnVec = calcUniqueCombn( data$group )




  ### SAMPLE DIFFERENCES

  plotTbl = tibble( split = combnVec ) %>%
    mutate( split1 = split ) %>%
    separate( split1, into = c( "y", "x" ) ) %>%
    group_by( x, y, split ) %>%
    summarise( origStat = calcStat( respVec[ groupVec == x ] ) - calcStat( respVec[ groupVec == y ] ) ) %>%
    ungroup()

  # if only sample differences required
  if( is.null( fdr ) ) return( plotTbl )



  ### P-VALUES AND DIRECTIONS OF DIFFERENCES

  if( !is.null( fdr ) ){ # if we need to calculate significant differences

    # bootstrap stats for each group
    origBootStatTbl = ldply( split( respVec, groupVec ), function(x) boot( x, calc_boot_stat, B )$t ) %>%
      rename( group = .id, resp = `1` )

    # bootstrap stats of differences between each combn
    diffBootStatTbl = ldply( combnVec, function(x){
      splitVec = str_split( x, "_" )[[1]]
      level1 = splitVec[1]; level2 = splitVec[2]
      bootVec1 = origBootStatTbl %>% filter( group == level1 ) %>% `[[`("resp")
      bootVec2 = origBootStatTbl %>% filter( group == level2 ) %>% `[[`("resp")
      tibble( split = x, x = level1, y = level2, bootStat = bootVec2 - bootVec1 )
    } ) %>%
      as_tibble()

    # bootstrap stats

    if( method == 'perc' ) bootStatTbl = calcPercTbl( diffBootStatTbl, alpha, 0, "both" )

    if( method == 'bca' )  bootStatTbl = calcBCaTbl( diffBootStatTbl, plotTbl, alpha, 0,
      calcStat, "both" )


    # add
    plotTbl = full_join( plotTbl, bootStatTbl, by = c( "split" ) )


  }

}
MiguelRodo/ggboot documentation built on May 1, 2018, 12:22 a.m.