[R] contour for plotting confidence interval on scatter plot of bivariate normal distribution

[drflxms]

Thanx a lot Greg for the hint and letting me not alone with this!

Tried ellipse and it works well. But I am looking for something more
precise. The ellipse fits the real border to the nearest possible
ellipse. I want the "real" contour, if possible.

Meanwhile I found an interesting function named draw.contour
http://quantitative-ecology.blogspot.com/2008/03/plotting-contours.html
provided by "Forester", an "Assistant Professor in the Department of
Fisheries". - I love it how statistics brings people from completely
different specialities together. - I am a neurologist.

This function results not exactly in the same curve, I got with the code
I posted last, but is in deed very close by (parallel). So I think both
solutions are probably true and differ only because of rounding and
method (..!?).

Still I am not completely happy with

1. the numeric solution of setting to 68% to get 1SD. I'd prefer a
symbolic, formula driven solution instead.
2. me not comprehending the code completely.

While 2. will be hopefully solved soon by me delving into the code, 1
remains...

Best, Felix


Am 03.03.12 17:28, schrieb Greg Snow:
> Look at the ellipse package (and the ellipse function in the package)
> for a simple way of showing a confidence region for bivariate data on
> a plot (a 68% confidence interval is about 1 SD if you just want to
> show 1 SD).
> 
> On Sat, Mar 3, 2012 at 7:54 AM, drflxms <drflxms at googlemail.com> wrote:
>> Dear all,
>>
>> I created a bivariate normal distribution:
>>
>> set.seed(138813)
>> n<-100
>> x<-rnorm(n); y<-rnorm(n)
>>
>> and plotted a scatterplot of it:
>>
>> plot(x,y)
>>
>> Now I'd like to add the 2D-standard deviation.
>>
>> I found a thread regarding plotting arbitrary confidence boundaries from
>> Pascal Hänggi
>> http://www.mail-archive.com/r-help@r-project.org/msg24013.html
>> which cites the even older thread
>> http://tolstoy.newcastle.edu.au/R/help/03b/5384.html
>>
>> As I am unfortunately only a very poor R programmer, the code of Pascal
>> Hänggi is a myth to me and I am not sure whether I was able to translate
>> the recommendation of Brain Ripley in the later thread (which provides
>> no code) into the the correct R code. Brain wrote:
>>
>> You need a 2D density estimate (e.g. kde2d in MASS) then compute the
>> density values at the points and draw the contour of the density which
>> includes 95% of the points (at a level computed from the sorted values
>> via quantile()). [95% confidence interval was desired in thread instead
>> of standard deviation...]
>>
>> So I tried this...
>>
>> den<-kde2d(x, y, n=n) #as I chose n to be the same as during creating
>> the distributions x and y (see above), a z-value is assigned to every
>> combination of x and y.
>>
>> # create a sorted vector of z-values (instead of the matrix stored
>> inside the den object
>> den.z <-sort(den$z)
>>
>> # set desired confidence border to draw and store it in variable
>> confidence.border <- quantile(den.z, probs=0.6827, na.rm = TRUE)
>>
>> # draw a line representing confidence.border on the existing scatterplot
>> par(new=TRUE)
>> contour(den, levels=confidence.border, col = "red", add = TRUE)
>>
>> Unfortunately I doubt very much this is correct :( In fact I am sure
>> this is wrong, because the border for probs=0.05 is drawn outside the
>> values.... So please help and check.
>> Pascal Hänggis code seems to work, but I don't understand the magic he
>> does with
>>
>> pp <- array()
>> for (i in 1:1000){
>>        z.x <- max(which(den$x < x[i]))
>>        z.y <- max(which(den$y < y[i]))
>>        pp[i] <- den$z[z.x, z.y]
>> }
>>
>> before doing the very same as I did above:
>>
>> confidencebound <- quantile(pp, 0.05, na.rm = TRUE)
>>
>> plot(x, y)
>> contour(den, levels = confidencebound, col = "red", add = TRUE)
>>
>>
>> My problems:
>>
>> 1.) setting probs=0.6827 is somehow a dirty trick which I can only use
>> by simply knowing that this is the percentage of values inside +-1sd
>> when a distribution is normal. Is there a way doing this with "native"
>> sd function?
>> sd(den.z) is not correct, as den.z is in contrast to x and y not normal
>> any more. So ecdf(den.z)(sd(den.z)) results in a percentile of 0.5644 in
>> this example instead of the desired 0.6827.
>>
>> 2.) I would like to have code that works with any desired confidence.
>> Unfortunately setting probs to the desired confidence would probably be
>> wrong (?!) as it relates to den.z instead of x and y, which are the
>> underlying distributions I am interested in. To put it short I want the
>> confidence of x/y and not of den.z.
>>
>>
>> I am really completely stuck. Please help me out of this! Felix
>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> 
>