3D plot from lm model with a nonlinear parameter

Question

I am trying to plot a 3Dplot of a linear regression with two variables.

My linear model:

ls = lm(mpg ~ disp + qsec, data = mtcars)

The minimum and maximum values to predict:

mmd = min(mtcars$disp):max(mtcars$disp)
mmq = min(mtcars$qsec):max(mtcars$qsec)

The linear function:

flm = function(x1, x2) ( (ls$coefficients[1]) + (ls$coefficients[2] * x1) + (ls$coefficients[3] * x2)) 

The outer matrix:

t = outer(X = mmd, Y = mmq, FUN = flm)

The plot:

persp3D(z = t, ticktype = "detailed")

Now I would like to control for the non-linearity of the qsec parameter.

Is this the right way to do it?

ls = lm(mpg ~ disp + poly(qsec, degree = 2), data = mtcars)
mmd = min(mtcars$disp):max(mtcars$disp)
mmq = min(mtcars$qsec):max(mtcars$qsec)

flm = function(x1, x2) ( (ls$coefficients[1]) + (ls$coefficients[2] * x1) +
                           (ls$coefficients[3] * x2) 
                         + (ls$coefficients[4] * x2)
                         ) 

t = outer(X = mmd, Y = mmq, FUN = flm)
persp3D(z = t, ticktype = "detailed")

Answer 1

If I understand correctly, shouldn't ls$coefficients[4] * x2 be ls$coefficients[4] * x2 * x2?

These details could be easily avoided if instead of writing your function to evaluate the fit with new data, you used predict.lm. This is what I think the code should look like, making other small improvements I will comment about below:

disp.density <- 100L
qsec.density <- 100L
disp.range <- with(mtcars, seq(from       = min(disp),
                               to         = max(disp),
                               length.out = disp.density))
qsec.range <- with(mtcars, seq(from       = min(qsec),
                               to         = max(qsec),
                               length.out = qsec.density))
data.grid <- expand.grid(disp = disp.range,
                         qsec = qsec.range)

fit <- lm(mpg ~ disp + qsec, data = mtcars)
predicted.vec <- predict(ls, data.grid)
predicted.mat <- matrix(predicted.vec, nrow = disp.density,
                                       ncol = qsec.density,
                                   dimnames = list(disp.range,
                                                   qsec.range))
library(???)
persp3D(z = predicted.mat, ticktype = "detailed")

And the same code should still work if you replace the formula mpg ~ disp + qsec with mpg ~ disp + poly(qsec, degree = 2). Some more notes:

1. It is recommended you avoid variable names like ls and t since they are very basic functions from R. Or you run the risk of shadowing them.
2. disp and qsec are not integers so using min(mtcars$disp):max(mtcars$disp) is a bit unpredictable. I feel it is better to use seq like I did, and with control over the grid density.
3. I was not able to test my code because you did not mention the package name for persp3D so it is possible it will need some adjustment. In the future, please make sure to give us everything we need to run your code.