confidence / prediction ellipse
Hi Giuseppe You've posted a series of questions on this topic, never with any code you've tried or data, and usually with some undefined references to 'a' and 'b'. In spite of this, a variety of people have tried to give you helpful replies, intuiting what it is you might have meant. This is wasteful of everyone's time. If you want help from R-help, please make your effort to formulate a precise question, preferably with code and data. Otherwise, you might post to R-mindreaders if it ever gets established.
On 2/7/2013 11:20 AM, Giuseppe Amatulli wrote:
Hi Rolf, sorry for this late answer and thanks for your kind explanation and relative R code. I really appreciate. In reality the concept that I'm trying to address is a bit more complex. I'm fitting a model y vs 6 predictors with MARS / RandomForest / Multiple Linear Regression Models having 140 observations. I have the prediction of each model and would like to delineate the prediction ellipses for 3 models, for the 95% probability, and plotting them together with the observation vs prediction. I think that the prediction-ellipses code that you provide to me is valid also for predictions derived by not-linear model (such as MARS and RF). Is it correct? or should i use an alternative solution ? Moreover, I was expecting that the abline (lm(b,a)) would be correspond to the main axis of the prediction ellipse, but is not this the case. why? Thanks in advance Giuseppe On 28 January 2013 19:04, Rolf Turner <rolf.turner at xtra.co.nz> wrote:
I believe that the value of "radius" that you are using is incorrect. If you
have a data
matrix X whose columns are jointly distributed N(mu,Sigma) then a
confidence
ellipse for mu is determined by
n * (x - Xbar)' S^{-1}(x - Xbar) ~ T^2
where Xbar is the mean vector for X and S is the sample covariance matrix,
and where "T^2" means Hotelling's T-squared distribution, which is equal to
(n-1)*2/(n-2) * F_{2,n-2}
the latter representing the F distribution on 2 and n-2 degrees of freedom.
Thus (I think) your radius should be
radius <- sqrt(2 * (npts-1) * qf(0.95, 2, npts-2)/(npts*(npts-2)))
where npts <- length(a). Note that it is qf(0.95,2,npts-2) and *NOT*
qf(0.95,2,npts-1).
To get the corresponding *prediction* ellipse simply multiply the foregoing
radius by sqrt(npts+1). By "prediction ellipse" I mean an ellipse such that
the probability that a new independent observation from the same population
will fall in that ellipse is the given probability (e.g. 0.95). Note that
this does
not mean that 95% of the data will fall in the calculated ellipse (basically
because
of the *dependence* between S and the individual observations).
These confidence and prediction ellipses are (I'm pretty sure) valid under
the assumption that the data are (two dimensional, independent) Gaussian,
and that you use the sample covariance and sample mean as "shape" and
"centre". I don't know what impact your robustification procedure of using
cov.trob() will/would have on the properties of these ellipses.
A script which does the ellipses for your toy data, using the sample
covariance
and sample mean (rather than output from cov.trob()) is as follows:
#
# Script scr.amatulli
#
require(car)
a <- c(12,12,4,5,63,63,23)
b <- c(13,15,7,10,73,83,43)
npts <- length(a)
shape <- var(cbind(a, b))
center <- c(mean(a),mean(b))
rconf <- sqrt(2 * (npts-1) * qf(0.95, 2, npts-2)/(npts*(npts-2)))
rpred <- sqrt(npts+1)*rconf
conf.elip <- ellipse(center, shape, rconf,draw = FALSE)
pred.elip <- ellipse(center, shape, rpred,draw = FALSE)
plot(pred.elip, type='l')
points(a,b)
lines(conf.elip,col="red")
cheers,
Rolf Turner
On 01/27/2013 10:12 AM, Giuseppe Amatulli wrote:
Hi, I'm using the R library(car) to draw confidence/prediction ellipses in a scatterplot.
From what i understood the ellipse() function return an ellipse based
parameters: shape, center, radius . If i read dataEllipse() function i can see how these parameters are calculated for a confidence ellipse. ibrary(car) a=c(12,12,4,5,63,63,23) b=c(13,15,7,10,73,83,43) v <- cov.trob(cbind(a, b)) shape <- v$cov center <- v$center radius <- sqrt(2 * qf(0.95, 2, length(a) - 1)) # radius <- sqrt(dfn * qf(level, dfn, dfd)) conf.elip = ellipse(center, shape, radius,draw = F) plot(conf.elip, type='l') points(a,b) My question is how I can calculate shape, center and radius to obtain a prediction ellipses rather than a confidence ellipse? Thanks in Advance Giuseppe
-- Giuseppe Amatulli Web: www.spatial-ecology.net
Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. & Chair, Quantitative Methods York University Voice: 416 736-2100 x66249 Fax: 416 736-5814 4700 Keele Street Web: http://www.datavis.ca Toronto, ONT M3J 1P3 CANADA