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

[Bert Gunter]

Jun:

You need to:

1. Read the docs
2. Consult a local statistician, as you seem to be out of your depth
statistically here.

Re: 1:

?predict.loess explicitly says:

"If newdata was the result of a call to expand.grid, the predictions
(and s.e.'s if requested) will be an **array** of the appropriate
dimensions." (emphasis added)

So change your assignment to:

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


## or maybe if plot3d accepts a matrix argument (I don't do 3d plots),
more directly:

plot3d(predict(fit.loess,newdata = df))  ## May Not Work!!

Re: 2

Your comment that:

" I can see the critical point here is to find a right function to
make the prediction. "

is what indicates to me that your "critical point" is that you have
insufficient knowledge and need help. Feel free to disagree, of
course.

Cheers,
Bert

On Thu, Dec 5, 2013 at 7:33 AM, Jun Shen <jun.shen.ut at gmail.com> 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)
>
>
>
> On Wed, Dec 4, 2013 at 4:37 PM, Bert Gunter <gunter.berton at gene.com> wrote:
>>
>> ... or, more simply
>>
>> lm(z ~ polym(x,y, degree=2) )
>>
>> ?polym
>>
>> Cheers,
>> Bert
>>
>>
>>
>>
>>
>> On Wed, Dec 4, 2013 at 10:30 AM, David Winsemius <dwinsemius at comcast.net>
>> wrote:
>> >
>> > 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
>> >
>> > ______________________________________________
>> > R-help at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> > http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>>
>>
>>
>> --
>>
>> Bert Gunter
>> Genentech Nonclinical Biostatistics
>>
>> (650) 467-7374
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

(650) 467-7374