I have the classic multiple regression model $Z=β_0+β_1X+β_2Y+β_3XY$ that produces output with a significant interaction term. There's a lot of advice out there (and the formal principle of marginality) that says I shouldn't interpret the main (marginal) effects with a significant interaction term. However, in my case both X and Y are continuous variables and I'm struggling to understand what the main effect terms mean in this case and why I shouldn't interpret them.
I do understand why interpreting the main effect (if an interaction is significant) doesn't make sense for categorical data: the "average category" is meaningless because categorical data isn't ordered and discussing an "average bin" is nonsensical. With continuous variables though, it does seem like the main effect term is still a meaningful value because continuous variables have meaningful averages.
This seems to hold up under simulation as well - I can specify a model with both main effect and interaction terms, run a linear regression on it, and recover reasonable estimates for the parameters of both the main effect and the interaction effect. R code below. This implies that the main effect does have some meaning and can be interpreted if the variables are continuous. I do also realize that I cannot report the main term alone, and my discussion of this would be something like "the effect of X alone was estimated as <value> but varied between <upper> and <lower> depending on its interaction with Y". Is this an acceptable way to conclude things from the main term even though there's a significant interaction effect? Or is there something I'm missing?
set.seed(123)
x_vals <- (-10):10
y_vals <- (-10):10
mr_data <- expand.grid(x_vals, y_vals)
mr_data <- setNames(mr_data, c("x_vals", "y_vals"))
beta_1 <- 2
beta_2 <- 3
beta_3 <- 1
mr_data$z_vals <- beta_1*mr_data$x_vals + # Main effect of x
beta_2*mr_data$y_vals + # Main effect of y
beta_3*mr_data$x_vals*mr_data$y_vals + # Interaction effect
rnorm(441) # Some noise so lm() doesn't complain
lm(z_vals~x_vals*y_vals, data = mr_data) # Coefs are ~same as betas above
Coefficients:
(Intercept) x_vals y_vals x_vals:y_vals
0.01216 1.99236 2.99894 0.99953