Polynomial regression in R with multiple independent variables

Question

I want to do a polynomial regression in R with one dependent variable y and two independent variables x1 and x2. In my mind the model should look as follows,

y=b0 + b1x1+ b2x2+ b3x1^2+ b4x2^2+ b5x1x2

I tried lm(y~x1+x2+poly(x1,2,raw=TRUE)+poly(x2,2,raw=TRUE)) and also lm(y~x1+x2+I(x1^2)+I(x2^2)). But this only gives the squares and not the product of the two variables.

I can of cause do lm(y~x1+ x2+ x1^2+ x2^2+ x1x2). But i would like to know whether there is a much easier way than to write the whole equation out. I also would like to do to the power of 3 and 4 models which is more lengthy.

Answer 1

you can use polym

y ~ polym(x1, x2, degree=2, raw=TRUE) # is equivalent to
y ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2

But be careful with the order of the coefficients they are not the same as the second formula.

If you use degree=3 then it will add interactions of higher order like this I(x1^2):x2 +I(x2^2):x1, thus you have to adapt your formula.

NB : polym is a wrapper for poly, so you can use this latter with the same call.