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?
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
$\endgroup$lm
. You can make the interaction orthogonal to the main effects and the $R^2$ goes to zero in your sense. Thernorm
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. $\endgroup$