# 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 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 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. – Roland Jan 7 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 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 – Javier Fajardo Jan 15 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 at 21:58