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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
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@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.
______________________________________________
R-help@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.
