5

I am trying to plot a bunch of thin-plate spline response surfaces for measurements related to two continuous variables plus one discrete variable. So far, I have been subsetting the data based on the discrete variable to generate pairs of plots, but it seems to me there should be a way to create some slick trellised plots instead. It seems like this could be done by faceting heatmaps in ggplot2 with geom_tile and geom_contour, but I am stuck on

(1) how to reorganize the data (or interpret the predicted surface data) for plotting with ggplot2?

(2) syntax for creating trellised heatmaps with base graphics? or

(3) ways to use graphics from rsm to accomplish this (rsm can cope with higher-order surfaces, so I could coerce things to some extent, but plots aren't fully trellised).

Here's an example of what I've been working with so far:

library(fields)
library(ggplot2)

sumframe<-structure(list(Morph = c("LW", "LW", "LW", "LW", "LW", "LW", 
"LW", "LW", "LW", "LW", "LW", "LW", "LW", "SW", "SW", "SW", "SW", 
"SW", "SW", "SW", "SW", "SW", "SW", "SW", "SW", "SW"), xvalue = c(4, 
8, 9, 9.75, 13, 14, 16.25, 17.25, 18, 23, 27, 28, 28.75, 4, 8, 
9, 9.75, 13, 14, 16.25, 17.25, 18, 23, 27, 28, 28.75), yvalue = c(17, 
34, 12, 21.75, 29, 7, 36.25, 14.25, 24, 19, 36, 14, 23.75, 17, 
34, 12, 21.75, 29, 7, 36.25, 14.25, 24, 19, 36, 14, 23.75), zvalue = c(126.852666666667, 
182.843333333333, 147.883333333333, 214.686666666667, 234.511333333333, 
198.345333333333, 280.9275, 246.425, 245.165, 247.611764705882, 
266.068, 276.744, 283.325, 167.889, 229.044, 218.447777777778, 
207.393, 278.278, 203.167, 250.495, 329.54, 282.463, 299.825, 
286.942, 372.103, 307.068)), .Names = c("Morph", "xvalue", "yvalue", 
"zvalue"), row.names = c(NA, -26L), class = "data.frame")

sumframeLW<-subset(sumframe, Morph=="LW")
sumframeSW<-subset(sumframe, Morph="SW")

split.screen(c(1,2))
screen(n=1)
surf.teLW<-Tps(cbind(sumframeLW$xvalue, sumframeLW$yvalue), sumframeLW$zvalue, lambda=0.01)
summary(surf.teLW)
surf.te.outLW<-predict.surface(surf.teLW)
image(surf.te.outLW, col=tim.colors(128), xlim=c(0,38), ylim=c(0,38), zlim=c(100,400), lwd=5, las=1, font.lab=2, cex.lab=1.3, mgp=c(2.7,0.5,0), font.axis=1, lab=c(5,5,6), xlab=expression("X value"), ylab=expression("Y value"),main="LW plot")
contour(surf.te.outLW, lwd=2, labcex=1, add=T)
points(sumframeLW$xvalue, sumframeLW$yvalue, pch=21)
abline(a=0, b=1, lty=1, lwd=1.5)
abline(a=0, b=1.35, lty=2)

screen(n=2)
surf.teSW<-Tps(cbind(sumframeSW$xvalue, sumframeSW$yvalue), sumframeSW$zvalue, lambda=0.01)
summary(surf.teSW)
surf.te.outSW<-predict.surface(surf.teSW)
image(surf.te.outSW, col=tim.colors(128), xlim=c(0,38), ylim=c(0,38), zlim=c(100,400), lwd=5, las=1, font.lab=2, cex.lab=1.3, mgp=c(2.7,0.5,0), font.axis=1, lab=c(5,5,6), xlab=expression("X value"), ylab=expression("Y value"),main="SW plot")
contour(surf.te.outSW, lwd=2, labcex=1, add=T)
points(sumframeSW$xvalue, sumframeSW$yvalue, pch=21)
abline(a=0, b=1, lty=1, lwd=1.5)
abline(a=0, b=1.35, lty=2)

close.screen(all.screens=TRUE)
Henrik
  • 65,555
  • 14
  • 143
  • 159
rebeccmeister
  • 300
  • 2
  • 12

2 Answers2

3

As noted in a comment, melt() can be used to reshape the Tps() output, then it can be reformatted a bit (to remove NA's), recombined into a single data frame, and plotted. Here are plots with ggplot2 and levelplot:

library(reshape)
library(lattice)

LWsurfm<-melt(surf.te.outLW)
LWsurfm<-rename(LWsurfm, c("value"="z", "Var1"="x", "Var2"="y"))
LWsurfms<-na.omit(LWsurfm)
SWsurfms[,"Morph"]<-c("SW")

SWsurfm<-melt(surf.te.outSW)
SWsurfm<-rename(SWsurfm, c("value"="z", "X1"="x", "X2"="y"))
SWsurfms<-na.omit(SWsurfm)
LWsurfms[,"Morph"]<-c("LW")

LWSWsurf<-rbind(LWsurfms, SWsurfms)

LWSWp<-ggplot(LWSWsurf, aes(x,y,z=z))+facet_wrap(~Morph)
LWSWp<-LWSWp+geom_tile(aes(fill=z))+stat_contour()
LWSWp

ggplot2 image

or: levelplot(z~x*y|Morph, data=LWSWsurf, contour=TRUE)

lattice levelplot image

rebeccmeister
  • 300
  • 2
  • 12
  • Note that the x and y axes are no longer correct, however. There's probably a straightforward way to recalibrate these, based on the original axes and realizing they've been re-scaled between 0 and 80. – rebeccmeister Sep 20 '13 at 14:41
2
require(rgl)
open3d()
plot3d
surface3d(surf.te.outSW$x, surf.te.outSW$y, surf.te.outSW$z, col="red")
surface3d(surf.te.outLW$x, surf.te.outLW$y, surf.te.outLW$z, col="blue")
decorate3d()
      rgl.snapshot("OutRGL.png")

enter image description here

Another version where I scaled the x and y values by a factor of 10 and rotated to "look through" the gap. If this were your choice you might want to look at ?scaleMatrix

enter image description here

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • That's useful in terms of direct comparisons on the same plot, with better perspective on relative shape; I'll take a look at the package. The crude convention for print publication-worthy images has been colored, flat images, so I'm still after 2-D answers. – rebeccmeister Sep 18 '13 at 22:00
  • 1
    Then `lattice::contourplot` or `lattice::levelplot` would be my advice. – IRTFM Sep 18 '13 at 22:03
  • Any tips/ideas for how to translate the `Tps()` output (list) for plotting with `contourplot` or `levelplot`? It does look nice, as you've noted here: [methods-for-doing-heatmaps](http://stackoverflow.com/questions/7851602/methods-for-doing-heatmaps-level-contour-plots-and-hexagonal-binning). I am thinking if I can reorganize the two sets of output and append with the discrete variable I may be able to answer my own question. The z-values generated are an 80x80 matrix with NA's, and this seems to be generating errors. – rebeccmeister Sep 19 '13 at 15:52
  • I got single plots without difficulty. Getting the data into an 80 x 80 x 2 array to use the array method of handling grouping for posed some problems in that the axes couldn't get done automagically. I suppose plotting separately and then using grid.arrange would have been the "simple method" for me. – IRTFM Sep 19 '13 at 17:35
  • It looks like `melt()` from `reshape2` gets this done, with `na.omit` to clean up the resulting data - nearly automagic. I should be able to follow up soon with a ggplot2 version. – rebeccmeister Sep 19 '13 at 19:00
  • Nice work. I would not have thought that removing NA would be needed. Usually R plotting functions just leave NA spots blank. – IRTFM Sep 19 '13 at 20:19