Как векторизовать функцию ecdf в R?

У меня есть фрейм данных, который выглядит так:

set.seed(42)
data <- runif(1000)    
utility <- sample(c("abc","bcd","cde","def"),1000,replace=TRUE)
stage <- sample(c("vwx","wxy","xyz"),1000,replace=TRUE)
x <- data.frame(data,utility,stage)
head(x)
   data utility stage
1 0.9148060     def   xyz
2 0.9370754     abc   wxy
3 0.2861395     def   xyz
4 0.8304476     cde   xyz
5 0.6417455     bcd   xyz
6 0.5190959     abc   xyz

и я хочу создать кумулятивные функции распределения для уникальных комбинаций полезности и стадии. В моем реальном приложении я создам около 100 файлов cdf, но эти случайные данные будут иметь 12 (4x3) уникальных комбинаций. Но я буду использовать каждый из этих cdf тысяч раз, поэтому я не хочу вычислять cdf каждый раз на лету. Функция ecdf () работает именно так, как мне хотелось бы, за исключением того, что мне нужно ее векторизовать. Следующий код не работает, но это суть того, что я пытаюсь сделать:

ecdf_multiple <- function(x)
{
    i=0
    utilities <- levels(x$utilities)
    stages <- levels(x$stages)
    for(utility in utilities)
    {
        for(stage in stages)
        {
            i <- i + 1
            y <- ecdf(x[x$utilities == utility & x$stage == stage,1])
            # calculate ecdf for the unique util/stage combo
            z[i] <- list(y,utility,stage)
            # then assign it to a data element (list, data frame, json, whatever) note-this doesn't actually work
        }
    }
    z # return value
}

поэтому после запуска ecdf_multiple и присвоения его переменной я бы как-то сослался на эту переменную, передав значение (для которого мне нужен был cdf), утилиту и сцену.

Есть ли способ векторизовать функцию ecdf (или использовать / построить другую), чтобы я мог выводить результат несколько раз, не требуя создания распределений снова и снова?

------- Добавлено в ответ на отличное предложение @Pascal .-------

Как можно расширить это до более общего случая принятия «n» измерений категорий? Это мой удар, основанный на двухмерном случае Паскаля. Обратите внимание, как я пытался присвоить «y»:

set.seed(42)
data <- runif(1000)    
utility <- sample(c("abc","bcd","cde","def"),1000,replace=TRUE)
stage <- sample(c("vwx","wxy","xyz"),1000,replace=TRUE)
openclose <- sample(c("open","close"),1000,replace=TRUE)
x <- data.frame(data,utility,stage,openclose)
numlabels <- length(names(x))-1
y <- split(x, list(x[,2:(numlabels+1)]))
l <- lapply(y,function(x) ecdf(x[,"data"]))

#execute
utility <- "abc"
stage <- "xyz"
openclose <- "close"
comb <- paste(utility, stage, openclose, sep = ".")
# call the function
l[[comb]](.25)

Во время присвоения «y» выше я получаю следующее сообщение об ошибке:

"Error in sort.list(y) : 'x' must be atomic for 'sort.list'
Have you called 'sort' on a list?"