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

David Winsemius dwinsemius at comcast.net
Wed Dec 4 19:30:45 CET 2013


On Dec 4, 2013, at 8:56 AM, Duncan Murdoch wrote:

> On 04/12/2013 11:36 AM, Jun Shen wrote:
>> Hi,
>> 
>> I have a dataset with two independent variables (x, y) and a response
>> variable (z). I was hoping to generate a response surface by plotting x, y,
>> z on a three dimensional plot. I can plot the data with rgl.points(x, y,
>> z). I understand I may not have enough data to generate a surface. Is there
>> a way to smooth out the data points to generate a surface? Thanks a lot.
> 
> There are many ways to do that.  You need to fit a model that predicts z from (x, y), and then plot the predictions from that model.
> An example below follows yours.
>> 
>> Jun
>> 
>> ===========================
>> 
>> An example:
>> 
>> x<-runif(20)
>> y<-runif(20)
>> z<-runif(20)
>> 
>> library(rgl)
>> rgl.points(x,y,z)
> 
> Don't use rgl.points, use points3d() or plot3d().  Here's the full script:
> 
> 
> x<-runif(20)
> y<-runif(20)
> z<-runif(20)
> 
> library(rgl)
> plot3d(x,y,z)
> 
> fit <- lm(z ~ x + y + x*y + x^2 + y^2)
> 

 Newcomers to R may think they would be getting a quadratic in x and y. But R's formula interpretation will collapse x^2 to just x and then it becomes superfluous and is discarded. The same result is obtained with z ~ (x + y)^2).   I would have thought that this would have been the code:

fit <- lm(z ~ poly(x,2) +poly(y,2) + x:y )


> 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, newdata=df)
> 
> surface3d(xnew, ynew, df$z, col="red")

With the modified fitting formula one sees a nice saddle (for that particular random draw) using rgl.snapshot().


The result with the earlier formula is a more restrained:




Continued thanks to you Duncan for making this great tool available.

-- 
David.


> Duncan Murdoch
>> 

David Winsemius
Alameda, CA, USA



More information about the R-help mailing list