### R code from vignette source 'chap9.Rnw'

###################################################
### code chunk number 1: chap9.Rnw:5-6
###################################################
options(show.signif.stars=FALSE,digits=4,width=80)


###################################################
### code chunk number 2: spm
###################################################
data(banknote, package="alr3")
colMeans(banknote)
print(cor(banknote)[-7, -7], digits=3)
library(car)
print(scatterplotMatrix(~ . - Y | Y, data=banknote, smooth=FALSE, reg.line=FALSE,
   ellipse=TRUE, levels= c(0.9), by.groups=FALSE, diagonal="none"))


###################################################
### code chunk number 3: chap9.Rnw:26-27
###################################################
print(p1 <- prcomp(~ . - Y, data = banknote, scale=TRUE), digits=3)


###################################################
### code chunk number 4: chap9.Rnw:33-35
###################################################
n1 <- dim(banknote)[1] - 1
s <- svd(sqrt(1/n1) * scale(banknote[, -7]))


###################################################
### code chunk number 5: chap9.Rnw:38-39
###################################################
s$d


###################################################
### code chunk number 6: chap9.Rnw:45-46
###################################################
print(head(p1$x),digits=3)


###################################################
### code chunk number 7: chap9.Rnw:49-50
###################################################
summary(p1)


###################################################
### code chunk number 8: zero
###################################################
plot(p1)


###################################################
### code chunk number 9: one
###################################################
scatterplotMatrix(~PC1 + PC2 + PC3|factor(banknote$Y, labels=c("genuine", "fake")),
     data=p1$x,  diagonal="none", smooth=FALSE, reg.line=FALSE)


###################################################
### code chunk number 10: chap9.Rnw:71-72
###################################################
names(p1)


###################################################
### code chunk number 11: chap9.Rnw:77-79
###################################################
n <- dim(p1$x)[1]
(int <- log(p1$sdev[1]^2) + c(-1.96, 1.96)  * sqrt(2/(n-1)))


###################################################
### code chunk number 12: chap9.Rnw:82-83
###################################################
exp(int)


###################################################
### code chunk number 13: chap9.Rnw:92-99
###################################################
evals <- p1$sdev^2 # for convenience
psihat <- sum(evals[1:3]) / sum(evals)
betahat <- sum(evals[1:3]^2) / sum(evals^2)
trS <- sum(evals)
trS2 <- sum(evals^2)
omegahat2 <- (2 * trS2/(trS^2)) * (psihat^2 - 2*betahat * psihat + betahat^2)
psihat + c(-1, 1) * sqrt(omegahat2)/(sqrt(n-1))


###################################################
### code chunk number 14: chap9.Rnw:113-114
###################################################
p2 <- update(p1, scale=FALSE)


###################################################
### code chunk number 15: chap9.Rnw:117-121
###################################################
rXY <- diag(1/apply(banknote[,-7],2, sd)) %*% p2$rotation %*%
          diag(p2$sdev)
rownames(rXY) <- colnames(banknote)[-7]
print(rXY, digits=3)


###################################################
### code chunk number 16: chap9.Rnw:125-126
###################################################
(rXY <- with(p1, rotation %*% diag( sdev)))


###################################################
### code chunk number 17: pcfig1
###################################################
plot(rXY[,1:2], asp=1, xlim=c(-1, 1), ylim=c(-1, 1), pch=as.character(1:6))
abline(h=0, lty=2)
abline(v=0, lty=2)
require(plotrix)
draw.circle(0, 0, 1, lty=2, lwd=2)
draw.circle(0, 0, 0.75, lty=2)


###################################################
### code chunk number 18: onea
###################################################
loc <- "http://www.stat.umn.edu/~sandy/courses/8053/Data/uscomp1.dat"
head(uscomp <- read.table(url(loc),header=TRUE))
scatterplotMatrix(uscomp, id.n=3, smooth=FALSE, reg=FALSE)


###################################################
### code chunk number 19: chap9.Rnw:161-163
###################################################
print(pr1 <- prcomp(uscomp, center=TRUE, scale=TRUE), digits=3)
summary(pr1)


###################################################
### code chunk number 20: two
###################################################
y <- as.data.frame(pr1$x)
scatterplot(y$PC1, y$PC2, smooth=FALSE, reg=FALSE, box=F, id.n=4)


###################################################
### code chunk number 21: three
###################################################
print(pr2 <- prcomp(uscomp[-c(38, 40), ],center=TRUE, scale=TRUE), digits=3)
summary(pr2)
cor(pr2$x[, 1:2],uscomp[-c(38, 40), ])
par(mfrow=c(1, 2))
plot(pr2$x[, 1],pr2$x[, 2], main="2 PCs without 38, 40")
plot(pr2$x[, 1],pr1$x[-c(38, 40), 1], main="Comparison of first PCs")


###################################################
### code chunk number 22: four
###################################################
#pairs(uscomp)
library(car)
summary(powerTransform(uscomp,family="yjPower"))
print(pr3<-prcomp(~log(Assets) + log(Sales) + log(MarketValue) + Profits +
           CashFlow + log(Employees), data=uscomp, scale=TRUE),
           digits=3)
summary(pr3)
scatterplot(pr3$x[, 1], pr3$x[, 3], box=FALSE, reg=FALSE,
            smooth=FALSE, id.n=4)


###################################################
### code chunk number 23: five
###################################################
round(rXY <- with(pr3, rotation %*% diag( sdev)), 2)
plot(rXY[,1:2], asp=1, xlim=c(-1, 1), ylim=c(-1, 1), pch=as.character(1:6))
abline(h=0, lty=2)
abline(v=0, lty=2)
require(plotrix)
draw.circle(0, 0, 1, lty=2, lwd=2)
draw.circle(0, 0, .75, lty=2, lwd=1)