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)
1
  • I was thinking that pkg:rgl would be a nice environment for this effort.
    – IRTFM
    Sep 18, 2013 at 21:22

2 Answers 2

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

1
  • 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. Sep 20, 2013 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

6
  • 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. Sep 18, 2013 at 22:00
  • 1
    Then lattice::contourplot or lattice::levelplot would be my advice.
    – IRTFM
    Sep 18, 2013 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. 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. Sep 19, 2013 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, 2013 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. Sep 19, 2013 at 19:00

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.