[R] Contour Plot on a non Rectangular Grid

[David Winsemius]

On Oct 25, 2010, at 3:41 AM, Lorenzo Isella wrote:

> 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))

# Sorry. There was a single line omitted:

options(datadist="ddf")

> 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
>
> ############################################################

If you do not need a statistical approach to the contour plot, then  
Jim Lemon's more recent suggestion this morning may be more economical.

-- 
David.