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

Question

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)).

Entirely optional code section ahead!

---

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[1]] * dataset[, pair[2]]
    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) # [1] 0.0002497078

# Binomial data
sim_model(binom_data, 50) # [1] 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?

Answer 1

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)
[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).