Say I have a model with one response variable ($y$, a continuous variable) and two predictor variables ($x_1$ & $x_2$, both continuous). My model includes both the additive effects of these variables and their interaction.
$$ y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \delta(x_1 \cdot x_2) +\epsilon $$
The interpretation of $\delta$ in this case (i.e. the interaction coefficient) is the change in slope of $\beta_2$ given a unit change in $x_1$ (correct me if I'm wrong).
However, the data I'm working with (and modelling for the purpose of prediction rather than inference) exhibit a more complex interaction than this. The slope of $\beta_2$ seems to change as a quadratic function of $x_1$. Is there any way to incorporate this into my linear model?
The dataset I'm working with is too large to paste here, and if I could easily reproduce the problem I probably wouldn't be asking the question. I expect there is a simple solution and lack of coffee is just getting the better of me...