pca_scratch.r

####
# Eigens
#####

# How to calculate covariance matrix
# Great video: https://www.youtube.com/watch?v=9B5vEVjH2Pk
dat <- as.matrix(
  cbind(c(90, 90, 60, 30, 30),
        c(80, 60, 50, 40, 20),
        c(40, 80, 70, 70, 70))
)

apply(dat, MARGIN = 2, sum) # sum of scores for each subject
apply(dat, MARGIN = 2, mean) # mean score for each subject

# to calculate covariance matrix, we need to do
# A*At
#
# First, let's calculate the deviation matrix
#
# a = A - [1]*A*(1/5)
# In fact, A represents scores, and ( [1]*A*(1/5) ) represents the average
# scores, since [1]*A gives us the sum of each column in A, and dividing it
# by 1/5 gives us the average for each

# This is the unity matrix
unity <- matrix(1, 5, 5) # 5 x 5 matrix
unity %*% dat / dim(unity)[1] # 5x5 matrix times dat, divided by dat column lengths of 5
a = dat - unity %*% dat / dim(unity)[1] # deviation matrix (scores - avg score)

# these give us different results, one a student covariance matrix,
# and one a subject covariance matrix
a %*% t(a) # student cov
t(a) %*% a # subject cov

subjects = t(a) %*% a
eobj = eigen(subjects)

# eigenvalue % contributions
eobj$values / sum(eobj$values)
eobj$vectors

# to calculate loadings
eobj$vectors %*% sqrt(eobj$values) #referred to as P. P%*%t(P) = correlation matrix

# In fact, X = P L P'
eobj$vectors %*% diag(eobj$values) %*% t(eobj$vectors)
# This equals t(a) %*% a

# calculating scores
# first, "center" and normalize by sd. Then multiply by eigenvectors
apply((((dat - mean(dat)) / sd(dat)) %*% eobj$vectors), 2, var)
# the above should be the same as the eigenvalues


###

svdobj <- svd(a)

svdobj$d*svdobj$d / sum(svdobj$d*svdobj$d) # same as eigenvalues
t(svdobj$v) #identical to eigenvectors

# eigen(a %*% t(a))$vectors is the same as svdobj$u