title: Large fixed effects binomial regression in R
body: <p>I need to run a logistic regression on a relatively large data frame with 480.000 entries with 3 fixed effect variables. Fixed effect var A has 3233 levels, var B has 2326 levels, var C has 811 levels. So all in all I have 6370 fixed effects. The data is cross-sectional. If I can't run this regression using the normal <code>glm</code> function because the regression matrix seems too large for my memory (I get the message "<code>Error: cannot allocate vector of size 22.9 Gb</code>"). I am looking for alternative ways to run this regression on my Macbook Air (OS X 10.9.5 8GB RAM). I also have access to a server with 16GB RAM.</p>
  
  <p>I have tried solving the issue in a few different ways but so far none led to satisfactory results:</p>
  
  <p><strong>lfe/felm</strong>: 
  Using the felm regression function of the <a href="http://cran.r-project.org/web/packages/lfe/lfe.pdf" rel="noreferrer"><code>lfe</code> package</a> that subtracts fixed effects before running the regression. This works perfectly and allowed my to run the above regression as a normal linear model in just a few minutes. However, <code>lfe</code> does not support logistic regressions and glms. So felm was great for getting an idea about model fit for different models but doesn't work for the final logistic regression models.</p>
  
  <p><strong>biglm/bigglm</strong>: 
  I thought about using <a href="http://cran.r-project.org/web/packages/biglm/biglm.pdf" rel="noreferrer"><code>bigglm</code></a> to break my function into more manageable chunks. However, several sources (e.g. <a href="https://stackoverflow.com/questions/10502882/r-biglm-with-categorical-variables">link1</a>, <a href="https://stackoverflow.com/questions/24317418/model-run-with-glm-but-not-bigglm">link2</a>, <a href="http://blog.revolutionanalytics.com/2014/05/quick-history-2-glms-r-and-large-data-sets.html" rel="noreferrer">link3</a>) mention that in order for that to work, factor levels need be consistent across chunks, i.e. each chunk must contain at least one of each factor of each factor variable. Factor A and B contain levels that only appear once, so I can't split the sets into different chunks with consistent levels. If I delete 10 factors of fixed effect A and 8 factors of B (a minor change) I will only have factors with 4+ levels left, and splitting my data into 4 chunks will make it a lot more manageable already. However, then I still need to figure out how to sort my df in a way that will ensure that my 480.000 entries are sorted into 4 chunks in which each factor level of each of the 3 factors appearing at least once.</p>
  
  <p><strong>GlmmGS/glmgs</strong>:
  The <a href="http://cran.r-project.org/web/packages/glmmGS/glmmGS.pdf" rel="noreferrer"><code>glmmgs</code></a> function in the package with the same name performs a fixed-effects subtraction like the <code>lfe</code> package for logistic regressions using a "Gauss-Seidel" Algorithm. Unfortunately, the package is no longer being developed. Being relatively new to R and having no deep experience with statistics I can't make sense of the output and have no idea of how to transform it in a way that would give me the normal "effect size", "model fit", "significance interval" indicators that glm regression summaries provide.</p>
  
  <p>I sent a message to the package's authors. They kindly responded as follows:</p>
  
  <blockquote>
    <p>The package provides no output in the same format of a glm object. However, you
    can easily calculate most of the fit statistics (standard error of the
    estimates, goodness of fit) given the current output (in the CRAN
    version, I believe that the current output is a vector of estimate of
    coefficients, and the associated vector of standard errors; same for
    the covariance components, but you need not worry about them if you
    are fitting model without random effects). Only beware that the
    covariance matrices used to calculated the standard errors are the
    inverse of the diagonal blocks of the precision matrix associated with
    the Gauss-Seidel algorithm, and so they tend to underestimate the
    standard errors of the joint likelihood. I am not maintaining the
    package any longer and I do not have time to get into the specific
    details; the seminal theory behind the package can be found in the
    <a href="http://www.stat.missouri.edu/~guhasu/gauss_seidel.pdf" rel="noreferrer">paper referenced in the manual</a>, everything else needs to be worked out
    by you with pen and paper :).</p>
  </blockquote>
  
  <p>If anyone can explain how to "easily calculate most of the fit statistics" in a way that someone without any education in statistics can understand it (might be impossible) or provide R code that shows on example of how this could be done I would be much obliged!</p>
  
  <p><strong>Revolution Analytics</strong>:
  I installed revolution analytics enterprise on a virtual machine that simulates Windows 7 on my Mac. The program has a function called <code>RxLogit</code> that is optimized for large logistic regressions. Using the <code>RxLogit</code> function I get <code>the error (Failed to allocate 326554568 bytes. Error in rxCall("RxLogit", params) : bad allocation)</code>, so that function also seems too run into memory issues. However, the software enables me to run my regression on a distributed computing cluster. So I could just "kill the problem" by purchasing computing time on a cluster with lots of memory. However, I wonder if the revolution analytics program provides any formulas or methods that I don't know off that would allow me to do some kind of <code>lfe</code>-like fixed-effects subtraction operation or <code>bigglm</code>-like chunking operation that takes factors into account.</p>
  
  <p><strong>MatrixModels/glm4</strong>:
  One person suggested I use the <code>glm4</code> function of the <code>MatrixModels</code> package with the <code>sparse = TRUE</code> attribute to speed up calculation. If I run a <code>glm4</code> regression with all fixed effects I get an <code>"Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed</code>" error. If I run it only with the fixed effect variables B OR A and C, the calculation works and returns a <code>"glpModel"</code> object. As with <code>glmmGS</code> I have some issues in turning that output into a form that makes sense to me since the standard <code>summary()</code> method does not seem to work on it.</p>
  
  <p>I would be happy for advice on any of the issues mentioned above or also completely different approaches for running logistic regressions with multiple large fixed effects in R with memory constraints.</p>