# Multivariable regression with binomial covariates and interaction terms. Why is the predictive power of linear regression w/o interactions so high?

TL;DR: How could the $$R^2$$ of a linear model approach 1 when the outcome is solely determined by interaction terms (i.e., there's no underlying "main effect")? It doesn't seem to happen when covariates are normally distributed, only with binomials. Looking for insights.

Hopefully the motivation here is simple enough: imagine we have a number of covariates, plus some interaction terms, and things are set up such that the value of our outcome variable is solely determined by the interaction terms. For some strange reason, we're interested in seeing what a linear model with no interactions comes out to anyway. In other words, we could say we have a data generating process roughly like so:

y = b0 + b1*x1 + b2*x2 + b3*x1*x2


where b1 and b2 are set to 0.

The question is what happens we regress y on x1 and x2 anyway, without modelling the interaction. My intuition was that the $$R^2$$ of such a model would approximate to 0. Indeed, this is what I observe when my covariates are drawn at random from a standard normal.

Something interesting happens when I draw the values of the covariates from a binomial distribution. Depending on the binomial probability and the total number of interactions, the $$R^2$$ of the linear model can approach 1 according to my simulations. See this chart of the $$R^2$$ vs number of interactions and the binomial probability (this shows just the effect of binomial probability)).

I can share my full simulation code if anyone is interested, but to keep things simple if you want to convince yourself of my basic point you can run this R code:

sim_model <- function(dataset, n_interactions) {
n_samples <- nrow(dataset)
n_vars <- ncol(dataset)

vars_with_interaction <- sample(n_vars, n_interactions)

# Specify pairwise interactions
pair_inds <- matrix(nrow = 0, ncol = 2)
for (var in vars_with_interaction) {
pair <- c(var, sample(setdiff(vars_with_interaction, var), 1))
pair_inds <- rbind(pair_inds, pair)
}

# Generate a matrix of the interaction terms and their values
interaction_matrix <- matrix(0, nrow = n_samples, ncol = 0)
for (i in seq_len(n_interactions)) {
pair <- pair_inds[i, ]
interaction <- dataset[, pair] * dataset[, pair]
interaction_matrix <- cbind(interaction_matrix, interaction)
colnames(interaction_matrix)[i] <- paste0("interaction_", i)
}

# Generate the outcome variable based on the sum over the interactions
# (some noise is added too)
dataset <- as.data.frame(dataset)
dataset$outcome <- rowSums(interaction_matrix) + rnorm(n_samples) # Model on the original dataset (no interaction terms) & return R2 model <- lm(outcome ~ ., data = dataset) summary(model)$adj.r.squared
}

# Generate the datasets
norm_data <- matrix(rnorm(100000 * 100), nrow = 100000, ncol = 100)
binom_data <- matrix(rbinom(100000 * 100, 2, 0.5), nrow = 100000, ncol = 100)

set.seed(123)

# Normal data
sim_model(norm_data, 50) #  0.0002497078

# Binomial data
sim_model(binom_data, 50) #  0.9062993


Unless something's gone terribly wrong, you should see that with the binomial dataset the $$R^2$$ is ~0.91 and ~0 for the standard normal dataset.

So the question is: Why would this be the case? Why should the main effects be so good at predicting the outcome in the binomial case?

• Note you can include pictures in posts, which will make it easier for readers to get your point quickly. Jul 11 at 15:12
• The code is extremely complex. I was hoping to see a single model (no simulation) with two covariates to show your point. Jul 12 at 11:54
• @FrankHarrell, I am inclined to agree and have provided one in my answer below. Jul 12 at 13:29
• @FrankHarrell Sorry about that. This was in the context of a debate w/ my advisor about a project we are working on so I already had to work hard to abstract out domain specific logic. With a bit longer thinking about this as a pure stats problem I could indeed conjure simpler e.g.s: r a_bin <- rbinom(n, 1, 0.5) b_bin <- rbinom(n, 1, 0.5) o_bin <- a_bin * b_bin summary(lm(o_bin ~ a_bin + b_bin))$adj.r.squared # 0.658 a_norm <- rnorm(n) b_norm <- rnorm(n) o_norm <- a_norm * b_norm summary(lm(o_norm ~ a_norm + b_norm))$adj.r.squared # 0.000  Jul 14 at 7:07
• You are modeling the interaction term instead of using an interaction term (+ main effects) to model an outcome. So the result is not surprising. Product terms are collinear with their component variables. Also use a binary model not lm. You can make the interaction orthogonal to the main effects and the $R^2$ goes to zero in your sense. The rnorm version of what you did centers the variables so everything is orthogonal. Center the binary variable and you'll get the same result. But not sure why this is interesting. Jul 14 at 11:53

First, indeed, we would expect zero (adjusted) $$R^2$$ when the slope coefficients of the regression are zero. Hence, when $$R^2$$ is far from zero, that must be related to nonzero coefficient estimates. Since the true main effects are zero, nonzero estimates thereof ought to stem from omitted variable bias.

In a misspecified regression model, the omitted variable bias is related to the covariance of the included regressors with the omitted ones (see e.g. https://www.econometrics-with-r.org/6-1-omitted-variable-bias.html for the case of a simple regression)

In the standard normal case, we have, for $$x_1$$ (same logic for $$x_2$$) $$Cov(x_1,x_1x_2)=E(x_1^2x_2)-E(x_1)E(x_1x_2)=E(x_1^2x_2)=E(x_1^2)E(x_2)=0,$$ using the mean-zero property of the standard normal and independence of the regressors. Hence, in this case, there does not seem to be OVB.

In turn, the binomial regressors do not have mean zero (take the moments from https://en.wikipedia.org/wiki/Binomial_distribution#Properties), so that the above argument that included and omitted regressors have zero covariance does not go through. Hence, the coefficient estimates on the main effect will converge to something else than their true value of zero, leading to nonzero $$R^2$$.

Hence, the issue is less related to the particular distribution than to the mean zero property of your normal regressors. You "switch on" the normal non-zero mean regressors in my (or your, of course) snippet below to illustrate this point.

In a bit more detail, we regress on (let us omit a constant from the regression for simplicity; as we know, that makes computing an $$R^2$$ a bit questionable, but does not affect the general idea of the argument) $$X=(x_{1i}\quad x_{2i})_{i=1,\ldots,n}$$ such that our OLS estimator becomes $$\check\beta=(X'X)^{-1}X'y$$ Plugging in the DGP $$y_i=\beta x_{1i} x_{2i}+u_i\equiv\beta \tilde x_{i}+u_i$$ gives, in matrix notation, $$\check\beta=(X'X)^{-1}X'(\tilde X\beta+u)=(X'X)^{-1}X'\tilde X\beta+(X'X)^{-1}X'u$$ This will, assuming a law of large numbers to hold and exogeneity assumptions on the errors $$u_i$$, converge to $$\check\beta=\left(\frac{X'X}{n}\right)^{-1}\frac{X'\tilde X}{n}\beta+\left(\frac{X'X}{n}\right)^{-1}\left(\frac{X'u}{n}\right)\to_p\Sigma_{x_1,x_2}^{-1}\Sigma_{x,\tilde x}\beta,$$ where $$\Sigma_{x_1,x_2}=E(xx')=\begin{pmatrix}E(x_1^2)&E(x_1x_2)\\E(x_1x_2)&E(x_2^2)\end{pmatrix}$$ for $$x=(x_1,x_2)'$$ and $$\Sigma_{x,\tilde x}=E(x\tilde x)=\begin{pmatrix}E(x_1x_1x_2)\\E(x_2x_1x_2)\end{pmatrix}.$$

Now, in the case of independent N(0,1) regressors, $$E(xx')=\begin{pmatrix}E(x_1^2)&E(x_1)E(x_2)\\E(x_1)E(x_2)&E(x_2^2)\end{pmatrix}$$ is just the identity matrix and $$\Sigma_{x,\tilde x}=\begin{pmatrix}E(x_1x_1x_2)\\E(x_2x_1x_2)\end{pmatrix}=\begin{pmatrix}E(x_1^2x_2)\\E(x_1x_2^2)\end{pmatrix}=\begin{pmatrix}E(x_1^2)E(x_2)\\E(x_1)E(x_2^2)\end{pmatrix}=\begin{pmatrix}0\\0\end{pmatrix}$$ When (slightly adapting your example) regressors are independent Bernoulli with success probability $$p$$ (and hence expectation $$E(x_1)=E(x_2)=p$$), we have from the above link that $$E(x_1^2)=E(x_2^2)=p(1-p)+p^2=p$$ so that $$\begin{eqnarray*} \Sigma_{x_1,x_2}^{-1}&=&\frac{1}{p^2-p^4}\begin{pmatrix}p&-p^2\\-p^2&p\end{pmatrix}\\&=&\frac{1}{p-p^3}\begin{pmatrix}1&-p\\-p&1\end{pmatrix} \end{eqnarray*}$$ and $$\Sigma_{x,\tilde x}=\begin{pmatrix}p^2\\p^2\end{pmatrix}$$ Hence, $$\begin{eqnarray*} \check\beta&\to_p&\Sigma_{x_1,x_2}^{-1}\Sigma_{x,\tilde x}\beta\\ &=&\frac{1}{p-p^3}\begin{pmatrix}p^2-p^3\\p^2-p^3\end{pmatrix}\beta\\ &=&\frac{1}{(p+1)(p-p^2)}\begin{pmatrix}p(p-p^2)\\p(p-p^2)\end{pmatrix}\beta\\ &=&\frac{p}{p+1}\begin{pmatrix}1\\1\end{pmatrix}\beta \end{eqnarray*}$$

n <- 50000
p <- 0.5
x1 <- rnorm(n, mean = 2) # try these regressors instead of binomial below to see nonzero-mean normal regressors give same "story"
x2 <- rnorm(n, mean = 3)

x1 <- rbinom(n, 1, p)
x2 <- rbinom(n, 1, p)

y <- 2*x1*x2 + rnorm(n)

> summary(lm(y~x1+x2-1))

Call:
lm(formula = y ~ x1 + x2 - 1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
x1 0.668841   0.008413   79.50   <2e-16 ***
x2 0.672735   0.008389   80.19   <2e-16 ***
---

> 2*p/(p+1)
 0.6666667


In a similar way, you could analyze the plim in the normal case with $$E(x_j)=\mu\neq0$$ and $$E(x_j^2)=Var(x_j)+\mu^2$$ (and of course, if desired, allow for different $$p_j$$ or $$\mu_j$$ for the different regressors).

• Nice one, that's a great writeup. I originally took it to be that my binomial case is the strange one, but quick experiments show this is not true. Intuitively, it makes sense that the sum of two random variables will be in some sense correlated with the product. I then actually lose the intuition for why they oughtn't in the 0 mean case, even though I follow your calculation. BTW, my conclusion was that $Cov[(A \cdot B), (A + B)] =Var(A) \cdot E[B] + Var(B) \cdot E[A]$, implying that both means need to be 0 for there to be 0 covariance in this case. I assume our results are consistent though. Jul 14 at 7:01
• Thanks. Your result about the covariance of product with sum looks relevant indeed, and indeed points towards necessity. My analysis omits the constant, so there will be slight differences, but as you say, zero mean appears to play a central role. Jul 16 at 13:54