title: How to plot interaction effects from extremely large data sets (esp. from rxGlm output)
body: <p>I am currenlty computing <code>glm</code> models off a huge data data set. Both <code>glm</code> and even <code>speedglm</code> take days to compute.</p>
  
  <p>I currently have around 3M observations and altogether 400 variables, only some of which are used for the regression. In my regression I use 4 integer independent variables (<code>iv1</code>, <code>iv2</code>, <code>iv3</code>, <code>iv4</code>), 1 binary independent variable as factor (<code>iv5</code>), the interaction term (<code>x * y</code>, where <code>x</code> is an integer and <code>y</code> is a binary dummy variable as factor). Finally, I have fixed effects along years <code>ff1</code> and company ids <code>ff2</code>. I have 15 years and 3000 conmpanies. I have introduced the fixed effects by adding them as factors. I observe that especially the 3000 company fixed effects make the computation very slow in <code>stats</code> <code>glm</code> and also <code>speedglm</code>.</p>
  
  <p>I therefore decided to try Microsoft R's <code>rxGlm</code> (RevoScaleR), as this can address more threads and processor cores. Indeed, the speed of analysis is a lot faster. Also, I compared the results for a sub-sample to the one of standard <code>glm</code> and they matched.</p>
  
  <p>I used the following function:</p>
  
  <pre><code>mod1 &lt;- 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)
  </code></pre>
  
  <p>However, I am facing a problem when trying to plot the interaction term using the <code>effects</code> package. Upon calling the following function, I am receiving the following error:</p>
  
  <pre><code>&gt; plot(effect("x*y", mod1))
  Error in terms.default(model) : no terms component nor attribute
  </code></pre>
  
  <p>I assume the problem is that <code>rxGlm</code> does not store the data needed to plot the interaction. I believe so because the <code>rxGlm</code> object is a lot smaller than the <code>glm</code> oject, hence likely containing less data (80 MB vs several GB).</p>
  
  <p>I now tried to convert the <code>rxGlm</code> object to <code>glm</code> via <code>as.glm()</code>. Still, the <code>effects()</code> call does not yield a result and results in the following error messages:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>If I now compare an original glm to the "converted glm", I find that the converted glm contains a lot less items. E.g., it does not contain <code>effects</code> and for contrasts it states only <code>contr.treatment</code> for each variable.</p>
  
  <p>I am now looking primarily for a way to transpose the <code>rxGlm</code> output object in a format so I can use if with the <code>effect()</code> function. If there is no way to do so, how can I get an interaction plot using functions within the <code>RevoScaleR</code> package, e.g., <code>rxLinePlot()</code>? <code>rxLinePlot()</code> also plots reasonably quick, however, I have not yet found a way how to get typical interaction effect plots out of it. I want to avoid first calculating the full <code>glm</code> model and then plot because this takes very long.</p>