rm(list = ls()) # clear the memory ########################### X = as.matrix(read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T8-4.DAT")) colnames(X) = c("JPM", "Citi", "WellsF", "Shell", "Exxon") out = prcomp(X, scale = T) # out\$center = colMeans(X) out\$center # out\$scale = vector of standard deviations of the columns of X out\$scale # out\$sdev = square roots of the eigenvalues of R (assuming that scale = T); # thus sum(out\$sdev^2) = p out\$sdev Gamma = out\$rotation # matrix whose columns are the eigenvectors of R scores = out\$x # columns are the sample principal components # cov(scores) = diag(out\$sdev^2) = diag matrix of the eigenvalues of R # check: compare round(cov(scores),4) with eigen(cor(X)) plot(out, type='l', main = "Scree plot") # A 'scree' plot indicating the contributions of the sample pc's # The first two account for 3/4 of the variation; the first three for 87%. cumsum(out\$sdev^2)/5 # [1] 0.4874546 0.7688572 0.8689597 0.9489660 1.0000000 lam1 = out\$sdev[1]^2 ci = c(log(lam1) - sqrt(2)*qnorm(.975)/sqrt(nrow(X)), log(lam1) + sqrt(2)*qnorm(.975)/sqrt(nrow(X))) ci1 = exp(ci) lam2 = out\$sdev[2]^2 ci = c(log(lam2) - sqrt(2)*qnorm(.975)/sqrt(nrow(X)), log(lam2) + sqrt(2)*qnorm(.975)/sqrt(nrow(X))) ci2 = exp(ci) ############################################################# rm(list = ls()) # clear the memory X = as.matrix(read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T5-8.DAT")) colnames(X) = c("legal", "extraordinary", "holdover", "COA", "meeting") #X0 = X-as.vector(rep(1,16))%*%t(as.vector(colMeans(X))) # R will always center the data at least, so the result using X0 will be the same as using X # thus the sample scores will always have an average of 0 ## Without scaling: par(mfrow = c(3,1)) out = prcomp(X, scale = F) plot(out, type = 'l', main = "") title(main = "pca; without scaling", outer = T, line = -2) round(cumsum(out\$sdev^2)/sum(out\$sdev^2),3) out ## Plot first two pca's and confidence ellipsoid scores = out\$x p1 = scores[,1] p2 = scores[,2] r = sqrt(qchisq(.95,2)) phi = 2*pi*seq(0, 1, length=200) z = rbind(r*cos(phi), r*sin(phi)) plot(out\$sdev[1]*z[1,],out\$sdev[2]*z[2,], type = 'l', ylim = c(-4000, 4000), xlab = "first pc", ylab = "second pc") points(p1, p2) text(p1[11], p2[11], "#11", adj = 1.5) y345 = scores[,3:5] lam345 = diag(out\$sdev^2)[3:5,3:5] t = diag(y345%*%solve(lam345,t(y345))) plot(t, type = 'l', ylim=c(0,8.2), xlab = "pay period") abline(h = qchisq(.95, 3)) text(12, t[12], "#12", adj = 1.5) text(13, t[13], "#13", adj = -1) ## How much of the variation of the individual variables is explained by the correlatons with the first two pcs? RHO = diag(1/sqrt(diag(cov(X))))%*%out\$rotation%*%diag(out\$sdev) RHO^2 # check: rowSums(RHO^2) #################################### ## With scaling (but then does the asymptotic chi-square approximation remain valid?): out = prcomp(X, scale = T) dev.new() par(mfrow = c(3,1)) plot(out, type = 'l', main = "") title(main = "pca; with scaling", outer = T, line = -2) round(cumsum(out\$sdev^2)/sum(out\$sdev^2),3) out ## Plot first two pca's and confidence ellipsoid scores = out\$x p1 = scores[,1] p2 = scores[,2] r = sqrt(qchisq(.95,2)) phi = 2*pi*seq(0, 1, length=200) z = rbind(r*cos(phi), r*sin(phi)) plot(out\$sdev[1]*z[1,],out\$sdev[2]*z[2,], xlab = "first pc", ylab = "second pc", ylim = c(-3,3), type = 'l') points(p1, p2) text(p1[11], p2[11], "#11", adj = 1.5) y345 = scores[,3:5] lam345 = diag(out\$sdev^2)[3:5,3:5] t = diag(y345%*%solve(lam345,t(y345))) plot(t, type = 'l', ylim=c(0,8.2), xlab = "pay period") abline(h = qchisq(.95, 3)) text(12, t[12], "#12", adj = 1.5) text(13, t[13], "#13", adj = -1) ## How much of the variation of the individual variables is explained by the correlatons with the first two pcs? RHO = out\$rotation%*%diag(out\$sdev) # sigma[i,i] = 1 when prcomp based on R RHO^2 # check: rowSums(RHO^2)