Я работаю с ggstatsplot, чтобы получить визуальное представление моего статистического анализа.
У меня есть множество наборов данных, очень похожих по составу. Некоторые работают нормально, а другие нет. data1 - рабочий пример, а data2 не работает.
data1 <- structure(list(
treatment = structure(c(1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L),
.Label = c("negative_ctrl", "positive_ctrl", "treatmentA", "treatmentB", "treatmentC", "treatmentD"), class = "factor"),
value = c(1.74501, 2.04001, 1.89501, 1.84001,
1.89501, 9.75001, 8.50001, 8.80001, 11.50001, 10.25001, 7.90001,
9.25001, 11.45001, 7.75001, 7.75001, 7.55001, 8.70001, 8.20001,
6.95001, 6.60001, 7.40001, 7.15001, 8.25001, 9.20001, 8.95001,
6.45001, 6.05001, 5.40001, 7.95001, 6.80001, 4.65001, 6.40001,
6.40001, 6.70001, 5.40001, 3.20001, 2.70001, 4.30001, 4.10001,
3.60001, 4.00001, 3.00001, 4.70001, 3.10001, 3.50001, 6.45001,
5.45001, 4.90001, 7.25001, 4.55001, 4.70001, 6.25001, 5.65001,
6.00001, 5.10001)),
row.names = c(NA, -55L), class = c("tbl_df", "tbl", "data.frame"))
data2 <- structure(list(
treatment = structure(c(1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L),
.Label = c("negative_ctrl", "positive_ctrl", "treatmentA", "treatmentB", "treatmentC", "treatmentD"), class = "factor"),
value = c(1.00001, 1.00001, 1.00001, 1.00001, 1.00001, 6.77501,
5.68751, 5.99201, 8.24501, 7.01251, 4.79501, 5.99126, 8.26276,
5.35376, 5.38751, 4.60251, 5.38901, 4.85201, 4.44401, 5.20501,
6.20701, 5.77001, 4.05201, 3.65126, 3.02401, 4.68351, 3.90001,
2.56951, 3.70001, 3.61901, 3.96401, 2.93601, 1.53901, 1.40801,
2.05601, 2.08501, 1.89701, 1.79501, 1.50001, 2.09151, 1.53551,
1.57501, 3.88851, 3.09151, 2.75501, 4.40626, 2.42001, 2.60951,
3.83501, 3.37151, 3.70001, 2.92701)),
row.names = c(NA, -52L), class = c("tbl_df", "tbl", "data.frame"))
Я называю самый простой анализ для обоих наборов данных:
library(Rmpfr)
library(ggstatsplot)
ggstatsplot::ggbetweenstats(
data = data1,
x = treatment,
y = value,
messages = FALSE )
ggstatsplot::ggbetweenstats(
data = data2,
x = treatment,
y = value,
messages = FALSE )
Для data1 я получаю это:
для данных2 я получаю:
> Error in stats::optim(par = 1.1 * rep(lambda, 2), fn = function(x) { : non-finite value supplied by optim
Сначала я подумал, что проблема может заключаться в нескольких нулях, которые я передал в отрицательном контроле, но сначала я увеличил их на небольшую величину, а затем на 1, чтобы убедиться, что диапазон значений не является проблемой. Единственное несоответствие, которое я вижу, это то, что у меня есть только 7 вместо 10 измерений для обработки A (уровень 3) в данных 2, но 10 в данных 1 (пришлось удалить несколько NA из-за отказа образца). Однако в обоих случаях отрицательный контроль (уровень 1) имеет только 5 значений, и я не думаю, что в этом типе анализа есть проблема с разными размерами выборки между группами.
В этих случаях рекомендуется попробовать базовые графики, например, выделить ящики:
Итак, сравнивая два набора данных:
boxplot(value ~ treatment, data=data1)
boxplot(value ~ treatment, data=data2)
data2
имеет лечение без вариабельности ("negative_ctrl"
), 0 SD. Я предполагаю, что эта функция выполняет некоторые тесты, требующие изменений. Вам нужно будет прочитать документацию для функции, чтобы увидеть, поднимается ли она, но вы можете получить представления, либо удалив эти процедуры, либо заставив очень небольшое количество вариаций, например
# run without negative_ctrl
ggstatsplot::ggbetweenstats(
data = data2[data2$treatment != "negative_ctrl",],
x = treatment,
y = value,
messages = FALSE )
# add some tiny fake variation to force it through (this is a hack)
data3 <- data2
data3[data3$treatment= = "negative_ctrl",][1,][["value"]] <- 1.0001
ggstatsplot::ggbetweenstats(
data = data3,
x = treatment,
y = value,
messages = FALSE )
Интересно, что если вы запустите урезанный пример только с 5 образцами для каждой из двух групп (взятых из моих наборов данных выше), в случае CV = 0 вы получите t is large; approximation invoked.
, а анализ и график отобразятся просто отлично. Поэтому я думаю, что есть ошибка, которая не сразу распознает необходимость аппроксимации в некоторых особых случаях (например, n (группы)> 2).
И если вам интересно, я отправил сообщение о проблеме на странице ggstatsplot github. Я получил записку от Индраджита Патила. Он смог воспроизвести и проследить проблему: «Это проблема с самим пакетом статистики, который не может выполнить дисперсионный анализ Уэлча. Я ничего не могу здесь сделать».
Спасибо. Это сводило меня с ума. Я пытался все больше и больше изолировать, в чем проблема, но я сосредоточился на значениях/образцах и забыл взглянуть на резюме. И да, достаточно просто добавить немного случайного шума, чтобы обойти это. Я думаю, что на самом деле это ошибка в системе, потому что во многих анализах значения могут быть установлены на нижний или верхний пределы обнаружения, а не удалены, поэтому это может естественным образом привести к группам с CV = 0.