[R] Confidence ellipse for correlation

Prof Brian Ripley ripley at stats.ox.ac.uk
Tue Oct 28 20:20:43 CET 2003


On Tue, 28 Oct 2003, Tanya Murphy wrote:

> Hello,
> 
> SAS' point and click interface has the option of produce a scatterplot with a 
> superimposed confidence ellipse for the correlation coefficient. Since I 
> generally like R so much better, I would like to reproduce this in R. I've 
> been playing with the ellipse package. In order to have the points and the 
> ellipse on the same graph I've done the following. 
> (Load ellipse package...)
> > data(Puromycin)
> > attach(Puromycin)
> > my<-mean(rate)
> > mx<-mean(conc)
> > sdy<-sd(rate)
> > sdx<-sd(conc)
> > r<-cor(conc,rate)
> > plot(ellipse(r,scale=c(sdx,sdy),centre=c(mx,my)),type='l')
> > points(conc,rate)
> 
> 1) Is my use of 'scale' and 'centre' theoretically correct?

Depends on whose theory you have in mind!  This is not `a confidence
ellipse for the correlation coefficient', as confidence ellipses are for
pairs of parameters, not variables.  It seems to be a plot of a contour of
the fitted bivariate normal.

> 2) Is there a more efficient way to get the 5 parameters? (I guess I could 
> write a little function, but has it already been done?)

You could do things like
mxy <- mean(Puromycin[c("rate", "conc")])
sxy <- sapply(Puromycin[c("rate", "conc")], sd)

> The non-linear relationship between these variables brings up another point: 
> Is there a way to plot a contour (empirical?) containing, say, 95% of the 
> values.

Yes.  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()).

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595




More information about the R-help mailing list