# R formula for higher order polynomials and interactions, only allowing polynomial of degree 1 to interact

I am trying to build a (mixed) model using several predictor variables, and including some interactions and potentially higher degree polynomial versions of the continuous variables. The model formula looks like this:

y ~ factor * poly(x, 2)   # where x is a continuous variable


However, the resulting model ends up being quite complex in terms of coefficients, because in addition to the factor level I also have a complex structure of nested random effects.

There, I was thinking about only testing the interaction with x of degree 1, not with the quadratic term. I would write that formula as this:

y ~ factor * x + poly(x, 2)   # where x is a continuous variable


However, when I run lmer() using that formula, I get the following warning:

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient


I believe that the warning is produced because the formula states twice to use x as a predictor -- once as part of the interaction, and once as part of the polynomial term.

I am aware that the warning() does not necessarily imply that this is problematic. But, is it in this case?

Also, I am curious about whether there is a clean way to write the formula where I can achieve what I was looking for (including quadratic terms, but not as part of the interaction). Thanks!

• Is there some reason why you are using polynomials in this model? Restricted cubic splines can fit well and models using them can generalize better to new samples.
– EdM
Jan 6, 2020 at 15:04
• Just drop the first column of poly(x,2). You should be able to do that with poly(x,2)[,-1].
– whuber
Jan 6, 2020 at 17:12
• @whuber I think factor : x + poly(x, 2) would be even simpler. Although, I'm not sure mixing raw and orthogonal polynomials is sensible. Jan 7, 2020 at 7:01
• I have reopened this question because I believe it was mistakenly closed: although couched in terms of R formulas, it concerns statistical issues, as Dimitris Rizopoulos' answer indicates.
– whuber
Jan 7, 2020 at 18:21
• @EdM I did not think about restricted cubic splines, just because I am more familiar with polynomials. But I have been exploring them in more details since your suggestion, and they seem a better approach, as also suggested in the response that I have accepted below Jan 15, 2020 at 9:49

The reason why you get this warning is indeed because the term factor * x expands to factor + x + factor:x, and poly(x, 2) is equivalent (but not the same because it uses orthogonal polynomials) with x + I(x^2). Hence, you have two times the linear term for x. As a solution, you could try specifying the following formula y ~ factor + poly(x, 2) + factor:x.

With that being said, it is not very logical to only use the interaction with the linear term. For example, say that factor is sex with levels male and female. Then you postulate that for males and females the quadratic effect is the same but the linear effect is different. This is a quite strong assumption for which you should have some prior indication that it holds.

As suggested in one of the comments, it would be better to work with splines instead of polynomials. These are more stable numerically and they can also quite easily be used in R. For example, if you load the splines package, you can use specify the formula as y ~ factor * ns(x, 2). This will assume a restricted cubic spline for x with one internal knot.

• Re "not very logical:" I think it's pretty common to see models of the form y ~ x + z + x*z + x^2, but such models are exactly the same as the one in the question (merely expressed differently).
– whuber
Jan 6, 2020 at 21:58
• @whuber Could you explain whether it is appropriate to use models in the form y ~ x + z + x*z + x^2 despite the "not very logical" reason suggested by @DimitrisRizopoulos? A reviewer suggested I should also include the interaction x^2*z if the polynomial x^2 has been included. May 23, 2022 at 14:02
• @Niamh The degree of a polynomial has an intrinsic meaning in the sense that it cannot increase whenever the variables are subjected to a linear transformation. Introducing a degree-three term, like $x^2z,$ into a degree-two model thereby fundamentally changes it mathematical nature. One usually introduces terms of the next higher degree only for diagnostic purposes such as checking goodness of fit. The reviewer implicitly takes the position that $x^2$ is just another variable, but that ignores its polynomial relationship with $x,$ which becomes less tenable as more terms are included.
– whuber
May 23, 2022 at 14:09
• (continued) My extended analysis at stats.stackexchange.com/a/408855/919 situates these ideas within a standard framework (orthogonal polynomials), thereby grounding them in conventional practice. The part of especial interest for you begins partway through at "At this point, parts of the univariate program break down... ."
– whuber
May 23, 2022 at 14:10
• @whuber I'm sure this is a perfectly accurate answer but wondering if could you simplify further (non-mathematical explanation)? Are you saying it's generally acceptable to include the x*z interaction and not x^2*z because all sequential polynomials are orthogonal to one another? In practise, does this mean that the interaction term in y ~ x + z + x^2 + x*z is "similar enough" to the interaction term here y ~ x + z + x^2 + x*z^2 ? May 23, 2022 at 16:27