demo/princomp.R

In yihui/MSG: Data and Functions for the Book Modern Statistical Graphics

library(grid)
library(mvtnorm)
set.seed(123)
sigma = matrix(c(5, 2, 2, 3), ncol = 2)
x = rmvnorm(n = 200, mean = c(1, 2), sigma = sigma)
range.v1 = range(x[, 1])
range.v2 = range(x[, 2])
ext.r = max(diff(range.v1), diff(range.v2)) / 2
scale.v1 = mean(range.v1) + c(-1.05, 1.05) * ext.r
scale.v2 = mean(range.v2) + c(-1.05, 1.05) * ext.r
vp = viewport(
  x = 0.58, y = 0.58, width = 0.8,
  height = 0.8,
  xscale = scale.v1,
  yscale = scale.v2
)
grid.newpage()
pushViewport(vp)
grid.rect()
grid.xaxis()
grid.yaxis()
grid.text("V1", y = unit(-3, "lines"))
grid.text("V2", x = unit(-3, "lines"), rot = 90)
grid.points(x[, 1], x[, 2],
  pch = 19,
  gp = gpar(col = "blue", alpha = 0.3)
)
grid.points(mean(x[, 1]), mean(x[, 2]),
  pch = 19,
  gp = gpar(col = "red")
)
grid.clip()
pr = princomp(x)
ran.pr = apply(pr$scores, 2, function(x) max(abs(x)))
scale.pr = rbind(-1.5 * ran.pr, 1.5 * ran.pr)
size.pr = 3 * ran.pr / c(diff(scale.v1), diff(scale.v2))
angle.pr = asin(pr$loadings[1, 2]) / pi * 180
vp.xaxis = viewport(
  x = unit(mean(x[, 1]), "native"),
  y = unit(mean(x[, 2]), "native"),
  width = size.pr[1],
  just = c("center", "bottom"),
  angle = angle.pr,
  xscale = scale.pr[, 1],
  gp = gpar(col = "red")
)
pushViewport(vp.xaxis)
grid.lines(c(-2, 2), c(0, 0))
grid.xaxis()
popViewport()
vp.yaxis = viewport(
  x = unit(mean(x[, 1]), "native"),
  y = unit(mean(x[, 2]), "native"),
  height = size.pr[2],
  just = c("left", "center"),
  angle = angle.pr,
  yscale = scale.pr[, 2],
  gp = gpar(col = "red")
)
pushViewport(vp.yaxis)
grid.lines(c(0, 0), c(-2, 2))
grid.yaxis()