pdf("rpd2_15.pdf")


X1 <- matrix(
c(11.1,25.2,33.0,1576,12.6,26.0,48.2,1767,
9.3,27.3,27.3,1580,9.8,26.7,34.3,1578,
5.8,25.7,27.3,1222,12.1,28.3,20.4,1926),byrow=T,ncol=4)

Xmn <- colMeans(X1)

Xmeans <- rbind(t(Xmn),t(Xmn),t(Xmn),t(Xmn),t(Xmn),t(Xmn))

Xss <- c()

for (i in 1:ncol(X1)) {
Xss[i] <- (nrow(X1)-1)*var(X1[,i])
}

Xsse <- rbind(t(Xss),t(Xss),t(Xss),t(Xss),t(Xss),t(Xss))
X <- (X1-Xmeans)/sqrt(Xsse)

svdX <- svd(X)

U <- svdX$u
Z <- svdX$v
lambdasr <- svdX$d

W <- U%*%diag(lambdasr);

W

site <- seq(1:nrow(X))

plot(-W[,1],W[,2],type="n",xlab='First Principal Component',ylab='Second Principal Component',
xlim=c(-1.2,1.2),ylim=c(-1.2,1.2))
text(-W[,1],W[,2],labels=site)
abline(v=0)
abline(h=0)

# We plot -W[,1] since R gives negative first and last principal component vectors

dev.off()