Scatter3D add multiple surfaces

I am trying to add multiple surfaces to one scatter3d plot in R. Here is an example from scatter3d documentation:

library(plot3D)

with (mtcars, {
  
  # linear regression
  fit <- lm(mpg ~ wt + disp)
  
  # predict values on regular xy grid
  wt.pred <- seq(1.5, 5.5, length.out = 30)
  disp.pred <- seq(71, 472, length.out = 30)
  xy <- expand.grid(wt = wt.pred, 
                    disp = disp.pred)
  
  mpg.pred <- matrix (nrow = 30, ncol = 30, 
                      data = predict(fit, newdata = data.frame(xy), 
                                     interval = "prediction")[,1])
  
  # fitted points for droplines to surface
  fitpoints <- predict(fit) 
  
  scatter3D(z = mpg, x = wt, y = disp, pch = 18, cex = 2, 
            theta = 20, phi = 20, ticktype = "detailed",
            xlab = "wt", ylab = "disp", zlab = "mpg",  
            surf = list(x = wt.pred, y = disp.pred, z = mpg.pred,  
                        facets = NA, fit = fitpoints),
            main = "mtcars")
  
})

which produces the following plot:

Let's say I wanted to add two more surfaces that are :

surf1<-mpg.pred+5
surf2<-mpg.pred-5

How would I go about doing that?

You can use surf3D with add = TRUE to add the extra surfaces:

library(plot3D)

with (mtcars, {
  
  # linear regression
  fit <- lm(mpg ~ wt + disp)
  
  # predict values on regular xy grid
  wt.pred <- seq(1.5, 5.5, length.out = 30)
  disp.pred <- seq(71, 472, length.out = 30)
  xy <- expand.grid(wt = wt.pred, 
                    disp = disp.pred)
  
  mpg.pred <- matrix (nrow = 30, ncol = 30, 
                      data = predict(fit, newdata = data.frame(xy), 
                                     interval = "prediction")[,1])
  
  # fitted points for droplines to surface
  fitpoints <- predict(fit) 
  
  p1 <- scatter3D(z = mpg, x = wt, y = disp, pch = 18, cex = 2,
               theta = 20, phi = 20, ticktype = "detailed",
               xlab = "wt", ylab = "disp", zlab = "mpg",
               surf = list(x = wt.pred, y = disp.pred, z = mpg.pred,
                           facets = NA, fit = fitpoints),
               main = "mtcars")
  p2 <- surf3D(x = matrix(wt.pred, 30,30),
            y = matrix(disp.pred, 30,30,byrow=TRUE), z = mpg.pred + 5, facets = NA, add = TRUE)
  p3 <- surf3D(x = matrix(wt.pred, 30,30),
            y = matrix(disp.pred, 30,30,byrow=TRUE), z = mpg.pred - 5, facets = NA, add = TRUE)
})

^{Created on 2022-04-27 by the reprex package (v2.0.1)}

The legend and color schemes are a little bit wrong; you could perhaps fix them, but I don't know how. Using colkey=list(plot=FALSE) in the two surf3D calls helps a bit, but the color scheme is still wrong for the extra planes.