Plotting a 3D surface in R

Question

I have a lm of the form lm(log(z) ~ x*y).

I can make a 3D scatterplot of the data using s3d <- scatterplot3d(x,y,log(z)).

How do I plot the surface of my regression function?

I know I can do s3d$plane3d(lm(log(z) ~ x + y) for a simple linear model, but this doesn't work for nonplanar surfaces; e.g. s3d$plane3d(lm(log(z) ~ x * y) won't work.

Answer 1

So first, you are much more likely to get help if you provide a sample dataset (or your actual data, which would be better) and a reproducible example of your code. Otherwise you are forcing us to do it - which I think is why your question was ignored for 6 hours. See this link for more details.

Second, using lm(...) like that, with formula=log(z)~x*y is an exceptionally bad idea. Linear modeling is based on the assumption that the errors in the response (log(z) in your case) are normally distributed with constant variance. If your z-data has normally distributed error with constant variance, then log(z) certainly will not. This is a classic mistake; the right way to do this is using generalized linear modeling (see the glm package, among others), with family=poisson.

Finally, to your question. The code below creates a 3D scatter plot overlaid onto a response surface. It uses the rgl package, which generates rotatable 3D plots.

Here I've colored the points so that those below the surface are red and those above are green, and added dropline from each point to the surface.

The tricky bit is that in surface3d(...), the x and y arguments are vectors corresponding to a grid, and z is a matrix with one row for every x value and one column for every y value.

For your real data you may need to tweak the scale=... argument in open3d(...).

# create sample dataset - you have this already,,,
set.seed(1)   # for reproducible example
df <- data.frame(x=sample(1:50,50)/50,y=sample(1:50,50)/50)
df$z <- with(df,exp(4*x + 2*y - 6*x*y + 6)+rnorm(50,sd=500))
fit <- lm(log(z) ~ x*y,df)   # bad, bad, bad - don't do this!!!

# you start here...
df$pred <- predict(fit)
# set up matrix of z-values
x <- seq(min(df$x),max(df$x),len=100)
y <- seq(min(df$y),max(df$y),len=100)
plot.df <- expand.grid(x=x,y=y)
plot.df$z <- predict(fit,newdata=plot.df)
library(reshape2)
z <- dcast(plot.df,x~y,value.var="z")[-1]

# plot the points, the fitted surface, and droplines
library(rgl)
colors <- 2.5+0.5*sign(residuals(fit))
open3d(scale=c(1,1,0.2))
points3d(df$x,df$y,log(df$z),col=colors)
surface3d(x,y,as.matrix(z),col="blue",alpha=.2)
apply(df,1,function(row)lines3d(rep(row[1],2),rep(row[2],2),c(log(row[3]),row[4]),col=colors))
axes3d()
title3d(xlab="X",ylab="Y",zlab="log(Z)")