5
$\begingroup$

This question is not meant to be a software question, but I will illustrate the issue using R a bit.

My Understanding of the Simple Case

If I have a simple linear model with a categorical variable on $k$ levels, my understanding of the approach is to replace that categorical variable with $(k-1)$ dummy variables.

So if $A_i$ is the categorical variable for subject $i$ that can take three levels then one way of doing the dummy variables is:

$Y=\beta_1 X_1 + \beta_2 X_2 + c + \epsilon$

Where $X_1$ and $X_2$ are $1$ or $0$ to flag if $A_i$ is in the second or third level (and both are $0$ to indicate the first level).

This means the design matrix is $3 \times N$.

Of course, with just one category, I can just as well split the data up by category and average each set independently, thus avoiding the linear algebra completely. This clearly scales as big as I want without problems, since I can even scan through my data if $N$ is truly gigantic.

The Problem

The problem comes if I have two (or more) categorical variables, each with a large number of possible levels.

Consider the following R code:

library(data.table)
N = 4000000
C1 = 2000
C2 = 500
DT = data.table(A = factor(sample(C1, N, rep = T)), B = factor(sample(C2, N, rep=T)), Y = rnorm(N))
modelFit <- lm(Y~ A+ B, data = DT)

This toy example is a 4m row dataset with two categorical variables, one of 2000 levels and the other of 500. That should result in a total of 2500 parameters, but fitted with 4m datapoints, so this should be sufficient data.

The problem is by the method I described above this needs a $2500 \times 4000000$ design matrix. Even at 1 byte per cell, that matrix would take 9.3gb!

It seems that lm in R hits a problem like this, since when I run the above code I get:

Error: cannot allocate vector of size 74.5 Gb
In addition: Warning messages:
1: In model.matrix.default(mt, mf, contrasts) :
  Reached total allocation of 8087Mb: see help(memory.size)

The Question

So, is there an alternative approach, like in the single category case, that can get the result without this huge design matrix? Some sort of clever averaging procedure?

It seems like such a simple model having only categorical variables.

Improve this question
$\endgroup$
14
  • $\begingroup$ This question appears to be off-topic because it is only about software. Morever, the syntax is incorrect--the arguments to lm have been reversed--and the problem cannot be reproduced (this machine completed the task in 18 seconds using less than 1GB RAM). $\endgroup$
    – whuber ♦
    Commented Nov 13, 2014 at 18:52
  • 1
    $\begingroup$ Do you need all the bells & whistles of, say, lm(), or do you only need point estimates from the fit? Can you load all of the data into your software, or do you need to load & removed parts of your data at a time? $\endgroup$
    – gung - Reinstate Monica
    Commented Nov 13, 2014 at 19:50
  • 1
    $\begingroup$ @whuber ok so the computing solution is a sparse matrix package for the design matrix, the $(X^T X)^{-1} X^T Y$ should calculate ok... Still feel like there might be another statistical way to view it though, but may be not... $\endgroup$
    – Corvus
    Commented Nov 13, 2014 at 21:32
  • 1
    $\begingroup$ Yes, that's precisely what I calculated. As you point out, $X$ itself contains only about ten million nonzero entries and so requires only a few hundred million bytes of storage. The calculation of $X^\prime X$ can occur in-place, creating a matrix with less than two million entries. You don't actually want to find $(X^\prime X)^{-1}$, so you just solve the system $X^\prime X \beta = X^\prime Y$, which is what took most of the two seconds. $\endgroup$
    – whuber ♦
    Commented Nov 13, 2014 at 21:36
  • 1
    $\begingroup$ No, that's the hard way! It takes extensive computational effort to invert a $p\times p$ matrix: it amounts to solving $p$ separate systems of linear equations, whereas you need to solve only one such system. In your case that would be $1300$ times as much work as you really need to do :-). Another nice aspect of this is that you can solve these normal equations even when $X^\prime X$ is not invertible--but in such cases there are entire linear families of solutions (all of which are equally good fits). A good linear equation solver won't be bothered by such complications at all. $\endgroup$
    – whuber ♦
    Commented Nov 13, 2014 at 21:57

1 Answer 1

Reset to default
4
$\begingroup$

I would consider the following two approaches: dimensionality reduction and sparse regression. The traditional methods for the dimensionality reduction are projective and manifold, as mentioned in my answer here and the reference within, as well as latent variables modeling (LVM) methods, such as exploratory factor analysis (EFA) and confirmatory factor analysis (CFA).

Recently I've run across two R packages - ClustOfVar and clere - that essentially implement dimensionality reduction, but differently than traditional methods like PCA, and, in my view, are closer in their nature to the LVM approach. This alternative approach is referred to as variable clustering and you can find more details in the packages' JSS vignettes: this paper and this paper (unpublished yet?), correspondingly.

Additionally, since high-dimensional data is frequently sparse, the second main approach is to use sparse regression. R ecosystem offers many packages, useful for sparse regression analysis, such as Matrix, SprseM, MatrixModels, glmnet and flare. For links and more relevant resources, please see my related answer on DS SE site: https://datascience.stackexchange.com/a/918/2452.

For some overview and more specific examples of support for categorical explanatory variables by PCA, MDS and MCA methods as well as latent variable modeling approaches, please see this paper, these presentation slides (starting from slide 15), this paper, this paper and this paper.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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