R: перекрестная проверка в наборе данных с факторами

Часто я хочу выполнить кросс-проверку в наборе данных, который содержит некоторые фактор-переменные, и после запуска некоторое время процедура перекрестной проверки не выполняется с ошибкой: factor x has new levels Y.

Например, используя пакет boot:

library(boot)
d <- data.frame(x=c('A', 'A', 'B', 'B', 'C', 'C'), y=c(1, 2, 3, 4, 5, 6))
m <- glm(y ~ x, data=d)
m.cv <- cv.glm(d, m, K=2) # Sometimes succeeds
m.cv <- cv.glm(d, m, K=2)
# Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) : 
#   factor x has new levels B

Обновить. Это пример игрушки. Такая же проблема возникает и с более крупными наборами данных, где имеется несколько вхождений уровня C, но ни один из них не присутствует в обучающем разделе.

Функция createDataPartition из пакета caret делает стратифицированную выборку для переменных результата и правильно предупреждает:

Кроме того, для 'createDataPartition, очень малых размеров классов (< = 3) классы могут не отображаться как в учебных, так и в тестовых данных.

Есть два решения, которые spring:

Сначала создайте подмножество данных, сначала выбрав один случайный образец каждого factor level, начиная с самого редкого класса (по частоте), а затем жадно удовлетворяя следующий редкий класс и так далее. Затем с помощью createDataPartition на остальной части набора данных и слияния результатов для создания нового набора данных поезда, который содержит все levels.
Использование createDataPartitions и выполнение отбраковки.

До сих пор вариант 2 работал у меня из-за размеров данных, но я не могу не думать, что должно быть лучшее решение, чем ручная работа.

В идеале, я бы хотел, чтобы решение, которое просто работает для создания разделов и не срабатывает раньше, если нет способа создать такие разделы.

Есть ли фундаментальная теоретическая причина, почему пакеты этого не предлагают? Предлагают ли они это, и я просто не мог их заметить из-за слепого пятна? Есть ли лучший способ сделать эту стратифицированную выборку?

Пожалуйста, оставьте комментарий, если я задам этот вопрос на stats.stackoverflow.com.

Обновление

Это то, что мое ручное решение (2) выглядит следующим образом:

get.cv.idx <- function(train.data, folds, factor.cols = NA) {

    if (is.na(factor.cols)) {
        all.cols        <- colnames(train.data)
        factor.cols     <- all.cols[laply(llply(train.data[1, ], class), function (x) 'factor' %in% x)]
    }

    n                   <- nrow(train.data)
    test.n              <- floor(1 / folds * n)

    cond.met            <- FALSE
    n.tries             <- 0

    while (!cond.met) {
        n.tries         <- n.tries + 1
        test.idx        <- sample(nrow(train.data), test.n)
        train.idx       <- setdiff(1:nrow(train.data), test.idx)

        cond.met        <- TRUE

        for(factor.col in factor.cols) {
            train.levels <- train.data[ train.idx, factor.col ]
            test.levels  <- train.data[ test.idx , factor.col ]
            if (length(unique(train.levels)) < length(unique(test.levels))) {
                cat('Factor level: ', factor.col, ' violated constraint, retrying.\n')
                cond.met <- FALSE
            }
        }
    }

    cat('Done in ', n.tries, ' trie(s).\n')

    list( train.idx = train.idx
        , test.idx  = test.idx
        )
}

Ответ 1

Все согласны с тем, что это оптимальное решение. Но лично, я бы просто try вызов cv.glm, пока он не будет работать, используя while.

m.cv<- try(cv.glm(d, m, K=2)) #First try
class(m.cv) #Sometimes error, sometimes list
while ( inherits(m.cv, "try-error") ) {
m.cv<- try(cv.glm(d, m, K=2))
}
class(m.cv) #always list

Я пробовал его с 100 000 строк в data.fame, и это занимает всего несколько секунд.

library(boot)
n <-100000
d <- data.frame(x=c(rep('A',n), rep('B', n), 'C', 'C'), y=1:(n*2+2))
m <- glm(y ~ x, data=d)

m.cv<- try(cv.glm(d, m, K=2))
class(m.cv) #Sometimes error, sometimes list
while ( inherits(m.cv, "try-error") ) {
m.cv<- try(cv.glm(d, m, K=2))
}
class(m.cv) #always list

Ответ 2

Когда я вызываю трассировку, я получаю следующее:

> traceback()
9: stop(sprintf(ngettext(length(m), "factor %s has new level %s", 
       "factor %s has new levels %s"), nm, paste(nxl[m], collapse = ", ")), 
       domain = NA)
8: model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels)
7: model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels)
6: predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type == 
       "link", "response", type), terms = terms, na.action = na.action)
5: predict.glm(d.glm, data[j.out, , drop = FALSE], type = "response")
4: predict(d.glm, data[j.out, , drop = FALSE], type = "response")
3: mean((y - yhat)^2)
2: cost(glm.y[j.out], predict(d.glm, data[j.out, , drop = FALSE], 
       type = "response"))
1: cv.glm(d, m, K = 2)

И глядя на функцию cv.glm, вы получаете:

> cv.glm
function (data, glmfit, cost = function(y, yhat) mean((y - yhat)^2), 
    K = n) 
{
    call <- match.call()
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
    n <- nrow(data)
    out <- NULL
    if ((K > n) || (K <= 1)) 
        stop("'K' outside allowable range")
    K.o <- K
    K <- round(K)
    kvals <- unique(round(n/(1L:floor(n/2))))
    temp <- abs(kvals - K)
    if (!any(temp == 0)) 
        K <- kvals[temp == min(temp)][1L]
    if (K != K.o) 
        warning(gettextf("'K' has been set to %f", K), domain = NA)
    f <- ceiling(n/K)
    s <- sample0(rep(1L:K, f), n)
    n.s <- table(s)
    glm.y <- glmfit$y
    cost.0 <- cost(glm.y, fitted(glmfit))
    ms <- max(s)
    CV <- 0
    Call <- glmfit$call
    for (i in seq_len(ms)) {
        j.out <- seq_len(n)[(s == i)]
        j.in <- seq_len(n)[(s != i)]
        Call$data <- data[j.in, , drop = FALSE]
        d.glm <- eval.parent(Call)
        p.alpha <- n.s[i]/n
        cost.i <- cost(glm.y[j.out], predict(d.glm, data[j.out, 
            , drop = FALSE], type = "response"))
        CV <- CV + p.alpha * cost.i
        cost.0 <- cost.0 - p.alpha * cost(glm.y, predict(d.glm, 
            data, type = "response"))
    }
    list(call = call, K = K, delta = as.numeric(c(CV, CV + cost.0)), 
        seed = seed)
}

Кажется, проблема связана с вашим чрезвычайно небольшим размером выборки и категориальным эффектом (со значениями "A", "B" и "C" ). Вы соответствуете glm с 2 эффектами: "B: A" и "C: A". В каждой итерации CV вы загружаете образец выборки и устанавливаете новую модель d.glm. Учитывая размер, загрузочные данные гарантированно будут иметь 1 или более итераций, в которых значение "C" не будет выбрано, поэтому ошибка исходит из вычисления установленных вероятностей из модели начальной загрузки из данных обучения, в которых данные валидации имеют Уровень "С" для х не наблюдается в данных обучения.

Фрэнк Харрелл (часто на stats.stackexchange.com) написал в "Стратегии моделирования регрессии", что нужно одобрять проверку образца образца, когда размер выборки мал и/или некоторые подсчеты клеток невелики при анализе категориальных данных. Сингулярность (как вы видите здесь) является одной из многих причин, почему я считаю, что это правда.

Учитывая небольшой размер выборки здесь, вы должны рассмотреть некоторые альтернативные варианты проверки перекрестной выборки, такие как тест перестановки или параметрический бутстрап. Еще одно важное соображение - именно то, почему вы считаете, что основанный на модели вывод неверен. Когда Туки сказал о бутстрапе, он хотел бы назвать это дробовиком. Это сдует голову от любой проблемы, если вы захотите собрать детали.