Search code examples
rstatistics-bootstrap

boot.ci() can't compute confidence intervals


I am trying compute confidence intervals for eta-squared for the Kruskal-Wallis test using the boot package in R. But it doesn't seem to work and I am not sure why.

Getting the eta-squared (works!):

Taking small steps, I start with a custom function that return eta-squared value. So this works.

# set up
set.seed(123)
library(tidyverse)
library(PMCMRplus)

# custom function to get eta-squared value
kw_eta_h <- function(data,
                     x,
                     y) {
  # creating a dataframe from entered data
  data <- dplyr::select(
    .data = data,
    x = !!rlang::enquo(x),
    y = !!rlang::enquo(y)
  ) %>%
    dplyr::filter(.data = ., !is.na(x), !is.na(y)) %>%
    tibble::as.tibble(x = .)

  # running the function
  fit <-
    PMCMRplus::kruskalTest(
      formula = y ~ x,
      data = data,
      dist = "KruskalWallis"
    )

  # calculating the eta-squared estimate using the H-statistic
  # ref. http://www.tss.awf.poznan.pl/files/3_Trends_Vol21_2014__no1_20.pdf
  effsize <-
    (fit$statistic[[1]] - fit$parameter[[1]] + 1) /
      (fit$parameter[[3]] - fit$parameter[[1]])

  # return the value of interest: effect size
  return(effsize[[1]])
}

# using the function
kw_eta_h(iris, Species, Sepal.Length)
#> [1] 0.6458329

Getting the eta-squared (doesn't work):

Now I use the custom function that was just used in conjunction with the boot package but it produces identical values for the eta-squared and so no confidence intervals are calculated. What am I doing wrong here?

# function to get confidence intervals
kw_eta_h_ci <- function(data,
                        x,
                        y,
                        nboot = 100,
                        conf.level = 0.95,
                        conf.type = "norm",
                        ...) {
  # creating a dataframe from entered data
  data <- dplyr::select(
    .data = data,
    x = !!rlang::enquo(x),
    y = !!rlang::enquo(y)
  ) %>%
    dplyr::filter(.data = ., !is.na(x), !is.na(y)) %>%
    tibble::as.tibble(x = .)

  # eta-squared value
  eta_sq_H <- kw_eta_h(
    data = data,
    x = x,
    y = y
  )

  # function to obtain 95% CI for for eta-squared
  eta_h_ci <- function(data, x, y, indices) {
    # allows boot to select sample
    d <- data[indices, ]

    # running the function
    fit <-
      kw_eta_h(
        data = data,
        x = x,
        y = y
      )

    # return the value of interest: effect size
    return(fit)
  }

  # save the bootstrapped results to an object
  bootobj <- boot::boot(
    data = data,
    x = x,
    y = y,
    statistic = eta_h_ci,
    R = nboot,
    parallel = "multicore",
    ...
  )

  # get 95% CI from the bootstrapped object
  bootci <- boot::boot.ci(
    boot.out = bootobj,
    conf = conf.level,
    type = conf.type
  )

  # extracting ci part
  if (conf.type == "norm") {
    ci <- bootci$normal
  } else if (conf.type == "basic") {
    ci <- bootci$basic
  } else if (conf.type == "perc") {
    ci <- bootci$perc
  } else if (conf.type == "bca") {
    ci <- bootci$bca
  }

  # preparing a dataframe out of the results
  results_df <-
    tibble::as_data_frame(x = cbind.data.frame(
      "eta_sq_H" = eta_sq_H,
      ci,
      "nboot" = bootci$R
    ))

  # returning the results
  return(results_df)
}

# using the function
kw_eta_h_ci(iris, Species, Sepal.Length)
#> [1] "All values of t are equal to  0.645832897963594 \n Cannot calculate confidence intervals"
#> Error in data.frame(..., check.names = FALSE): arguments imply differing number of rows: 1, 0

Created on 2018-11-16 by the reprex package (v0.2.1)


Solution

  • Inside eta_h_ci you create d as the new sample, but then you call the unsampled data within kw_eta_h. This corrects the behavior on my end.

    eta_h_ci <- function(data, x, y, indices) {
    # allows boot to select sample
    d <- data[indices, ]
    
    # running the function
    fit <-
      kw_eta_h(
        data = d, # d instead of data
        x = x,
        y = y
      )
    
    # return the value of interest: effect size
    return(fit)
    }