Suppose I have the following function, representing a linear model with an interaction term: $$ f(x, y) = \beta_{1} x + \beta_{2} y + \beta_{3} xy. $$
Now I want to see how the function changes if both $x$ and $y$ change. I can calculate it exactly by doing: $$ \Delta f = f(x + \Delta x, y + \Delta y) - f(x, y), \\[0.5em] \phantom{abcdefghijklm}= (\beta_{1} + \beta_{3}y)\Delta x + (\beta_{2} + \beta_{3}x)\Delta y + \beta_{3}\Delta x \Delta y. $$ How can I approximate this result using derivatives? If I apply the total derivative, i.e summing the partial derivatives I am missing the last term: $$ \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} = (\beta_{1} + \beta_{3}y) + (\beta_{2} + \beta_{3}x) $$ The last term only appears when adding $\partial f/\partial x \partial y$, but I don't understand why I'd have to do that.
Thanks in advance for any useful hints!