# using R for a principal component analysis

help(prcomp)

# review chromatography data from last class

load("Chrom3D.RData")

View(chrom3d)

library(plot3D)

persp3D(x = time, y = wavelength, z = as.matrix(chrom3d), xlab = "time", ylab = "wavelength", zlab = "absorbance", ticktype = "detailed", phi = 0, theta = 45)

contour2D(z = as.matrix(chrom3d), x = time, y = wavelength, xlab = "time", ylab = "wavelength")

plot(x = time, y = chrom3d$'282', type = "l", lwd = 2, col = "blue", xlab = "time", ylab = "absorbance at 282 nm")

plot(x = wavelength, y = chrom3d[8, ], type = "l", lwd = 2, col = "blue", xlab = "wavelength", ylab = "absorbance at 8 minutes")

# run a PCA on chromatography data and examine output

chrom.pca = prcomp(chrom3d)

summary(chrom.pca)

head(chrom.pca$x)

head(chrom.pca$rotation)

plot(chrom.pca$x, pch = 19, col = "blue", type = "b")

identify(chrom.pca$x, labels = time)

plot(time, chrom.pca$x[ , 1], type = "l", col = "blue", lwd = 2)

plot(time, chrom.pca$x[ , 2], type = "l", col = "blue", lwd = 2)

plot(chrom.pca$rotation, pch = 19, col = "blue", type = "b")

identify(chrom.pca$rotation, labels = wavelength)

plot(wavelength, chrom.pca$rotation[ , 1], type = "l", col = "blue", lwd = 2)

plot(wavelength, chrom.pca$rotation[ , 2], type = "l", col = "blue", lwd = 2)

# centering and scaling

test = c(22, 21, 24, 23, 25)

mean(test)

var(test)

center = test - mean(test)

center

mean(center)

var(center)

scale = test/sd(test)

scale

mean(scale)

var(scale)

c.and.s = (test - mean(test))/sd(test)

c.and.s

mean(c.and.s)

var(c.and.s)

# using PCA for pattern recognition

load("Elements.RData")

View(elements)

elem.pca = prcomp(elements[ , 3:7], center = TRUE, scale = TRUE)

summary(elem.pca)

old.par = par(mfrow = c(1, 2))

plot(elem.pca$x, pch = 19, col = elements$group)

identify(elem.pca$x, labels = elements$element)

plot(elem.pca$rotation, pch = 19, col = "blue")

elem.pca$rotation

par(old.par)

elem.pca = prcomp(elements[, 3:6], center = TRUE, scale = TRUE)

summary(elem.pca)

old.par = par(mfrow = c(1, 2))

plot(elem.pca$x, pch = 19, col = elements$group)

identify(elem.pca$x, labels = elements$element)

plot(elem.pca$rotation, pch = 19, col = "blue")

elem.pca$rotation

plot(elem.pca$x[ , 1], elem.pca$x[ , 3], pch = 19, col = elements$group)

identify(elem.pca$x, labels = elements$element)

plot(elem.pca$rotation[ , 1], elem.pca$rotation[ , 3], pch = 19, col = "blue")

elem.pca$rotation

par(old.par)

scatter3D(x = elem.pca$x[ , 1], y = elem.pca$x[ , 2], z = elem.pca$x[ , 3], pch = 19, col = elements$group, phi = 15, theta = 45, colvar = NULL, ticktype = "detailed", xlab = "pc1", ylab = "pc2", zlab = "pc3")