Tidyverse решение для группировки и циклического перебора данных, чтобы выполнить dunn.test

Question

У меня есть пример данных ниже.

example.df <- data.frame( 
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))

Сначала я хочу сгруппировать по treated, а затем перебрать все числовые переменные в данных, чтобы выполнить тест Данна (pairw.kw из пакета asbio), используя species какпояснительную переменную, извлеките объект summary dataframe и свяжите столбцы из подсписков yes и no в новый объект dataframe.

Я уже получил частичное решение здесь и здесь с использованием частично "аккуратного" подхода, который работает довольно хорошо.Я просто ищу более элегантное решение Tidyverse, чтобы помочь мне научиться быть лучшим R-пользователем.

Любая помощь приветствуется.

Редактировать: у меня есть вывод из кода в частично "аккуратного" решения.

structure(list(var1.Diff = structure(1:2, .Label = c("-7.05229", 
"-2.25"), class = "factor"), var1.Lower = structure(1:2, .Label = 
c("-13.23198", 
"-8.25114"), class = "factor"), var1.Upper = structure(1:2, .Label = 
c("-0.87259", 
"3.75114"), class = "factor"), var1.Decision = structure(1:2, .Label = 
 c("Reject H0", 
 "FTR H0"), class = "factor"), var1.Adj..P.value = structure(1:2, .Label = 
 c("0.025305", 
 "0.462433"), class = "factor"), resp.Diff = structure(1:2, .Label = 
 c("1.10458", 
 "0"), class = "factor"), resp.Lower = structure(1:2, .Label = c("-5.07512", 
 "-6.00114"), class = "factor"), resp.Upper = structure(1:2, .Label = 
 c("7.28427", 
 "6.00114"), class = "factor"), resp.Decision = structure(c(1L, 
 1L), .Label = "FTR H0", class = "factor"), resp.Adj..P.value = 
 structure(1:2, .Label = c("0.726092", 
 "1"), class = "factor"), effect.Diff = structure(1:2, .Label = 
 c("-1.27451", 
 "-0.5625"), class = "factor"), effect.Lower = structure(1:2, .Label = 
 c("-7.4542", 
 "-6.56364"), class = "factor"), effect.Upper = structure(1:2, .Label = 
 c("4.90518", 
 "5.43864"), class = "factor"), effect.Decision = structure(c(1L, 
 1L), .Label = "FTR H0", class = "factor"), effect.Adj..P.value = 
 structure(1:2, .Label = c("0.686047", 
 "0.85424"), class = "factor")), row.names = c("No", "Yes"), class = 
 "data.frame")

Answer 1

Вот аккуратный подход к одновременному запуску нескольких тестов.

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

Расширить этот подход достаточно просто, чтобы запустить 6 тестов, по одному для каждой комбинации группы лечения и переменной ответа (var1, value1 или resp).Например, мы можем преобразовать тиббл из широкоформатного формата (три столбца ответа) в узкий формат (три ответа сгруппированы по строкам), а затем продолжить работу, как описано выше.

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

^{Создано в 2019-03-04 * представьте пакет (v0.2.1)}

А что, если вы хотите быть гибкими в отношении количества и названия ответов?Затем укажите список ответов:

responses <- c("value1", "var1", "resp")

И используйте аккуратную оценку в выражении gather следующим образом:

example.df %>%
  # From wide format to narrow format, with tidy evaluation
  gather(varname, y, !!!rlang::syms(responses))
  # Rest of computation here.....