2

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!

Dimitris Rizopoulos
  • 21,024
  • 2
  • 20
  • 47
Javier Fajardo
  • 123
  • 2
    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 06 '20 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 06 '20 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 07 '20 at 07: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 07 '20 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 '20 at 09:49
  • @Roland, that was my concern, I wouldn't like to mix them. It's the same it would happen with whuber suggestion to drop the first column (although thanks for the hint, that's a way of dropping an element in poly() that I had not considered! – Javier Fajardo Jan 15 '20 at 09:51

1 Answers1

2

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.

Dimitris Rizopoulos
  • 21,024
  • 2
  • 20
  • 47
  • 3
    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 06 '20 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. – regressless May 23 '22 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 '22 at 14:09
  • (continued) My extended analysis at https://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 '22 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 xz interaction and not x^2z 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 ? – regressless May 23 '22 at 16:27
  • @regress This has little or nothing to do with orthogonality. Indeed, what might be 'acceptable' ultimately must be determined by the nature of the data and the objective of the analysis. I am only pointing to some principles and considerations that would support not including $x^2z$ in your model, as well as to suggest there is no general reason why including $x^2$ and $z$ requires one also to include $x^2z.$ – whuber May 23 '22 at 17:24