Interpreting main effect with significant interaction term in continuous by continuous multiple regression 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



With continuous variables though, it does seem like the main effect term is still a meaningful value because continuous variables have meaningful averages.

A problem is that continuous predictors might not be centered around their mean values when introduced into the model. The main effect” reported for X is typically its association with outcome when Y = 0 in the data presented to the model, and vice-versa. Thus re-centering Y affects the apparent "main effect" of X, and vice-versa. From the modeling perspective there thus isn't a unique "meaningful value" for the "main effect" of a continuous predictor involved in interaction terms. See this page for a worked-through example.
If you keep that in mind and center your data so that the mean of each continuous predictor equals 0, you can then interpret a "main effect" coefficient as the association with outcome when the interacting predictors are at their mean values. With an interaction, however, the effects of predictors on each other's associations with outcome and scenarios other than associations at mean values may be of more interest.
