[R] rgl: ellipse3d with axes

[Duncan Murdoch]

On 24/09/2008 12:32 PM, Michael Friendly wrote:
> Thanks Duncan (& others)
> 
> Here is a function that does what I want in this case, and tries to do 
> it to work generally
> with ellipse3d.  (Note that I reverse the order of centre and scale 
> 'cause I was bitten
> by trying ellipse3d.axes(cov, mu))
> 
> # draw axes in the data ellipse computed by ellipse3d
> ellipse3d.axes <-
> function (x, centre = c(0, 0, 0), scale = c(1, 1, 1), level = 0.95,
>     t = sqrt(qchisq(level, 3)), which = 1:3, ...)
> {
>     stopifnot(is.matrix(x)) # should test for square, symmetric
>     cov <- x[which, which]
>     eigen <- eigen(cov)
>     # coordinate axes, (-1, 1), in pairs
>     axes <- matrix(
>       c(0, 0, -1,   0, 0, 1,
>         0, -1, 0,   0, 1, 0,
>        -1, 0, 0,    1, 0, 0),  6, 3, byrow=TRUE)
> 
>     # transform to PC axes
>     axes <- axes %*% sqrt(diag(eigen$values)) %*% t(eigen$vectors)
>     result <- scale3d(axes, t, t, t)
>     if (!missing(scale))
>         if (length(scale) != 3) scale <- rep(scale, length.out=3)
>         result <- scale3d(result, scale[1], scale[2], scale[3])
>     if (!missing(centre))
>         if (length(centre) != 3) scale <- rep(centre, length.out=3)
>         result <- translate3d(result, centre[1], centre[2], centre[3])
>     segments3d(result, ...)
>     invisible(result)
> }
> 
> Test case:
> 
> library(rgl)
> data(iris)
> iris3 <- iris[,1:3]
> cov <- cov(iris3)
> mu <- mean(iris3)
> col <-c("blue", "green", "red")[iris$Species]
> plot3d(iris3, type="s", size=0.5, col=col, cex=2, box=FALSE)
> plot3d( ellipse3d(cov, centre=mu, level=0.68), col="gray", alpha=0.2,  
> add = TRUE)
> 
> axes <- ellipse3d.axes(cov, centre=mu)
> 
> One thing I can't explain, compared to your example is why the my axes 
> extend outside the ellipse,
> whereas yours didn't.

That's just because you specified level in the ellipse3d call, but not 
in the ellipes3d.axes call.

One thing that looks a little strange is that the PC axes don't appear 
to be orthogonal:  this is because the scaling is not the same on all 
coordinates.  It might look better doing the first plot as

plot3d(iris3, type="s", size=0.5, col=col, cex=2, box=FALSE, aspect="iso")

Duncan Murdoch

> 
> One final remark-  I knew that axes %*% chol(cov) did not give the 
> orthogonal PC axes I wanted,
> but at least it gave me something on the right scale and location. But 
> these axes also turn out to be
> useful for visualizing multivariate scatter and statistical concepts.  
> chol() gives the factorization of
> cov that corresponds to the Gram-Schmidt orthogonalization of a data 
> matrix -- orthogonal axes
> in the order x1,  x2|x1,  x3|x1, x2, ..., and vector length and 
> orientation in this coordinate system
> correspond to Type I SS in linear models.
> Thus, I could see generalizing my ellipse3d.axes function further to 
> allow a type=c("pca", "chol")
> argument.
> 
> -Michael
> 
> 
> 
> Duncan Murdoch wrote:
>> That's easy, but it doesn't give you the principal axes of the 
>> ellipse.  Just use
>>
>> axes %*% chol(cov)
>>
>> If you start with a unit sphere, this will give you points on its 
>> surface, but not the ones you want.  For those you need the SVD or 
>> eigenvectors.  This looks like it does what you want:
>>
>> axes <- matrix(
>>     c(0, 0, 0, # added origin
>>        0, 0, -1,   0, 0, 1,
>>        0, -1, 0,   0, 1,  0,
>>        -1, 0, 0,   1, 0, 0),  7, 3, byrow=TRUE)
>> axes <- axes[c(1,2,1,3,1,4,1,5,1,6,1,7),]  # add the origin before each
>>
>> cov <- cov(trees)
>> eigen <- eigen(cov)
>> shade3d(ellipse3d(cov, t=1, alpha=0.2, col='red'))
>> segments3d(axes %*% sqrt(diag(eigen$values)) %*% t(eigen$vectors))
>>
>> Duncan Murdoch
> 
>