# Creating Optimal “Smoothed” Contour-Type Plots (X,Y,Z) Using Mathematica 8.04

### Update: data here

Here is a contour plot produced in the SAS program:

In this plot, variable 23405 forms the X axis, variable 20201 forms the Y axis, and outcome variable 50502 forms a pseudo Z axis via color and contour lines. As one can readily see, the red patch marks the "sweet spot" as far as combining 23405 and 20201.

THE DATA: X and Y variables are continuous, with the range shown on both axes. The original Z data is continuous as well. However what is shown in the plot are not the original data. Z was "massaged" or "smoothed," based on the original Z data, so as to allow for the creation of the smooth-looking contour plot. In this sense, the surface shown itself does not pass through (precisely) the original Z points, but is an estimation. (Analagous to a single linear regression line not passing through all the data points.) This approximating process done by the particular program used (SAS vs. Mathematica vs. other) is what I'm particularly curious about. Bottom line -- will the supposed "sweet spot" arrived at be any different?

Using the same data, I would like to create a similar -- or even better -- contour plot in Mathematica. At the very least, to compare it with the one from SAS, to see if the conclusions appear to be the same.

There seem to be various ways of creating such plots in Mathematica, a few of the methods just recently added to the program in version 8.

Some possible ways in Mathematica:

And other ways as well.

I'm hoping someone here has had experience creating such contour plots in Mathematica, and can share insights as to which approach has worked well for him/her. Also, any added tricks or parameters needed to arrive at an optimal output?

-
Can you define "optimal"? Also, are the data available anywhere? –  Verbeia Jul 29 '12 at 8:09
I think the down voters have responded a bit too harshly to this question. Yes, it could benefit from some additional clarity re: "optimal" (perhaps "quantifying" the term in some way); a description of the structure of his data; and showing something the OP has tried in Mathematica. Still, the OP has provided a good example from SAS and has identified a few different ways of doing this in Mathematica. For my part I'd rather see questions like this (which I think do appear thoughtful) gently nudged towards greater clarity rather than down voted and possibly discouraged. –  Jagra Jul 29 '12 at 14:06
nk, some further description of the data and your assumptions would help. There are many contouring algorithms. Perhaps the most fundamental distinction is between algorithms for regular (or gridded) data and algorithms for "irregular" data. Yours could be the latter if there are missing values; otherwise, it's probably the former. Next, there's the question of how smooth you think the variation of the Z variable really is. SAS seems to think it's incredibly smooth, but is that really the case? The more you can say about these variables, the more specific and applicable the answers can be. –  whuber Jul 29 '12 at 20:25
@nk One perhaps does not have to have full access to the data, but nevertheless I would answer your query of 13 hours ago with an unqualified, emphatic, yes: the characteristics of the data strongly suggest and limit the appropriate methods of contouring them. This is a rich subject which we can't possibly get into within comments, but I can assure you my position is well-established and not something idiosyncratic--although I recognize, from the existence of many popular automatic contouring packages, that this position is not very well known or appreciated, either. –  whuber Jul 30 '12 at 13:06
nkormanik, thanks for providing the data. I put it on pastebin.com and linked it in your question. Please do not use rapidgator again as it uses pornographic advertisements. –  Mr.Wizard Jul 31 '12 at 13:14

This is a provisional answer, since without a definition of “optimal”, it’s hard to know what you are looking for.

The choice between ListContourPlot and SmoothDensityHistogram depends on the data you have.

• If you have lists of triplets, $x$, $y$ and $z$, then use ListContourPlot.
• If you have lists of $x$ and $y$ value pairs, some of which may be repeated or close to each other, use SmoothDensityHistogram.

If you had x-y pairs and wanted to convert them into x-y-z space, then you could use Tally to do this. Because the output will then be in the form {{{x1, y1}, 3},{{x2,y2}, 2}...} you will need to Map a function like Flatten or Join to get it into the form you want (there are several serviceable ways to do this, but that is not really germane to this question).

You could also use SmoothKernelDistribution to convert x-y pairs to the equivalent of an x-y-z data set.

In terms of getting “nice” output, you will need to consider several parameters and the options that access them:

-
SmoothDensityHistogram not appropriate for this problem. X,Y,Z in data file. But need to smooth Z. –  nkormanik Jul 31 '12 at 4:18
@all Yes, triplets. Assume an .xls spreadsheet with three columns, X,Y,Z. What would be the code to bring it in and use ListContourPlot? Wonder if it will look similar to the plot at top of page here. Curious to see it. Thanks for sticking with this little project. –  nkormanik Aug 3 '12 at 0:09

This has really turned into a statistical question, but here are some quick insights Mathematica can provide.

First, let's import the data (from a downloaded copy) and plot their supports (that is, their point locations in $(x,y)$ coordinates):

nk = Import["F:/temp/nk.mtp"];
{y, x, z} = First /@ nk /. nk (* Note the x-y reversal to match the map *);
data = Transpose[{x, y, z}];
ListPlot[Most /@ data]


The sparseness of the points at the right (large $x$ value) and their concentrations at the top and bottom (extreme $y$ values) indicate first, that any interpolation in the middle right will be uncertain; and second, that preliminary re-expressions of the coordinates are likely to help in the analysis. To accomplish this, let's look at their distributions:

TableForm[{Histogram /@ {x, y, z}}]


There are rigorous ways to identify good re-expressions, but experience and a little experimentation suggest a cube root will work well for the first and third variables and a scaled arcsine for the second. The objective is to obtain distributions that are somewhere between uniform and symmetric unimodal:

TableForm[{Histogram /@ {x0 = x^(1/3), y0 = ArcSin[(y - 5.5)/110], z0 = z^(1/3)}}]


That's about as good as we could hope for, absent any indication of what these variables mean or what their natural ranges might be. From now on, the analyses will be conducted in terms of the re-expressed variables. (The results can always be cast back into the original variables, because all three re-expressions are monotonic, one-to-one transformations.)

With these (conventional, necessary) preliminaries out of the way, let's start looking at these data in three dimensions. But no contouring yet: we need to see the data.

pointPlot = ListPointPlot3D[Transpose[{x0, y0, z0}],
ColorFunction -> Function[{x, y, z}, Hue[z, .8, .8]], AxesLabel -> {"x", "y", "z"}]


How about that! The contour plot has papered over a truly large and messy scatter of the data.

This plot has been saved for later: it will show up again below.

What we need now is some heavy-duty statistical modeling. One of the best approaches (ultimately) would be a Generalized Additive Model. I don't believe Mathematica is quite at that level, so at this point--just to make a little progress--I dropped the data into R for further analysis:

Export["F:/temp/nk.csv", Transpose[{x0, y0, z0}]];


For the record, in R I used a quick-and-dirty local smoother that acts much like a GAM; here are the commands to get the data in, process them, and write them back out for further work in Mathematica.

data <- read.table("f:/temp/nk.csv", sep=",", col.names=c("x","y","z"))
fit <- loess(z ~ x + y, data=data, span=.25)
write.table(fit$fitted, "f:/temp/fit.csv", row.names=FALSE, col.names=FALSE)  To display more detail, I reduced the default value of span from$0.75$to$0.25$: this causes the smooth to follow the "local wiggles" in the data more closely. (Usually, one experiments with the amount of smoothing to explore possible patterns in the data; the default is only a start, not a guide.) Let's get this back into Mathematica and visualize it. Because the smooth is a set of points (with the same support as the original data) intended to represent a continuous surface, a different approach to visualization is appropriate. fit = Transpose[Import["F:/temp/fit.csv"]] // First; fitPlot = ListSurfacePlot3D[Transpose[{x0, y0, fit}], ColorFunction -> Function[{x, y, z}, Hue[z, .8, .8]], MeshStyle -> Opacity[.1], PlotStyle -> Opacity[.8], AxesLabel -> {"x", "y", "z"}]  This is one of infinitely many ways to approximate, interpolate, or otherwise model these data. It makes little sense on its own: it needs to be understood in the context of the original points. We saved a plot of them earlier, so let's show the data together with their smooth: Show[pointPlot, fitPlot]  It is now evident that the apparent rise in$z$values for larger$x$is relatively small compared to the scatter of the$z\$ values themselves. As we saw at the outset, this rise is not supported by many points at all: it would be prudent to view it as merely suggestive of an actual trend and possibly as a spurious artifact.

There are many directions in which this analysis could now go, but they depend on knowing what these data mean and the purposes of the analysis. Because these are not mentioned in the question (which is probably appropriate, because this site is not devoted to statistical questions), here is a good place to stop. I hope this short introductory tour of some of Mathematica's visualization capabilities has helped with the present question and will help future readers with similar inquiries.

-
A closely related question appeared on stats an hour ago: stats.stackexchange.com/questions/33402/…. Some good ideas are already emerging in the comments and answers. –  whuber Jul 31 '12 at 14:15
Really top-notch analysis. I'll read and re-read it over many times. And carefully follow the steps in Mathematica. Surely your treatment here will be helpful to others as well. –  nkormanik Aug 1 '12 at 1:15
My only additional requests would be to show how the smoothing can be done DIRECTLY within Mathematica 8.04, as best one can, without using any other program (R, etc.). And to produce a contour plot similar to the one shown at top, so that a comparison can be made. –  nkormanik Aug 1 '12 at 1:16
Your insights into what's going on behind the data are great. I especially appreciate the line, "This is one of infinitely many ways to approximate, interpolate, or otherwise model these data." Thank you for pointing this out. The "sweet spot" is so elusive. –  nkormanik Aug 1 '12 at 1:16
@nk, I am very sympathetic with your request for a solution entirely within Mathematica, but I just don't know of any one that is handy. And why code one? Because you are already working in SAS, you could do all this analysis directly in SAS. Mathematica is great for simulation, visualization, and theoretical analyses while SAS (and other stats packages) shine in the areas of data management, specialized capabilities, and support for those capabilities (such as diagnostic statistics and plots, help files, and so on). You might as well capitalize on the strengths of each. –  whuber Aug 1 '12 at 15:32