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
Michael Friendly wrote:> 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" >The normals component contains the surface normals. It is used to help in rendering the surface, but isn't much use for your purposes. Unfortunately, I'm not familiar enough with the internals of princomp to tell you how to get the axes you want. Duncan Murdoch> -Michael > >
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.