0

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

PMouton
  • 1
  • 4
  • In the formula try `poly(x, y, degree = 2, raw = TRUE)`. – Rui Barradas May 02 '21 at 11:22
  • @RuiBarradas: Thank you for the suggestion. It gives different coefficients and the formula when using raw=TRUE does seem to work. What I still don't understand however, is that the constructed trend surface in the example above is correct, but that somehow the coefficients cannot translated into the proper equation that is the source of this trend surface. – PMouton May 02 '21 at 13:42
  • Without `raw=TRUE` the polynomials returned, quoting the documentation, *are all orthogonal to the constant polynomial of degree 0*. What you need is the *raw* polynomials. – Rui Barradas May 02 '21 at 14:38
  • @RuiBarradas: I edited my question to show that there is a difference between the raw=T and raw=F trend surfaces. As I like the raw=F surface better, is there a way to still make sense of the coefficients that accompany this model? – PMouton May 02 '21 at 14:52

1 Answers1

0

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

Elia
  • 2,210
  • 1
  • 6
  • 18
  • Thanks for your answer. As showed in the edited question, I like the results better when 'raw = F'. If I use your function p, the results show strange numbers not corresponding to the output of the lm+poly function. How should I interpret the coefficients when poly(..., raw=F)? – PMouton May 03 '21 at 06:38