[R] How to generate a smoothed surface for a three dimensional dataset?

[Duncan Murdoch]

On 05/12/2013 10:33 AM, Jun Shen wrote:
> Hi Federico/Duncan/David/Bert,
>
> Thank you for your thoughtful comments and it's a great learning
> experience. I can see the critical point here is to find a right function
> to make the prediction. So I was thinking to start with "loess". However
> the predict.loess gave me an error as follows
>
> Error in `$<-.data.frame`(`*tmp*`, "z", value = c(0.417071766265867,
> 0.433916401753023,  :
>    replacement has 20 rows, data has 400
>
> Here is the code I tried. Thank you for your help again!
>
> Jun
> =====================================
>
> x<-runif(20)
> y<-runif(20)
> z<-runif(20)
>
> library(rgl)
> plot3d(x,y,z)
>
> loess(z~x+y,control=loess.control(surface='direct'),span=.5,degree=2)->fit.loess
>
> xnew <- seq(min(x), max(x), len=20)
> ynew <- seq(min(y), max(y), len=20)
>
> df <- expand.grid(x = xnew, y = ynew)
>
> df$z<-predict(fit.loess,newdata=df)

After the error, use traceback() to find which function called 
`$<-.data.frame`.  It shows that it was your final assignment

df$z<-predict(fit.loess,newdata=df)


which causes the error, because the predict function returns a matrix.  
So you can get the plot using

surface3d(xnew, ynew, predict(fit.loess,newdata=df), col="gray")

You may want

aspect3d(1,1,1)

afterwards; loess isn't so good at extrapolation.  Or you may want to 
set predictions to NA outside the convex hull of your data.  (I'm not 
sure which function is easiest to calculate that, but here's one way:

hullx <- x[chull(x,y)]
hully <- y[chull(x,y)]
keep <- sp::point.in.polygon(df$x, df$y, hullx, hully)
znew <- predict(fit.loess,newdata=df)
znew[!keep] <- NA
plot3d(x,y,z)
surface3d(xnew, ynew, znew, col="gray")
aspect3d(1,1,1)