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

[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.
>



-- 
Gregory (Greg) L. Snow Ph.D.
538280 at gmail.com