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