Thanks again, Duncan. Please allow me to ask one more question. Is it
possible to generate a contour plot overlaying with the plot3d() plot?
Jun
On Thu, Dec 5, 2013 at 11:08 AM, Duncan Murdoch <murdoch.dun...@gmail.com>wrote:
> 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)
>
>
>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
