Interpretation of coefficients 2nd degree polynomial in R

Question

I made a trend surface using the lm and poly functions in R. The summary of my results looks as follows:

> summary(tls.lm)

Call:
lm(formula = DEM_small_domain_TLS_5cm ~ poly(x, y, degree = 2), 
    data = tlsDF)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.09543 -0.01378 -0.00087  0.01260  0.10404 

Coefficients:
                           Estimate Std. Error  t value Pr(>|t|)    
(Intercept)                2.08e+01   4.98e-05 416945.4   <2e-16 ***
poly(x, y, degree = 2)1.0 -1.51e+01   2.00e-02   -754.0   <2e-16 ***
poly(x, y, degree = 2)2.0  6.79e-01   2.00e-02     33.9   <2e-16 ***
poly(x, y, degree = 2)0.1 -2.01e+01   2.00e-02  -1007.3   <2e-16 ***
poly(x, y, degree = 2)1.1 -7.17e+02   8.04e+00    -89.2   <2e-16 ***
poly(x, y, degree = 2)0.2 -6.02e-01   2.00e-02    -30.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02 on 161446 degrees of freedom
Multiple R-squared:  0.908, Adjusted R-squared:  0.908 
F-statistic: 3.19e+05 on 5 and 161446 DF,  p-value: <2e-16

How do I translate the coefficients into an equation? I tried the following:

Intercept + coeff1.0 * x + coeff2.0 * x^2 +coeff1.1 * x * y + coeff 0.1 * y + coeff 0.2 * y^2

I also tried to replace x and y in the above formula, but the results from filling in these equations with some extracted x and y values are completely different from the (correctly constructed) trend surface result.

Edit In the image below the differences between the constructed trend surfaces for poly(..., raw=F) and poly(...,raw=T) are shown. It is recommended that I use raw=T, but I like the raw=F surface better. Is there still a way to make sense of the coefficients above? Poly_raw_difference

Answer 1

As suggested by @Rui you have to add raw=T as an argument to poly: here an example with the built-in dataset mtcars.

data("mtcars")
fit <- lm(formula = mpg ~ poly(disp, hp, degree = 2,raw = T), 
   data = mtcars)
summary(fit)
#coefficient of regression
coef <- coef(fit)
#use coefficient for predict outcome
p <- function(coef, x, y) {
  coef[1]+coef[2]*x+coef[3]*x^2+coef[4]*y+coef[5]*x*y+coef[6]*y^2
}

p1 <- round(p(coef,mtcars$disp,mtcars$hp),11)
p2 <- round(unname(predict(fit,mtcars[,3:4])),11)
identical(p1,p2)

the 2 vectors are identical up to 11 decimal places