Функция хорошо работает на фиктивных данных, на реальных данных "Ошибка: фактор группировки должен иметь ровно 2 уровня"?

Question

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

Error in wilcox.test.formula(tab[[dependent]] ~ as.factor(tab$group),  : 
      grouping factor must have exactly 2 levels

и предупреждающие сообщения:

In wilcox.test.default(x = c(11.2558701380866, 31.8401548036613,  :   cannot compute exact p-value with ties

Итак, «порог» в моей функции , похоже, неправильно разделяет реальные данные на две группы . Кроме того, подмножество реальных данных неверно . Но я не понимаю, почему ?? Структура фиктивных и реальных таблиц выглядит одинаково:

---

Структура фиктивных и реальных данных:

пустышки:

> str(tab)
'data.frame':   80 obs. of  3 variables:
 $ infGrad    : num  14.15 12.53 3.03 9.21 16.36 ...
 $ distance   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ uniqueGroup: Factor w/ 2 levels "x","y": 1 2 1 2 1 2 1 2 1 2 ...

Реальный:

> str(tab)
'data.frame':   142 obs. of  10 variables:
 $ distance   : num  100 100 100 100 100 100 100 100 100 100 ...
 $ infGrad    : num  11.3 17.4 31.8 11.1 47.8 ...
 $ uniqueGroup: Factor w/ 6 levels "x",..: 5 2 5 2 5 5 5 5 3 6 ...

---

Я обнаружил, что NAs может вызвать эти проблемы или спецификацию формулы wilcox.test(y ~ x).

Итак, я попытался добавить na.omit к своей функции и вместо wilcox.test(y~x) использовать wilcox.test(y, x). Ни один из них не сработал.

Есть ли у вас идеи как заставить мою функцию работать или как сделать ее более надежной для приема моих реальных данных? Ваша помощь высоко ценится.

---

Что делает код:

1. классифицирует данные по двум группам по "скользящему порогу"
2. проверить статистические различия между этими группами.

Я запускаю функцию с вложенным lapply, чтобы варьировать мои пороги и разные подмножества данных.

Мои фиктивные данные:

set.seed(10)

infGrad <- c(rnorm(20, mean=14, sd=8),
            rnorm(20, mean=13, sd=5),
            rnorm(20, mean=8, sd=2),
            rnorm(20, mean=7, sd=1))
distance <- rep(c(1:4), each = 20)
uniqueGroup <- rep(c("x", "y"), 40)

tab<-data.frame(infGrad, distance, uniqueGroup)


# Create moving threshols function
movThreshold <- function(th, tab, dependent, ...) {
  tab<-na.omit(tab)

  # Classify data 
  tab$group<- ifelse(tab$distance < th, "a", "b") # does not WORK on REAL data

  # Calculate wincoxon test 
  test<-wilcox.test(tab[[dependent]] ~ as.factor(tab$group),  # specify column name
                    data = tab)

  # Put results in a vector 
  c(th, dependent, round(test$p.value, 3))

}

# Define two vectors to run through
# unique group
gr.list<-unique(tab$uniqueGroup)

# unique threshold
th.list<-c(2,3,4)


# apply function over threshols and subset
res<-lapply(gr.list, function(x) lapply(th.list,
                                        movThreshold,
                                        tab = tab[uniqueGroup == x,], # does not work on REAL data
                                        dependent = "infGrad"))

То, что не работает с реальными данными:

1. Классификация групп внутри функции

    tab$group<- ifelse(tab$distance < th, "a", "b")
   

2. Подмножество данных во вложенном цикле lapply

    subsetting: tab = tab[uniqueGroup == x,]
   

Answer 1

Проблема, вероятно, возникает из-за одной группы значений. Вы можете воспроизвести ошибку, например, добавив высокое значение к th.list.

# unique threshold
th.list<-c(2,3,4,100)

Самый простой способ избежать этого - проверить длину tab$group перед выполнением теста. Это изменение в функции должно быть достаточно:

movThreshold <- function(th, tab, dependent, ...) {
  tab<-na.omit(tab)

  # Classify data 
  tab$group<- ifelse(tab$distance < th, "a", "b") # does not WORK on REAL data

  # Check there are two groups
  if(length(unique(tab$group))<2){return(NA)}
  # Calculate wincoxon test 
  test<-wilcox.test(tab[[dependent]] ~ as.factor(tab$group),  # specify column name
                    data = tab)

  # Put results in a vector 
  c(th, dependent, round(test$p.value, 3))

}