[R] rgl: ellipse3d with axes

[John Fox]

Dear Michael,

Maybe this is irrelevant, since you appear to have a satisfactory solution
now, but here's an approach (from a figure that I drew in a recent book)
that computes the axes directly. This example is in 2D but I think that it
wouldn't be hard to generalize it:

------- snip --------

library(car)
library(MASS)

set.seed(12345)
Sigma <- matrix(c(1, .8, .8, 1), 2, 2)
Z <- mvrnorm(50, mu=c(0,0), Sigma=Sigma, empirical=TRUE)

eqscplot(Z, axes=FALSE, xlab="", ylab="")
abline(h=0, v=0)

ellipse(c(0,0), Sigma, radius=1, center.pch=FALSE, 
    col="black", segments=500)
    
E <- eigen(Sigma)
L <- E$vectors
lam <- sqrt(E$values)

lines(c(1, -1)*lam[1]*L[1,1], c(1, -1)*lam[1]*L[2,1], lwd=2)
lines(c(1, -1)*lam[2]*L[1,2], c(1, -1)*lam[2]*L[2,2], lwd=2)

------- snip --------

Some notes: eqscplot() from the MASS package does the analog of what Duncan
mentioned -- using equal units for both horizontal and vertical axes;
ellipse() from the car package draws the ellipse by transforming a circle,
but the axes are drawn directly using the eigenvalues and vectors of the
covariance (here, correlation) matrix.

Regards,
 John


------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
> Behalf Of Michael Friendly
> Sent: September-24-08 9:24 AM
> To: R-Help
> Subject: [R] rgl: ellipse3d with axes
> 
> Last week I asked about data ellipses with rgl:::ellipse3d() with lines
> showing the principal axes.
> (The goal is a visual demonstration of PCA as a rotation of variable
> space to component space.)
> I was trying, unsuccessfully, to use princomp() to generate the PCA axes
> and plot them using
> segments3d:
> 
> > > PC <- princomp(trees)
> > > sdev <- PC$sdev         # component standard deviations
> > > sd <- sqrt(diag(cov))   # variable standard deviations
> > >
> > > # vectors in variable space of principal components
> > > vec <- matrix(mu,3,3, byrow=TRUE) + diag(sd) %*% PC$loadings
> > >
> > > for (j in 1:3) {
> > >  mat <- rbind(mu, vec[j,])
> > >  segments3d(mat, col="red")
> > > }
> However, it occurred to me that these axes are just the orthogonal axes
> of the unit sphere
> that is transformed (using chol()) in ellipse3d, so plotting the axes
> transformed in the
> same way would give me what I want.
> 
> Looking at the result returned by ellipse3d, I see a normals component,
> but I'm not sure if this
> represents what I want, or, if it is, how to use it to draw the ellipse
> major axes in the plot.
> 
>  > e1 <-ellipse3d(cov, centre=mu, level=0.68)
>  > str(e1)
> List of 6
>  $ vb           : num [1:4, 1:386] 4.95 2.64 2.03 1.00 6.74 ...
>  $ ib           : num [1:4, 1:384] 1 195 99 196 51 197 99 195 27 198 ...
>  $ primitivetype: chr "quad"
>  $ homogeneous  : logi TRUE
>  $ material     : list()
>  $ normals      : num [1:4, 1:386]  0.290 -0.902 -0.320  1.000  0.635 ...
>  - attr(*, "class")= chr "qmesh3d"
> 
> -Michael
> 
> --
> Michael Friendly     Email: friendly AT yorku DOT ca
> Professor, Psychology Dept.
> York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Street    http://www.math.yorku.ca/SCS/friendly.html
> Toronto, ONT  M3J 1P3 CANADA
> 
> ______________________________________________
> 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.