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
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.
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
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
OK, then let me repharase it. The axes of the unit sphere are like
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)
taken in pairs. I'd like to transform these coordinates the same was as
in ellipse3d() and add them to the plot.
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
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
OK, then let me repharase it. The axes of the unit sphere are like
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)
taken in pairs. I'd like to transform these coordinates the same was as
in ellipse3d() and add them to the plot.
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
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
OK, then let me repharase it. The axes of the unit sphere are like
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)
taken in pairs. I'd like to transform these coordinates the same was as
in ellipse3d() and add them to the plot.
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
Why did I add the origin? It's not needed, your axes were fine. Sorry.
Duncan Murdoch
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.
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
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
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
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.