As mentioned in the other answers, the existence of the random effects is irrelevant.
As a matter of R code, your function does not include
a*c. If that was intentional, be aware that
a*b*c is automatically expanded to:
a:b + a:c + b:c + a:b:c. That is, the term will be included. To deliberately exclude it, you should list the terms using the colon operator (
:). However, it is generally recommended not to leave out lower level terms, including lower level interaction terms, once a higher level term exists (see: Including the interaction but not the main effects in a model, and Do all interactions terms need their individual terms in regression model?). On the other hand, if that was an unintentional omission, note that your
a*b + b*c is redundant (the first term automatically expands, as just discussed).
You correctly note that "For each of a, b and c these are the expected increase in Y when that independent variable increases by 1 (assuming the other variables are unchanged)." However, it is important to recognize that "the other variables" cannot remain "unchanged", unless the other two main effects variables are exactly $0$. To understand this more fully, it may help to read my answer to: What does “all else equal” mean in multiple regression? Likewise, because you have a three-way interaction term in your model, the included two-way interactions have their clearest interpretation when the variable that is included in the three-way interaction but not included in the given two-way interaction is $0$. (For example,
a:b is the interaction between
At this point we can start to think about the interpretation of the beta for an interaction. As elsewhere, the standard interpretation is ceteris paribus ('all else held equal'). However, note that it is not possible to increment
a:b:c (or even
a:b) by $1$ without any of the other terms changing. This makes the situation difficult. There are some rules, of increasing complexity, for interpreting the sign of the beta as, say, superadditive ('more than the sum of its parts') or subadditive ('less than the sum of its parts'), etc., but these depend on the signs of other betas as well. I've never met anyone who has found those rules useful. You can compute the partial derivative and think of it as the rate of change in the slope / lower order interaction, if that helps you, but again, many people find that too much to wrap their heads around.
The best way is to try to interpret the model as a whole, and not the interaction term in isolation. To do that, solve for the simple effects of whichever variable has precedence in your mind at several values of the other interacting variables. You could then write out the equations of those lines, or you could plot them (cf., How to visualize a fitted multiple regression model?).