Ggtcf

I'm quite an R newbie and facing the following challange.
I'll share my code here but applied to a different dataframe since I cannot share the original dataframe.
This is my code:

fit = gam( carb ~ te(cyl, hp, k=c(3,4)), data = mtcars)
plot(fit,rug=F,pers=T,theta=45,main="test")

using my company's data, this generates a nice surface with the predicted values on the Z axes.
I would like to add the actual response values as red dots on Z axis so that I could see where predicted values are under/over estimating the actual reponse.
Would you know what parameter I should add to plot in order to do that?
Many thanks

As @李哲源 pointed out in the comments, you shouldn't use plot here, because it's not flexible enough. Here's a version based on the referenced question Rough thin-plate spline fitting (thin-plate spline interpolation) in R with mgcv.

# First, get the fit
library(mgcv)
fit <- gam( carb ~ te(cyl, hp, k=c(3,4)), data = mtcars)

# Now expand it to a grid so that persp will work
steps <- 30
cyl <- with(mtcars, seq(min(cyl), max(cyl), length = steps) )

hp <- with(mtcars, seq(min(hp), max(hp), length = steps) )
newdat <- expand.grid(cyl = cyl, hp = hp)
carb <- matrix(predict(fit, newdat), steps, steps)

# Now plot it
p <- persp(cyl, hp, carb, theta = 45, col = "yellow")

# To add the points, you need the same 3d transformation
obs <- with(mtcars, trans3d(cyl, hp, carb, p))
pred <- with(mtcars, trans3d(cyl, hp, fitted(fit), p))
points(obs, col = "red", pch = 16)

# Add segments to show where the points are in 3d
segments(obs$x, obs$y, pred$x, pred$y)

That produces the following plot:

You might not want to make predictions so far from the observed data. You can put NA values into carb to avoid that. This code does that:

exclude <- exclude.too.far(rep(cyl,steps), 
 rep(hp, rep(steps, steps)), 
 mtcars$cyl, 
 mtcars$hp, 0.15) # 0.15 chosen by trial and error
carb[exclude] <- NA
p <- persp(cyl, hp, carb, theta = 45, col = "yellow")
obs <- with(mtcars, trans3d(cyl, hp, carb, p))
pred <- with(mtcars, trans3d(cyl, hp, fitted(fit), p))
points(obs, col = "red", pch = 16)
segments(obs$x, obs$y, pred$x, pred$y)

That produces this plot:

Finally, you might want to use the rgl package to get a dynamic graph instead. After the same manipulations as above, use this code to do the plotting:

library(rgl)
persp3d(cyl, hp, carb, col="yellow", polygon_offset = 1)
surface3d(cyl, hp, carb, front = "lines", back = "lines")
with(mtcars, points3d(cyl, hp, carb, col = "red"))
with(mtcars, segments3d(rep(cyl, each = 2), 
 rep(hp, each = 2), 
 as.numeric(rbind(fitted(fit),
 carb))))

Here's one possible view:

You can use the mouse to rotate this one if you want to see it from a different angle. One other advantage is that points that should be hidden by the surface really are hidden; in persp, they'll plot on top even if they should be behind it.