Thanks to all those who replied to my query (don't you lot have anything
better to do on a Friday afternoon?). The consensus seems to be that I
should fit a kernel to the posteriors, and draw the contours around that
(the details are below). I'll now go away and have a play....
The replies in detail were:
Deb Niemeier:
Hi - the only way I have found to do this is using Mathematica.
Simon Jackman:
you will need to generate a 2d density estimate using the MCMC output as
"data", and then use contour, or Rob Hydnman's HDR (highest density
regions) routines, available from StatLib. for the density estimation,
try ksmooth or Clive Loader's locfit.
to see an example, check out
http://jackman.stanford.edu/fig3.ps
from a recent paper, also available from my web site.
Martyn Plummer:
Get the KernSmooth package for R and use the function bkde2D() to get a
two-dimensional density estimate for your output. Then use the
contour()
function to plot the limits of the confidence regions (they are not
constrained to be ellipses).
The only hard part is calculating the right level to plot for the
contour function. I'll give you that on Monday, but I want to go
home now.
You may also want to play with different bandwidths for the kernel
density estimate.
Jeff Rondinone:
I think this can be done with a spreadsheet like Microsoft Excel, if
you
have the patience.
Condense your posterior density output to a modest (maybe 25 - 50 bins)
histogram for each dimension. Make a square array of the joint
posterior by
lining up one histogram in a row and the second in a column, and
programming
each cell in the resulting array to be the product of the row and
column.
With the posterior arranged like this, Excel has a simple tool for
making a
contour plot in a number of formats.
I'm certain there are more elegant solutions, but my software budget is
limited.
John Sprague:
The best solution would be Trellis contour plots in S-PLUS or R
since you
have that running. The code for these displays is available on the
math-soft site where all of the code that Cleveland uses for the Trellis
(and other) figures is available along with the data as well.
Presumably
much of this is R compatible.
John P. Burkett:
If R is like most mathematical software, it can plot functions but not
correspondences
that give two y values for each x value. Such software can be used to
draw
an ellipse if
you split the ellipse into upper and lower contours. To identify the
left
and right join
points, take the derivative of x w.r.t. y, set it to zero, and compute
the
two roots.
Plot the upper contour between the join points; repeat for the lower
contour.
P.Malewski:
for an example look at "kde2d" in library "MASS":
library(MASS)
example(kde2d)
Again, many thanks.
Bob
--
Bob O'Hara
Metapopulation Research Group
Division of Population Biology
Department of Ecology and Systematics
PO Box 17 (Arkadiankatu 7)
FIN-00014 University of Helsinki
Finland
tel: +358 9 191 7382 fax: +358 9 191 7301
email: [log in to unmask]
To induce catatonia, visit:
http://www.helsinki.fi/science/metapop/
And a Happy New 1900 to you all!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|