R/pca3d.R In pca3d: Three Dimensional PCA Plots

Documented in pca3d

```## functions starting with the dot are, by my convention, considered internal
## for this file and the one specific function that is exported from it.

## radius depends on the scale in each direction

for(i in 1:3) {
r <- 2*max(abs(pca.coords[ , i ]))
radius[i] <- radius[i] * r / 100
}
}

## for drawing circles in 3D, precalculate some values
.sin.t <- sin(seq(0, 2 * pi, len= 10))
.cos.t <- cos(seq(0, 2 * pi, len= 10))

## draw labels for each centroid
.group.labels <- function(coords, group, ...) {
centr.coords <- calc.centroids(coords, group)
texts3d(centr.coords[,1], centr.coords[,2], centr.coords[,3], levels(group), col= "black", cex= 2, ...)
}

## add titles of axes
.show.axe.titles <- function(r, show.scale, axes.color, components, axe.titles) {
if(missing(axe.titles) | is.null(axe.titles)) {
axe.titles = paste("PC", components)
}

cex <- 1.2

texts3d(
c(cex *r[ 2,1 ], 0, 0),
c(0, cex *r[ 2,2 ], 0),
c(0, 0, cex *r[ 2,3 ]), col= "black", axe.titles)

if(show.scale) {
axis3d('x',pos=c(NA, 0, 0), col= axes.color)
axis3d('y',pos=c(0, NA, 0), col= axes.color)
axis3d('z',pos=c(0, 0, NA), col= axes.color)
}
}

## add the xy plane
.show.plane <- function(r) {
rbind(
c(r[1,1],0,r[1,3]),
c(r[1,1],0,r[2,3]),
c(r[2,1],0,r[2,3]),
c(r[2,1],0,r[1,3])
),
col= "grey",
alpha= 0.2
)
}

## add the axes
.show.axes <- function(axes.color, ranges) {
axes <- rbind(
c(ranges[1,1], 0, 0),
c(ranges[2,1], 0, 0),
c(0, ranges[1,2], 0),
c(0, ranges[2,2], 0),
c(0, 0, ranges[1,3]),
c(0, 0, ranges[2,3])
)
segments3d(axes, col= axes.color)

scale <- c(1, 1, 1)
if(! missing(ranges)) {
scale <- ranges[2,]
radius <- max(scale) / 50
}

arrows3d(axes, radius= radius, scale= scale,  col= axes.color)
}

## Plot ellipses surrounding each group
.group.ellipses.3D <- function(coords, group, group.col, ellipse.ci=0.95) {
centr.coords <- calc.centroids(coords, group)
lg <- levels(group)
for(i in 1:length(lg)) {
g   <- lg[i]
sel <- group == g
s   <- cov(coords[sel,,drop=FALSE])
cc  <- centr.coords[i,]
#lines(ellipse(s, centre=cc), col=group.col[i])
shade3d(ellipse3d(s, centre=cc, level=ellipse.ci), col=group.col[i], alpha=0.2)
}
}

## add group centroids
.centroids.3D <- function(coords, group, radius, group.shape, group.col, shape.functions) {

centr.coords <- calc.centroids(coords, group)

for(l in levels(group)) {
shape.functions[[ group.shape[l] ]](centr.coords[l,,drop= F], col= group.col[l], radius= radius, alpha= 0.7)
n.l <- length(which(group == l))
tmp <- coords[ group == l, ]

tmp <- rbind(tmp, t(sapply(1:nrow(tmp), function(x) centr.coords[ l, ])))
tmp <- tmp[ as.vector(rbind(1:n.l, (n.l + 1):(2*n.l))), ]
segments3d(tmp, col= group.col[ l ])
}
}

## drop shadows on the xy plane
.draw.shadows.3D <- function(pca.coords, xz.functions, shape, radius) {

seg <- NULL
for(i in 1:nrow(pca.coords)) {
tmp <- pca.coords[ i, ]
seg <- rbind(seg, tmp)
tmp[ 2 ] <- 0
seg <- rbind(seg, tmp)
#squarexz(tmp[1], tmp[2] + 0.001, tmp[3], radius= radius, col= "grey")
xz.functions[[ shape[i] ]](tmp[1], 0 + 0.001 * radius[2], tmp[3], radius= radius, col= "grey")
xz.functions[[ shape[i] ]](tmp[1], 0 - 0.001 * radius[2], tmp[3], radius= radius, col= "grey")
}
segments3d(seg, col= "grey", alpha= 0.5)
}

.squarexz <- function(x, y, z, radius= c(1, 1, 1), col= "grey", alpha= 0.5) {
tmp <- NULL
for(i in 1:length(x)) {
tmp <- rbind(tmp,
c(x[i] - r[1], y[i], z[i] - r[3]),
c(x[i] - r[1], y[i], z[i] + r[3]),
c(x[i] + r[1], y[i], z[i] - r[3]),
c(x[i] + r[1], y[i], z[i] + r[3]),
c(x[i] - r[1], y[i], z[i] + r[3]),
c(x[i] + r[1], y[i], z[i] - r[3]))
}
triangles3d(tmp, col= col, alpha= alpha)
}

.trianglexz <- function(x, y, z, radius= c(1, 1, 1), col= "grey", alpha= 0.5) {
tmp <- NULL
for(i in 1:length(x)) {
tmp <- rbind(tmp,
c(x[i]                            , y[i], z[i] + radius[3]    ),
c(x[i] + sqrt(3) * radius[1] / 2, y[i], z[i] - radius[3] / 2),
c(x[i] - sqrt(3) * radius[1] / 2, y[i], z[i] - radius[3] / 2)
)
}
triangles3d(tmp, col= col, alpha= alpha)

}

.circlexz <- function(x, y, z, radius= c(1, 1, 1), col= "grey", alpha= 0.5) {

r <- radius * 2 / 3
n <- length(.sin.t)
xv <- x + r[1] * .sin.t
yv <- rep(y, n)
zv <- z + r[3] * .cos.t

tmp <- NULL
for(i in 1:(n-1)) {
tmp <- rbind(tmp,
c(x, y, z),
c(xv[i],   yv[i],   zv[i]  ),
c(xv[i+1], yv[i+1], zv[i+1]))
}

triangles3d(tmp, col= col, alpha= alpha)
}

.myspheres3d <- function(coords, radius= c(1, 1, 1), ...) {
}

## create a 3D biplot
.biplot.3D <- function(biplot.coords, biplot.vars, ranges= NULL) {

nvar <- length(biplot.vars)

coords <- matrix(0, ncol= 3, nrow= nvar * 2)
tips <- biplot.coords[ biplot.vars, ]
coords[ seq(2, nvar * 2, by= 2), ] <- tips

mx <- max(abs(coords))

coords <- coords / mx
scale <- c(1, 1, 1)
if(! missing(ranges)) {
scale <- ranges[2,]
coords <- t(apply(coords, 1, function(x) x * scale))
}

arrows3d(coords, col= "red", scale= scale)

if(length(rownames(biplot.coords)) > 0)
names <- rownames(biplot.coords)[biplot.vars]
else
names <- biplot.vars

coords <- coords[ seq(2, nvar * 2, by = 2), ]
texts3d(coords[,1], coords[,2], coords[,3], names, col= "black")
}

.rgl.window.init <- function(bg) {
cat("Creating new device\n")
rgl.open()
#open3d(antialias= 3)
par3d(windowRect= 50 + c(0, 0, 640, 640))
rgl.viewpoint(theta= 45, phi= 30, fov= 60, zoom= 1)
bg3d(bg)
rgl.pop("lights")
light3d(theta=-45, phi=15, specular= "black")
light3d(theta=55, phi=55,  specular= "black")
#light3d(theta=15, phi=-15,  specular= "white")
}

## initialize the rgl window, clear etc.
.pca_plot_init3d <- function(new, show.plane, show.axes, show.axe.titles, show.scale,
axes.color, axe.titles, main.ranges, components, bg) {

# open new window if requested or not opened yet
if(new || rgl.cur() == 0)
.rgl.window.init(bg)

par3d(skipRedraw= TRUE)
rgl.clear()

#r <- apply(pca.coords, 2, function(x) c(-max(abs(range(x))), max(abs(range(x)))))
if(show.plane) .show.plane(main.ranges)
if(show.axes) .show.axes(axes.color, main.ranges * 1.1)
if(show.axe.titles) .show.axe.titles(main.ranges, show.scale, axes.color, components, axe.titles)
if(show.scale) {
axis3d('x',pos=c(NA, 0, 0))
axis3d('y',pos=c(0, NA, 0))
axis3d('z',pos=c(0, 0, NA))
}
}

## finalize the plot
.pca_plot_finish_3d <- function() {
aspect3d(c(1, 1, 1))
#aspect3d("iso")
par3d(skipRedraw= F)
rgl.bringtotop()
}

.plot.points <- function(shape, shape.functions, radius, col, pca.coords) {

for(s in unique(shape)) {
for(c in unique(col)) {
sel <- which(shape == s & col == c)
if(length(sel) > 0) {
shape.functions[[s]](pca.coords[ sel,, drop=F ],
specular="#111111",
col= as.character(c), alpha= 1, radius= radius, shininess=0 )
}
}
}

}

.labels.3D <- function(show.labels, pca.coords, labels.col, n.p) {
if(is.logical(show.labels) && length(show.labels) == 1) {
show.labels <- rownames(pca.coords)
}
if(is.null(show.labels)) show.labels <- 1:n.p
texts3d(pca.coords[,1:3], texts=show.labels, col=labels.col, adj=c(1, 1))
}

#' @rdname pca3d-package
#' @export
pca3d <- function(pca, components=1:3, col=NULL, title=NULL, new=FALSE,
axes.color ="grey",
bg="white",
group=NULL,
shape=NULL,
palette=NULL,
fancy=FALSE,
biplot=FALSE,
biplot.vars=5,
legend=NULL,
show.scale=FALSE,
show.labels=FALSE,
labels.col="black",
show.axes=TRUE,
show.axe.titles=TRUE,
axe.titles=NULL,
show.plane=TRUE,
show.centroids=FALSE,
show.group.labels=FALSE,
show.shapes=TRUE,
show.ellipses=FALSE,
ellipse.ci=0.95
) {

if(fancy) {
show.labels       <- TRUE
show.centroids    <- TRUE
show.group.labels <- TRUE
}

if(missing(palette))
palette <- defaultPalettePCA3D()

# prepare pca coordinates to plot
pca.coords <- get.pca.coords(pca, 3, components)
n.p <- nrow(pca.coords)
main.ranges <- apply(pca.coords, 2, function(x) c(-max(abs(range(x))), max(abs(range(x)))))

if(!missing(biplot)) {
biplot.coords <- get.biplot.coords(pca, biplot)
biplot.coords <- biplot.coords[ , components]
biplot <- TRUE
biplot.vars <- get.biplot.vars(biplot.coords, biplot.vars)
}

# ----------------------------------
# selecting shapes and colors
if(missing(group))
group <- rep("default", n.p)

group <- factor(group)

col <- autocol(n.p, col, group, palette, "col")
group.col <- getGroupVals(col, group)

shape <- autocol(n.p, shape, group, all.shapes, "shape")
shape <- matchShapes(shape)
group.shape <- getGroupVals(shape, group)

ret <- print.legend(group, group.col, group.shape)
# ----------------------------------

# other things we need for drawing in 3D
shape.functions <- list(sphere= .myspheres3d, tetrahedron= tetrahedrons3d, cube= cubes3d, octahedron=octahedrons3d)
xz.functions    <- list(sphere= .circlexz,    tetrahedron= .trianglexz,     cube= .squarexz, octahedron= .squarexz)

# initialize the plot
.pca_plot_init3d(new,
show.plane, show.axes, show.axe.titles, show.scale,
axes.color, axe.titles, main.ranges, components, bg)

if(show.centroids) .centroids.3D(pca.coords, group, 2 * radius, group.shape, group.col, shape.functions)
if(show.group.labels) .group.labels(pca.coords, group)
if(show.ellipses) .group.ellipses.3D(pca.coords, group, group.col, ellipse.ci)

# draw the samples
if(show.shapes) {
.plot.points(shape, shape.functions, radius, col, pca.coords)
}

# add the biplot
if(biplot)
.biplot.3D(biplot.coords, biplot.vars, main.ranges)

# label the samples
if(!missing(show.labels))
.labels.3D(show.labels, pca.coords, labels.col, n.p)

# finalize the plot
.pca_plot_finish_3d()

# draw legend
if(!is.null(legend) && !is.null(ret)) {
legend3d(legend,
ret\$groups,
col=ret\$colors,
pch=ret\$pch,
box.lty=0
)
}

return(invisible(ret))
}
```

Try the pca3d package in your browser

Any scripts or data that you put into this service are public.

pca3d documentation built on Oct. 23, 2020, 8:26 p.m.