Update 2022/03/16
The tidyverse has evolved and so should this solution.
library("asbio")
#> Loading required package: tcltk
library("tidyverse")
# It is good practice to set the seed before simulating random data
# Otherwise we can't compare
set.seed(12345)
example.df <-
tibble(
species = sample(c("primate", "non-primate"), 50, replace = TRUE),
treated = sample(c("Yes", "No"), 50, replace = TRUE),
gender = sample(c("male", "female"), 50, replace = TRUE),
var1 = rnorm(50, 100, 5), resp = rnorm(50, 10, 5), value1 = rnorm(50, 25, 5)
) %>%
# Make sure species is a factor as required by `pair.kw`
mutate(
species = factor(species)
)
example.df %>%
# We want to perform the test separately for each group
group_by(
treated
) %>%
group_modify(
# Perform test and extract summary
~ pairw.kw(.x$var1, .x$species, conf = 0.95)$summary
)
#> # A tibble: 2 × 6
#> # Groups: treated [2]
#> treated Diff Lower Upper Decision `Adj. P-value`
#> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 No -2.67949 -8.12299 2.76402 FTR H0 0.334663
#> 2 Yes 4.07071 -2.98047 11.12188 FTR H0 0.257843
It is straightforward to extend this approach to run 6 tests, one for each combination of a treatment group and a response variable (var1, value1 or resp). For example we can convert the tibble from wide format (three response columns) to narrow format (three responses stacked row-wise) and then proceed pretty much as above.
responses <- c("value1", "var1", "resp")
example.df %>%
pivot_longer(
all_of(responses),
names_to = "variable"
) %>%
group_by(
treated, variable
) %>%
group_modify(
~ pairw.kw(.x$value, .x$species, conf = 0.95)$summary
)
#> # A tibble: 6 × 7
#> # Groups: treated, variable [6]
#> treated variable Diff Lower Upper Decision `Adj. P-value`
#> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 No resp -1.46154 -6.90505 3.98197 FTR H0 0.598725
#> 2 No value1 4.62821 -0.8153 10.07171 FTR H0 0.095632
#> 3 No var1 -2.67949 -8.12299 2.76402 FTR H0 0.334663
#> 4 Yes resp 3.75758 -3.2936 10.80875 FTR H0 0.29627
#> 5 Yes value1 -2.97475 -10.02592 4.07643 FTR H0 0.408311
#> 6 Yes var1 4.07071 -2.98047 11.12188 FTR H0 0.257843
Created on 2022-03-16 by the reprex package (v2.0.1)
Old solution
Here is a tidy approach to running multiple tests simultaneously.
library("asbio")
#> Loading required package: tcltk
library("tidyverse")
# It is good practice to set the seed before simulating random data
# Otherwise can't compare
set.seed(12345)
example.df <- tibble(
species = sample(c("primate", "non-primate"), 50, replace = TRUE),
treated = sample(c("Yes", "No"), 50, replace = TRUE),
gender = sample(c("male", "female"), 50, replace = TRUE),
var1 = rnorm(50, 100, 5), resp=rnorm(50, 10,5), value1 = rnorm (50, 25, 5)) %>%
# Make sure species is a factor as required by `pair.kw`
mutate(species = factor(species))
example.df %>%
# Nest the data for each treatment group
nest(- treated) %>%
# Run the test on each treatment group
mutate(test = map(data, ~ pairw.kw(.$var1, .$species, conf = 0.95))) %>%
# There is no broom::tidy method for objects of class pairw
# but we can extract the summary ourselves
mutate(summary = map(test, pluck, "summary")) %>%
# Cast all the factor columns in the summary table to character, to
# avoid a warning about converting factors to characters.
mutate(summary = map(summary, mutate_if, is.factor, as.character)) %>%
unnest(summary)
#> # A tibble: 2 x 8
#> treated data test Diff Lower Upper Decision `Adj. P-value`
#> <chr> <list> <list> <chr> <chr> <chr> <chr> <chr>
#> 1 No <tibble [2~ <S3: pa~ -1.705~ -7.99~ 4.58~ FTR H0 0.595163
#> 2 Yes <tibble [2~ <S3: pa~ -1.145~ -6.45~ 4.16~ FTR H0 0.672655
It is straightforward to extend this approach to run 6 tests, one for each combination of a treatment group and a response variable (var1, value1 or resp). For example we can convert the tibble from wide format (three response columns) to narrow format (three responses stacked rowwise) and then proceed pretty much as above.
example.df %>%
# From wide format to narrow format
gather(varname, y, value1, var1, resp) %>%
# Nest the data for each treatment group and each variable
nest(- treated, - varname) %>%
# Run 6 tests = (number of treatments) * (number of response variables)
mutate(test = map(data, ~ pairw.kw(.$y, .$species, conf = 0.95))) %>%
# There is no broom::tidy method for objects of class pairw
# but we can extract the summary ourselves
mutate(summary = map(test, pluck, "summary")) %>%
# Cast all the factor columns in the summary table to character, to
# avoid a warning about converting factors to characters.
mutate(summary = map(summary, mutate_if, is.factor, as.character)) %>%
unnest(summary)
#> # A tibble: 6 x 9
#> treated varname data test Diff Lower Upper Decision `Adj. P-value`
#> <chr> <chr> <list> <list> <chr> <chr> <chr> <chr> <chr>
#> 1 No value1 <tibbl~ <S3: ~ 3.127~ -3.1~ 9.41~ FTR H0 0.329969
#> 2 Yes value1 <tibbl~ <S3: ~ -1.33~ -6.65 3.97~ FTR H0 0.622065
#> 3 No var1 <tibbl~ <S3: ~ -1.70~ -7.9~ 4.58~ FTR H0 0.595163
#> 4 Yes var1 <tibbl~ <S3: ~ -1.14~ -6.4~ 4.16~ FTR H0 0.672655
#> 5 No resp <tibbl~ <S3: ~ 1.421~ -4.8~ 7.71~ FTR H0 0.657905
#> 6 Yes resp <tibbl~ <S3: ~ 1.145~ -4.1~ 6.45~ FTR H0 0.672655
Created on 2019-03-04 by the reprex package (v0.2.1)
And what if you want to be flexible about the number and the name of responses? Then specify a list of responses:
responses <- c("value1", "var1", "resp")
And use tidy evaluation in the gather statement as follows:
example.df %>%
# From wide format to narrow format, with tidy evaluation
gather(varname, y, !!!rlang::syms(responses))
# Rest of computation here.....
