如何从超大的数据集绘制交互效应图(特别是从rxGlm输出)

Question

我目前正在计算一个巨大的数据集上的glm模型。无论是glm还是speedglm都需要几天的时间来计算。

我目前有大约3M个观察值和总共400个变量，其中只有一些用于回归。在我的回归中，我使用了4个整数自变量(iv1，iv2，iv3，iv4)，1个二进制自变量作为因子(iv5)，交互作用项(x * y，其中x是整数，y是作为因子的二进制虚拟变量)。最后，我对多年的ff1和公司ids ff2有固定的影响。我有15年的经验和3000个公司。我已经通过添加它们作为因子引入了固定效果。我观察到，特别是3000公司的固定效果使得stats，glm和speedglm的计算非常慢。

因此，我决定尝试微软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对象小得多，因此可能包含的数据更少(80MB对几GB)。

现在，我尝试通过as.glm()将rxGlm对象转换为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()函数中使用if。如果无法做到这一点，我如何使用RevoScaleR包中的函数(例如rxLinePlot() )获得交互作用图？rxLinePlot()也可以相当快速地绘制，但是，我还没有找到一种方法来获得典型的交互作用图。我希望避免先计算完整的glm模型，然后绘制，因为这需要很长时间。

Answer 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'))

​

​