# Multiple linear regression where explanatory variables almost always equal 0

I'm trying to do a multiple linear regression model to correlate a dependent $Y$ variable (normally distributed) against a set of 642 variables. These 642 variables codify for the presence or absence of a particular chemical group in a molecule, therefore they can assume only 1 or 0 as value. The size of the dataset is 1300 experiments (one $Y$) and 642 $X$ variables. The problem I have is that the 642 variables are set to 1 only 0.3% of the times, i.e., my 1300 x 642 X-matrix is made of 2600 1 and 832000 0. Nothing else (except the $Y$ of course). Does high frequency of 0 create problem for the MLR implementation in R? Is there a better statistical method more suitable to treat this case?

• I think you should be using NB or SVM instead. – sds Aug 14 '13 at 22:26

Provided each chemical group appears in at least two experiments you can be successful.

A good way to explore the issues is with a simulation. What matters is how much any random variation in the values of $Y$ compares to the sizes of the effects: when that variation is relatively small, you can pretty reliably estimate the effects and know they are real; when that variation is relatively large, forget about it.

### Let's explore using simulated data.

set.seed(17)
x.0 <- replicate(642, sample(c(rep(0, 1298), rep(1, 2))))
x.1 <- t(replicate(1300, sample(c(rep(0, 641), 1))))
x <- x.0 + x.1
sum(x) # Need to be around 2600
par(mfrow=c(1,2))
plot(apply(x, 1, sum), xlab="Case", ylab="Groups per case", cex=0.6)
plot(apply(x, 2, sum), xlab="Group", ylab="Cases per group", cex=0.6)


The plots confirm that each case includes at least one chemical group and each chemical group appears in at least two cases (so some assessment of its significance becomes possible). There are $1300$ cases, $642$ variables, and $2584$ ones in the dataset, close enough to the description in the problem.

The function sim generates the true values of $Y$ and adds some iid Normal variation of a specified standard deviation sd. It runs the multiple linear regression lm then it plots (a) a histogram of the "observed" values of $Y$ and (b) a histogram of the p-values for p-values between 0 and p.limit. (Its optional parameters are passed on to the two hist calls.) Reducing p.limit lets us zoom into the range of smaller p-values to see them better.

sim <- function(beta, sd, p.limit=1, ...) {
d <- data.frame(x, y=x %*% beta + rnorm(dim(x)[1], sd=sd))
fit <- lm(y ~ ., data=d)
hist(d$y, xlab="Y", main="Histogram of Y", ...) p <- summary(fit)$coefficients[, 4]
hist(p[p <= p.limit], xlab="p", main="Histogram of p-values", ...)
}


### No effects

It is now easy to explore scenarios. Begin with the case where there are no effects:

beta <- rep(0, dim(x)[2])
sim(beta, 1, freq=FALSE)


As we would expect, the p-values are approximately uniformly distributed throughout their possible range. This means you should expect around $.05 \times 642 = 32$ of them to be less than $5$% and around $.01 \times 642 = 6$ to be less than $1$% even though none of the variables has any relationship with $Y$ at all: you have a problem with false positives. One way to cope with that is a Bonferroni correction; it stipulates that only p-values less than $0.05/642$ should be considered significant at the overall $5$% level. (There are better methods, but this is an easy way to illustrate the role played by any multiple-comparisons correction.)

The effects in the next simulations are normally distributed with a standard deviation of $1$, centered around $0$.

### Strong effects

A standard deviation in the errors of only $0.1$ is quite small compared to these effects. Let's see how well ordinary least squares with Bonferroni correction works:

set.seed(17)
beta <- rnorm(dim(x)[2])
sim(beta, .1, p.limit=0.05/642)


It works marvelously well! Most of the effects are so large (compared to the random variation) that virtually all of them are considered significant.

### Weak effects

But don't assume your data will be this nice. Let's look at perhaps a more realistic situation where the effect sizes are comparable to the random variation:

sim(beta, 1, p.limit=0.05/642)


Sadly, only a couple dozen of the effects can be considered significant here. You just don't have the power to detect most of them.

### Conclusions

1. You must protect yourself against false positives by adopting some kind of multiple comparisons procedure.

2. If you have a great dataset with strong effects, ordinary least squares multiple regression can work beautifully. Otherwise, you may discover the data simply are not sufficient to detect the effects that might be there.

It will not create problems, although inclusion of all 642 variables in your regression will. You have only 1300 observations for each variable and if you include all variables you will be asking the data to answer 642 questions (asking whether this variable is significant or not) which is too many in my opinion. Something in the order of 10 variables would be more appropriate for data of this size.

Also important to note is that all your explanatory variables are dummy variables indicating the inclusion of a chemical group or not. This does not cause any problems but you should note that the coefficients estimated in your regression will give you the average Y for that particular variable. For example, if I run a regression of height on a constant and a dummy variable indicating whether someone is male or not (1=male, 0=female), the coefficient of the constant term will only tell me the average height for all females and the sum of the coefficient on the constant term and gender dummy will tell me the average height for all males in the sample. So if your Y is say cost of experiment and one of your dummies indicates the presence of benzine, your constant + coefficient will tell you the average cost of an experiment with benzine in it.

I would be mostly worried about performance of your model eroding very quickly as you add more and more explanatory variables. The ratio of zeros to ones is not a problem though.

• Welcome to our site! The assessment that more than 10 or so variables would not be appropriate needs more nuance and better justification. And could you be more precise about what you mean by "performance" in the last paragraph? – whuber Aug 14 '13 at 20:49

@JoeDanger has provided a good answer (+1); you should have some concern about the possibility of overfitting. A typical rule of thumb in multiple regression is that you want at least 10 observations per variable (see here: rules-of-thumb-for-minimum-sample-size-for-multiple-regression for more discussion of this topic). With 642 variables, that means you want 6,420 observations, whereas you have only 1,300 (~20%). This is only a rule of thumb and there is nothing magical about that number, but you should be cognizant of this.

In addition to his points, there isn't anything in the algorithm that R uses to fit the model that will be thrown off by having mostly 0s, if that's your main concern. However, you will have low power to detect the effects of those variables on the response. You can think of a regression analysis with a dummy variable as being a t-test, and thus you have a situation where your groups are of unequal size. I discussed this issue here: how-should-one-interpret-the-comparison-of-means-from-different-sample-sizes, which may help you to understand the issue more fully. A multiple regression like your situation will be analogous.

• In order to avoid the overfitting problem, would this be a suitable case for the lasso or elastic net (eg, using the glmnet package in R)? – DavidP Aug 15 '13 at 17:25
• I know Max Kuhn (author of "caret") at Pfizer has done work along these lines – mimenico Aug 19 '13 at 8:11