[R] Contour Plot on a non Rectangular Grid

[Lorenzo Isella]

On 10/25/2010 01:32 AM, David Winsemius wrote:

> You were advised to look at rms. Why have you dismissed this suggestion?
> Using your data setup below and packaging into a dataframe.
>
> require(rms)
> ddf <- datadist(xysf <- as.data.frame(xys))
> olsfit <- ols(V3~rcs(V1,3)+rcs(V2,3), data=xysf)
>
> bounds <- perimeter(xysf$V1, xysf$V2)
> plot(xysf$V1, xysf$V2) #demonstrates the extent of the data
> bplot(Predict(olsfit, V1,V2), perim=bounds) # a levelplot is the default
>
> bplot(Predict(olsfit, V1,V2), perim=bounds, lfun=contourplot)
> bplot(Predict(olsfit, V1,V2), perim=bounds, lfun=contourplot,
> xlim=c(-2.5,2.5))
> # to demonstrate that perimeter works
>
> # and as expected this shows very little variability d/t V1
> olsfit # note that
>  > anova(olsfit)
> Analysis of Variance Response: V3
>
> Factor d.f. Partial SS MS F P
> V1 2 0.01618738 8.093691e-03 19.47 <.0001
> Nonlinear 1 0.01618738 1.618738e-02 38.93 <.0001
> V2 2 470.67057254 2.353353e+02 566040.95 <.0001
> Nonlinear 1 470.67057254 4.706706e+02 1132081.91 <.0001
> TOTAL NONLINEAR 2 527.78127558 2.638906e+02 634723.80 <.0001
> REGRESSION 4 527.78127558 1.319453e+02 317361.90 <.0001
> ERROR 7663 3.18594315 4.157566e-04
> # most the the regression SS is in the V2 variable
> # Q.E.D.

Thanks David,
But I am experiencing some problems with your snippet.
When I run the code at the end of the email (saved as plot_circular.R), 
I get the following error

 > source('plot_circular.R')
Error in value.chk(at, which(name == n), NA, np, lim) :
   variable V1 does not have limits defined by datadist

which you clearly do not have on your machine. Have I left out some bits 
of your code?

Lorenzo

############################################################

require(rms)

R <- pi/2

n <- 100

x <- y <- seq(-R,R, length=n)

xys <- c()

temp <- seq(3)

for (i in seq(n)){

for (j in seq(n))

#check I am inside the circle
   if ((sqrt(x[i]^2+y[j]^2))<=R){

temp[1] <- x[i]
temp[2] <- y[j]
temp[3] <- abs(cos(y[j]))
xys <- rbind(xys,temp)

   }


}



  ddf <- datadist(xysf <- as.data.frame(xys))
  olsfit <- ols(V3~rcs(V1,3)+rcs(V2,3), data=xysf)

  bounds <- perimeter(xysf$V1, xysf$V2)
  plot(xysf$V1, xysf$V2) #demonstrates the extent of the data
  bplot(Predict(olsfit, V1,V2),  perim=bounds)  # a levelplot is the default

  bplot(Predict(olsfit, V1,V2),  perim=bounds, lfun=contourplot)
  bplot(Predict(olsfit, V1,V2),  perim=bounds, lfun=contourplot, 
xlim=c(-2.5,2.5))
# to demonstrate that perimeter works

############################################################