Как отображать эффекты взаимодействия с чрезвычайно большими наборами данных (например, из вывода rxGlm)

Я вычисляю вычисления glm с огромным набором данных. Для вычисления glm и даже speedglm требуется несколько дней.

В настоящее время у меня около 3M наблюдений и всего 400 переменных, только некоторые из них используются для регрессии. В моей регрессии я использую 4 целочисленные независимые переменные (iv1, iv2, iv3, iv4), 1 двоичная независимая переменная как фактор (iv5), член взаимодействия (x * y, где x - целое число, а y - двоичная фиктивная переменная как фактор). Наконец, у меня есть фиксированные эффекты в течение лет ff1 и идентификаторы компании ff2. У меня 15 лет и 3000 компаний. Я ввел фиксированные эффекты, добавив их в качестве факторов. Я отмечаю, что особенно фиксированные эффекты 3000 компаний делают вычисления очень медленными в stats glm, а также speedglm.

Поэтому я решил попробовать Microsoft R rxGlm (RevoScaleR), так как это может адресовать больше потоков и процессорных ядер. Действительно, скорость анализа намного быстрее. Кроме того, я сравнил результаты для суб-выборки с одним из стандартных glm, и они соответствовали.

Я использовал следующую функцию:

mod1 <- rxGlm(formula = dv ~ 
                      iv1 + iv2 + iv3+ 
                      iv4 + iv5 +
                      x * y +
                      ff1  + ff2,
                    family = binomial(link = "probit"), data = dat,
                    dropFirst = TRUE, dropMain = FALSE, covCoef = TRUE, cube = FALSE)

Тем не менее, я столкнулся с проблемой при попытке построить термин взаимодействия с помощью пакета effects. Вызывая следующую функцию, я получаю следующую ошибку:

> plot(effect("x*y", mod1))
Error in terms.default(model) : no terms component nor attribute

Я предполагаю, что проблема заключается в том, что rxGlm не хранит данные, необходимые для построения взаимодействия. Я считаю так, потому что объект rxGlm намного меньше, чем glm oject, поэтому, вероятно, содержит меньше данных (80 МБ против нескольких ГБ).

Теперь я попытался преобразовать объект rxGlm в glm через as.glm(). Тем не менее вызов effects() не дает результата и приводит к следующим сообщениям об ошибках:

Error in dnorm(eta) : 
  Non-numerical argument for mathematical function
In addition: Warning messages:
1: In model.matrix.default(mod, data = list(dv = c(1L, 2L,  :
  variable 'x for y' is absent, its contrast will be ignored

Если теперь сравнить оригинальный glm с "преобразованным glm", я обнаружил, что преобразованный glm содержит намного меньше элементов. Например, он не содержит effects, а для контрастов он указывает только contr.treatment для каждой переменной.

Теперь я ищу в основном способ переноса объекта вывода rxGlm в формате, поэтому я могу использовать, если с помощью функции effect(). Если нет способа сделать это, как мне получить график взаимодействия с использованием функций в пакете RevoScaleR, например, rxLinePlot()? rxLinePlot() также довольно разумно быстро, однако, я еще не нашел способ получить типичные эффекты взаимодействия. Я хочу избежать первого вычисления полной модели glm, а затем построить график, потому что это занимает очень много времени.

Ответы

Ответ 1

Если вы можете получить коэффициенты, вы не можете просто бросить свои собственные? Это не будет проблемой размера набора данных

# ex. data
n = 2000
dat <- data.frame( dv = sample(0:1, size = n, rep = TRUE), 
                   iv1 = sample(1:10, size = n, rep = TRUE),
                   iv2 = sample(1:10, size = n, rep = TRUE),
                   iv3 = sample(1:10, size = n, rep = TRUE),
                   iv4 = sample(0:10, size = n, rep = TRUE),
                   iv5 = as.factor(sample(0:1, size = n, rep = TRUE)),
                   x = sample(1:100, size = n, rep = TRUE),
                   y = as.factor(sample(0:1, size = n, rep = TRUE)),
                   ff1  = as.factor(sample(1:15, size = n, rep = TRUE)),
                   ff2  = as.factor(sample(1:100, size = n, rep = TRUE))
                   )

mod1 <- glm(formula = dv ~ 
                      iv1 + iv2 + iv3+ 
                      iv4 + iv5 +
                      x * y +
                      ff1  + ff2,
                    family = binomial(link = "probit"), data = dat)

# coefficients for x, y and their interaction
x1 <- coef(mod1)['x']
y1 <- coef(mod1)['y1']
xy <- coef(mod1)['x:y1']

x <- 1:100
a <- x1*x
b <- x1*x + y1 + xy*x

plot(a~x, type= 'line', col = 'red', xlim = c(0,max(x)), ylim = range(c(a, b)))
lines(b~x, col = 'blue')
legend('topright', c('y = 0', 'y = 1'), col = c('red', 'blue'))

here is how to make a reproduceable