I have a large regression problem with a lot of cases, but relatively few independent variables. One of them is a categorical factor with thousands of levels. Robust regression runs forever. In some cases the large number of dummy variables becomes too sparse to calculate with even "normal" lm.

What would usually make sense is to somehow calculate the average for each level of the factor, then adjust the dependent variable accordingly, and do the regression without the big factor. A colleague of mine could remember there is a two-letter R function that does that automatically, but he cannot remember the two letter combination.

Any help would be greatly appreciated.

  • $\begingroup$ Maybe you're thinking of plm? $\endgroup$ – Charlie Apr 16 '13 at 14:12
  • $\begingroup$ How would plm help? $\endgroup$ – Herbert Apr 16 '13 at 14:32
  • $\begingroup$ It sweeps the categorical fixed effects out of your model. $\endgroup$ – Charlie Apr 16 '13 at 14:56
  • $\begingroup$ That seems like a good idea! $\endgroup$ – Herbert Apr 16 '13 at 16:27
  • $\begingroup$ The package 'lfe' will also handle this problem. It implements a method for projecting out factors from OLS-regressions, the method is described in this forthcoming article: dx.doi.org/10.1016/j.csda.2013.03.024 $\endgroup$ – Simen Gaure Apr 17 '13 at 15:49

I would think lme4 would be highly appropriate for this. Treat your huge categorical factor as a practical random effect. I won't go into the theoretical definitions. Alternatively, use sparse.model.matrix() from Matrix to build the design frame and then pass that into glmnet() from glmnet package. (lme4 naturally builds the sparse design matrix so you don't need to use the sparse.model.matrix() before you go into it.)

If you really want to do the 'average for each level' trick, then be sure to calculate each observation's average excluding itself and include a few extra observations with each factor level at the population mean. Then use this derived variable as a feature in your models instead of the categorical variable. If the factor was the only feature, then this result would be identical to lme4 or glmnet (assuming you solved for how many average observations to add).

There are a few blog entries out there that call the 'average for each level' trick impact coding. Also from my experience, if there is a strong dense feature, you might want to fit a simple model on that feature and impact code the residuals by level of the huge categorical factor instead of the pure response.

As mentioned above, this is more practical advice. Other people will probably come along with some stronger theoretical advice.

| cite | improve this answer | |
  • $\begingroup$ I try to calculate each observation's category average excluding itself and use it in the model. It does not work so well, as out-of-sample prediction is poor (compared to simply using the big factor). Did I get your technique right? $\endgroup$ – Herbert Apr 18 '13 at 9:47
  • $\begingroup$ Did you do the second part and add a handful of fake observations at each level with the mean population response? This is analogous to ridging / L2-Penalties. You might have to tune how many pop-mean dummies to add. Having said that, I also highly recommend lme4. It's formula interface is very simple: y~x.other+(1|x.huge.factor). $\endgroup$ – Shea Parkes Apr 18 '13 at 12:54
  • $\begingroup$ Clarification: I do use sparse matrices currently with MatrixModels / lm.fit.sparse. The estimates I get from lme4 are very similar. However, my eventual goal is to get rid of the big factors to (also) be able to do more effective outlier detection or robust regression. So in such case, would you recommend on the "ridging" strategy? $\endgroup$ – Herbert Apr 18 '13 at 14:17
  • $\begingroup$ Well, in high enough dimensions, every point is an outlier. But in simple terms, if you want to do more complex analysis suited to dense design matrices, then yea, ridging the level averages would probably work out okay. $\endgroup$ – Shea Parkes Apr 18 '13 at 15:25
  • $\begingroup$ I just re-read this and saw reference to lm.fit.sparse(). I think I could safely say that glmnet would almost always be a better option than pure least-squares. Especially when you have such a sparse design matrix. $\endgroup$ – Shea Parkes Apr 21 '13 at 1:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.