Как выполнить стратифицированный хи-квадрат в R - PullRequest
Новая
       25

Как выполнить стратифицированный хи-квадрат в R

0 голосов
/ 28 июня 2019

Я пытаюсь выполнить тест хи-квадрат (или рыбаки) на нескольких слоях в толще.

Я попытался адаптировать ответ здесь , но меня интересует использование количества экземпляров в таблице, а не вычисление в процентах. Я не слишком знаком с используемыми формулами, поэтому не уверен, как адаптировать то, что там написано. Это было то, что я думал, были корректировки, но это не сработало 

set.seed(1)
df <- tibble(
  Drug = rep(c("Mero", "Vanco"), each=12),
  Service = rep(c("Cardiology", "Urology", "Neurology", "Pediatrics"), each=3, times=2),
  pre_post = rep (c("pre-aso", "post-aso", "post-int"), times = 8),
  LOT = round(runif(24, min=0, max=25), 0)) %>%
  mutate(`LOT >3` = if_else(LOT > 3, 1,  0),
         `LOT <=3` = if_else(LOT < 4, 1,  0)) %>%
  select(Drug, LOT, Service, pre_post, `LOT >3`, `LOT <=3`)

lst <- with(df, split(df, list(Service, pre_post)))
res <- lapply(lst, function(x) chisq.test(x[, -(1:4)]))
sapply(res, "[", "p.value")

В идеале, я бы имел значения хи-квадрат для числа LOT> 3 (против LOT <= 3) по сравнению с pre_post, стратифицированного как Drug, так и Service. Таким образом, многослойная таблица непредвиденных обстоятельств может выглядеть так (числа только заполнители): </p> 

Drug = Mero, Service = Cardiology

+----------+---------+----------+
|          | LOT > 3 | LOT <= 3 |
+----------+---------+----------+
| pre-aso  |       1 |        1 |
| post-aso |       1 |        1 |
| post-int |       1 |        1 |
+----------+---------+----------+

Мой код в настоящее время выдает либо бессмысленные числа, либо NaN. Я также подумал, можно ли использовать dplyr, как показано ниже, но это тоже не сработало. 

t1<- df %>%
  count(Service, Drug, `LOT <=3`) %>%
  table()

EDIT

Мне удалось выяснить один из возможных способов, используя вложенные тиблы, и натолкнулся на это , пытаясь это подтвердить. Однако предложенный способ dplyr дает мне ошибку, поэтому я отказался от нее перед переходом во вложенные списки.

Если я использую приведенный ниже код, я получаю эту ошибку, даже когда я указываю все как фактор и уровни. 

df %>%
  select(Service, Drug, pre_post, `LOT <=3`) %>%
  group_by(Service, Drug) %>%
  summarise(pvalue = chisq.test(pre_post, `LOT <=3`)$p.value)

# Error in summarise_impl(.data, dots) : 
#   Evaluation error: 'x' and 'y' must have at least 2 levels.

Однако это работает. 

df %>% 
  select(Service, Drug, pre_post, `LOT <=3`) %>%
  nest(-Service, -Drug) %>% 
  mutate(freq = map(data, ~table(.))) %>% 
  mutate(pvalue = map(freq, ~chisq.test(.)$p.value)) %>% 
  select(Service, Drug, pvalue) %>% 
  unnest()

Есть мысли, почему я не могу использовать более интуитивный метод dplyr?

...