Как вы запускаете chisq.test на двух разных выходных данных dplyr, а затем суммируете его в таблице?
Вопрос, который у меня есть, связан с тем, который я недавно опубликовал здесь
Jaap был фантастическим, потому что он помог мне создать этот замечательный вывод сводных таблиц подсчетов и частот (в процентах) категориальных переменных.
"Реальные данные", которые я анализирую, поступают из двух разных больниц, в каждой из которых разная частота употребления наркотиков, но не всегда одни и те же лекарства.
Резюме от Яапса func Функция отсюда выглядит следующим образом, а вся data.frame приводится ниже (больницы номер один и два):
id AB1 AB2 AB3 AB4 AB5 AB6 AB7 AB8 AB9 AB10 AB11 AB12 AB13 total perc
1 1st gen Cephalosporin 4 0 0 1 1 0 0 0 0 0 0 0 0 6 1.9
2 3rd gen Cephalosporin 44 7 8 1 3 2 0 0 0 0 0 0 0 65 20.5
3 4th gen Cephalosporin 3 3 0 1 2 1 0 0 0 0 0 0 0 10 3.2
Теперь я хотел бы запустить chisq.test (или же Fisher's если частота ниже 5) всех переменных (имен), найденных в id столбец с использованием общей частоты, найденной в total колонка, сравнивая больницу один против больницы два.
Итак, с точки зрения непрофессионала, я хочу ответить на следующий вопрос: "Цефалоспорины 1-го поколения давали чаще в больнице 1 по сравнению с больницей 2?" и т.п.
Поскольку некоторые переменные идентификаторы могут не совпадать в разных больницах, я ожидаю, что это может вернуть NULL расчет.
В идеале я бы хотел обобщить все эти выводы в таблице с соответствующим значением p, чтобы оно выглядело следующим образом:
id Hospital One Total Frequency Hospital Two Total Frequency p-value
xyz 15 30 0.01
Большое спасибо за Вашу помощь.
Все данные можно найти ниже.
ура
РЕДАКТИРОВАТЬ после поднятых вопросов Хашаа:
Это просто фиктивный вывод (в идеале, что я хотел бы иметь).
id Hospital One Total Frequency Hospital Two Total Frequency p-value
xyz n i x.xx
Как уже упоминалось, значение p должно быть получено из chisq.test или же fisher.test,
Я полагаю, что выходные данные должны быть сгенерированы таким образом, с больницей № 1 hosp1 и больница № 2 называется hosp2
# first take those columns of the dplyr output your interested in
hosp1_sel<-hosp1[,c("id","total")]
hosp2_sel<-hosp2[,c("id","total")]
#then merge the data.frames to one so you can perform analysis on one dataframe
new_df <- merge(hosp1_sel, hosp2_sel, by=0)
#this looks like this
> new_df
Row.names id.x total.x id.y total.y
1 1 1st gen Cephalosporin 6 3rd gen Cephalosporin 19
2 10 Trimethoprim 2 Polypeptide 1
3 11 Ureidopenicillin 46 Rifamycin 1
4 12 Carbapenem 19 Tetracycline 1
5 13 Fluorquinolone 17 Lincosamide 1
6 14 Nitromidazole 12 Quinolone 2
7 15 Antifungal 6 Sulfonamides 2
8 16 Oxazolidinone 2 Nitroimidazole 1
9 17 Rifamycin 1 Polymyxine 1
10 18 Polypeptide 1 Colistin 1
11 2 3rd gen Cephalosporin 65 Carbapenem 37
12 3 4th gen Cephalosporin 10 Fluoroquinolone 24
13 4 Aminoglycoside 31 Glycopeptide 32
14 5 Clindamycin 2 Penicillin 29
15 6 Glycopeptide 55 Ureidopenicillin 36
16 7 Macrolide 3 Lipopeptide 4
17 8 Penicillin 36 Macrolid 2
18 9 Tetracycline 2 Aminoglycoside 9
Вот где я застреваю. На мой взгляд, теперь я должен сделать это data.frame шире, чтобы потом иметь возможность запускать что-то вроде:
chisq.test(hosp1$Ureidopenicillin, hosp2$Ureidopenicillin)
Чтобы определить, чаще ли вводили "Уреидопенициллины" в больнице № 1 по сравнению с больницей № 2 и так далее.
Проблема в том, что это фактически сравнивает "подсчеты", а не "пропорции" из таблицы непредвиденных обстоятельств, хотя...
Есть идеи?
О.
Больница № 1 data.frame:
structure(list(id = structure(1:19, .Label = c("1st gen Cephalosporin",
"3rd gen Cephalosporin", "4th gen Cephalosporin", "Aminoglycoside",
"Clindamycin", "Glycopeptide", "Macrolide", "Penicillin", "Tetracycline",
"Trimethoprim", "Ureidopenicillin", "Carbapenem", "Fluorquinolone",
"Nitromidazole", "Antifungal", "Oxazolidinone", "Rifamycin",
"Polypeptide", "Lipopeptide "), class = "factor"), AB1 = c(4L,
44L, 3L, 1L, 1L, 7L, 1L, 7L, 2L, 1L, 12L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L), AB2 = c(0L, 7L, 3L, 7L, 0L, 16L, 2L, 9L, 0L, 0L,
9L, 1L, 2L, 6L, 0L, 0L, 0L, 0L, 0L), AB3 = c(0L, 8L, 0L, 5L,
1L, 13L, 0L, 5L, 0L, 0L, 12L, 4L, 1L, 2L, 0L, 0L, 0L, 0L, 0L),
AB4 = c(1L, 1L, 1L, 6L, 0L, 5L, 0L, 8L, 0L, 0L, 5L, 3L, 4L,
1L, 1L, 1L, 1L, 0L, 0L), AB5 = c(1L, 3L, 2L, 2L, 0L, 4L,
0L, 1L, 0L, 0L, 2L, 4L, 1L, 1L, 2L, 0L, 0L, 0L, 0L), AB6 = c(0L,
2L, 1L, 3L, 0L, 5L, 0L, 1L, 0L, 0L, 2L, 1L, 1L, 2L, 1L, 0L,
0L, 0L, 0L), AB7 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 3L, 0L,
0L, 2L, 2L, 2L, 0L, 0L, 1L, 0L, 1L, 0L), AB8 = c(0L, 0L,
0L, 3L, 0L, 1L, 0L, 0L, 0L, 0L, 2L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 1L), AB9 = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L,
0L, 2L, 2L, 0L, 0L, 0L, 0L, 0L, 0L), AB10 = c(0L, 0L, 0L,
1L, 0L, 2L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L), AB11 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 2L, 0L, 1L, 0L, 0L, 0L, 0L), AB12 = c(0L, 0L, 0L, 0L,
0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L
), AB13 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L), total = c(6, 65, 10, 31, 2,
55, 3, 36, 2, 2, 46, 19, 17, 12, 6, 2, 1, 1, 1), perc = c(1.9,
20.5, 3.2, 9.8, 0.6, 17.4, 0.9, 11.4, 0.6, 0.6, 14.5, 6,
5.4, 3.8, 1.9, 0.6, 0.3, 0.3, 0.3)), class = "data.frame", .Names = c("id",
"AB1", "AB2", "AB3", "AB4", "AB5", "AB6", "AB7", "AB8", "AB9",
"AB10", "AB11", "AB12", "AB13", "total", "perc"), row.names = c(NA,
-19L))
Больница № 2 data.frame:
structure(list(id = structure(1:18, .Label = c("3rd gen Cephalosporin",
"Carbapenem", "Fluoroquinolone", "Glycopeptide", "Penicillin",
"Ureidopenicillin", "Lipopeptide", "Macrolid", "Aminoglycoside",
"Polypeptide", "Rifamycin", "Tetracycline", "Lincosamide", "Quinolone",
"Sulfonamides", "Nitroimidazole", "Polymyxine", "Colistin"), class = "factor"),
AB1 = c(9L, 3L, 1L, 7L, 16L, 22L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L), AB2 = c(2L, 17L, 5L, 8L, 2L, 9L,
1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), AB3 = c(1L,
9L, 4L, 5L, 3L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 0L), AB4 = c(1L, 3L, 3L, 7L, 4L, 3L, 0L, 0L, 2L, 0L,
0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L), AB5 = c(3L, 1L, 4L, 1L,
4L, 1L, 2L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 0L),
AB6 = c(3L, 2L, 4L, 1L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 1L, 0L), AB7 = c(0L, 2L, 3L, 3L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L), total = c(19,
37, 24, 32, 29, 36, 4, 2, 9, 1, 1, 1, 1, 2, 2, 1, 1, 1),
perc = c(9.4, 18.2, 11.8, 15.8, 14.3, 17.7, 2, 1, 4.4, 0.5,
0.5, 0.5, 0.5, 1, 1, 0.5, 0.5, 0.5)), class = "data.frame", .Names = c("id",
"AB1", "AB2", "AB3", "AB4", "AB5", "AB6", "AB7", "total", "perc"
), row.names = c(NA, -18L))
1 ответ
Слияние по hospital <- left_join(hosp1, hosp2, by = "id") %>% select(id, total.x, total.y) привело к
#id total.x total.y
#1 1st gen Cephalosporin 6 NA
#2 3rd gen Cephalosporin 65 19
#3 4th gen Cephalosporin 10 NA
#4 Aminoglycoside 31 9
#5 Clindamycin 2 NA
#6 Glycopeptide 55 32
#7 Macrolide 3 NA
#8 Penicillin 36 29
#9 Tetracycline 2 1
#10 Trimethoprim 2 NA
#11 Ureidopenicillin 46 36
#12 Carbapenem 19 37
#13 Fluorquinolone 17 NA
#14 Nitromidazole 12 NA
#15 Antifungal 6 NA
#16 Oxazolidinone 2 NA
#17 Rifamycin 1 1
#18 Polypeptide 1 1
#19 Lipopeptide 1 NA
Странно что слишком много NAS производится для hosp2, При ближайшем рассмотрении, есть несоответствия между id переменные. Например, 14-й ряд в hosp1 является Nitromidazole тогда как 16-й ряд hosp2 является Nitroimidazoleи я не уверен, что они указывают на одно и то же лекарство.
Во всяком случае, хотя у меня есть некоторые сомнения по поводу вашего использования chisq.testжелаемый результат может быть получен следующим образом
pval <- function(x, y){
ifelse(!is.na(x) & !is.na(y), chisq.test(c(x, y))$p.value, NA)
}
p <- lapply(1:length(hospital$total.x),
function(i){
pval(hospital$total.x[i],hospital$total.y[i])
}
)
hospital$p_value <- unlist(p)
colnames(hospital) <- c("id", "Hospital One Total Frequency", "Hospital Two Total Frequency", "p-value")
Окончательный вывод выглядит
> hospital
# id Hospital One Total Frequency Hospital Two Total Frequency p-value
#1 1st gen Cephalosporin 6 NA NA
#2 3rd gen Cephalosporin 65 19 5.193805e-07
#3 4th gen Cephalosporin 10 NA NA
#4 Aminoglycoside 31 9 5.042182e-04
#5 Clindamycin 2 NA NA
#6 Glycopeptide 55 32 1.366852e-02
#7 Macrolide 3 NA NA
#8 Penicillin 36 29 3.852612e-01
#9 Tetracycline 2 1 5.637029e-01
#10 Trimethoprim 2 NA NA
#11 Ureidopenicillin 46 36 2.694564e-01
#12 Carbapenem 19 37 1.615693e-02
#13 Fluorquinolone 17 NA NA
#14 Nitromidazole 12 NA NA
#15 Antifungal 6 NA NA
#16 Oxazolidinone 2 NA NA
#17 Rifamycin 1 1 1.000000e+00
#18 Polypeptide 1 1 1.000000e+00
#19 Lipopeptide 1 NA NA