Example: two dimensions

Suppose we had two variables we observed.

library(ggplot2)
base1 <- rnorm(20)
x1 <- base1 + rnorm(20) * 0.4
y1 <- base1 + rnorm(20) * 0.4

dat <- data.frame(id = 1:20, x1, y1)

ggplot(dat, aes(x = x1, y = y1)) + geom_point(size = 5) + geom_text(aes(label = id), 
    col = "white") + geom_abline(slope = 1, intercept = 0, linetype = 2) + geom_abline(slope = -1, 
    intercept = 0, linetype = 2) + 
xlim(-3.5, 3.5) + ylim(-3.5, 3.5)

The data appear to be correlated, and falling along the line x=y. This is the vector of highest variance, and might be considered the first principal component. Since we have only two dimensions, the only remaining orthogonal dimension is the dimension along the negative diagonal–the second principal component. So, suppose we wanted to rotate the axes and redescribe the data based on these two dimensions. Any point would be re-described by its closest point on the first or second new dimension.

If we look at the the variability along this new first component, it looks like it will have a standard deviation of around 2; and the second dimension will have a standand devitaion much smaller–maybe less than 1. If we do an eigen decomposition on the covariance matrix, we lose the original data (everything is filtered through the covariance matrix), but we have redescribed the data in terms of how they map onto the new vectors.The eignevalues describe the variance

e <- eigen(cov(dat[, -1]))
print(e)

eigen() decomposition
$values
[1] 2.0937353 0.1449465

$vectors
          [,1]       [,2]
[1,] 0.6585969 -0.7524960
[2,] 0.7524960  0.6585969

We can use the eigenvector matrix to rotate the original data into this new set of axes. Notice how the arrangement is the same with an angular rotation; all the points are in the same configuration. Also, the variance in each dimension is the same as the eigenvalues.

rotated <- data.frame(id = dat$id, as.matrix(dat[, 2:3]) %*% e$vectors)

ggplot(rotated, aes(x = X1, y = X2)) + geom_point(size = 5) + geom_text(aes(label = id), 
    col = "white") + geom_abline(intercept = 0, slope = 0)

print(var(rotated$X1))

[1] 2.093735

print(var(rotated$X2))

[1] 0.1449465

So, what eigen decomposition is doing is rotating the entire space to a new set of orthogonal dimensions; so that the first axis is the one with the greatest variance, the next axis is the next-largest, and so on. In this case, because two variables were correlated, it picked the axis of greatest correlation as the first vector. If we have a lot more variables, the first most important dimension might strongly combine only some of the original val`ues, whereas a second dimension might pick out a second set of common variance for a different set of variables.

Example of Eigen Decomposition: Two latent factors

Here is an example where there are two ‘true’ underlying factors that the data are generated from–one set of three items based on one factor, one set based on a second factor, and a final set based on both. There is noise added to each dimension, so we wouldn’t be able to compress the data into just two dimensions, but it might have most of the variance on two.

library(ggplot2)
library(GGally)

base1 <- rnorm(100)
base2 <- rnorm(100)

x1 <- base1 + rnorm(100) * 0.6
x2 <- base1 + rnorm(100) * 0.6
x3 <- base1 + rnorm(100) * 0.6

y1 <- base2 - rnorm(100) * 0.8
y2 <- base2 + rnorm(100) * 0.7
y3 <- base2 - rnorm(100) * 0.7

z1 <- (base1 + base2)/2 + rnorm(100) * 0.5
z2 <- (base1 + base2)/2 + rnorm(100) * 0.5
z3 <- (base1 + base2)/2 + rnorm(100) * 0.5

## z1 <- rnorm(100)*.5 z2 <- rnorm(100)*.5 z3 <- rnorm(100)*.5

dat <- data.frame(x1, x2, x3, y1, y2, y3, z1, z2, z3)

ggpairs(dat)

Within each set of x/y, there is moderately strong correlations; and z vaus are likwise correlated with everything. Rather than treating each variable separately, what if we could make a linear combination of terms–just like regression, that best combine to account for a common pool of variance. But in this case, we have two independent sources of variance. It seems like an impossible task–how can we create a regression predicting a dependent variable we don’t know? And how can we do it when the values we know are some complex unlabeled mixture of these values?

It seems remarkable, but this is what eigen decomposition does. Eigen is a german term meaning its ‘own’, and it is also sometimes referred to as the ‘characteristic vector’. Eigen decomposition looks at the entire multivariate set of data and finds a new set of

e <- eigen(cor(dat))
print(round(e$values, 4))

[1] 4.1691 2.4475 0.4767 0.4139 0.3884 0.3537 0.3031 0.2437 0.2038

round(e$vectors, 4)

         [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
 [1,] -0.3431 -0.3137 -0.2655 -0.2815  0.3351  0.3796  0.5077  0.3462
 [2,] -0.3407 -0.3417 -0.3187 -0.0647 -0.2511  0.1772 -0.5371 -0.0929
 [3,] -0.3440 -0.3501 -0.0608  0.2928  0.1624 -0.1678  0.1839 -0.7171
 [4,] -0.2491  0.4453  0.3949  0.0804 -0.1052  0.5827  0.2094 -0.3110
 [5,] -0.1821  0.5023 -0.4082  0.0763  0.3112 -0.4517  0.1807 -0.0342
 [6,] -0.2500  0.4583 -0.4392 -0.1763 -0.0551  0.2259 -0.3097 -0.0700
 [7,] -0.4037  0.0202  0.4778  0.0418  0.5828 -0.0994 -0.4508  0.2230
 [8,] -0.3991  0.0258  0.2835 -0.6211 -0.4056 -0.4256  0.1593 -0.0354
         [,9]
 [1,] -0.0100
 [2,]  0.5195
 [3,] -0.2592
 [4,]  0.2932
 [5,]  0.4553
 [6,] -0.5904
 [7,] -0.0637
 [8,] -0.0408
 [ reached getOption("max.print") -- omitted 1 row ]

ggplot(data.frame(x = 1:9, values = e$values), aes(x = x, y = values)) + geom_bar(stat = "identity")

Now, the total variance in the first two dimensions is fairly large, and it trails off after that. What about the new dimensions: let’s look at just the first two dimensions:

dat <- data.frame(dim = as.factor(rep(1:9, 9)), vec = as.factor(rep(1:9, each = 9)), 
    val = as.vector(e$vectors))[1:18, ]

ggplot(dat, aes(x = dim, y = val, group = vec, col = vec)) + geom_point() + 
    geom_line()

Here, the first vector captures the shared variance across all questions–the last three questions are related to both groups. Any value close to 0 indicates that the question/item is near 0 on that component. Then, the second component gets coded as a binary factor that either explains items 1-3 or items 4-6. This sort of automatically creates a contrast coding between the first two dimensions.

Example of Eigen Decomposiion: network flows

Suppose you had a network with different nodes and transition probabilities between nodes. For example, you might have three economic classes: poverty, middle-class, and wealthy. Any person from one generation has a certain probability of ending up in each of the other classes, maybe with a structure like this:

transition <- rbind(c(0.5, 0.5, 0), c(0.25, 0.5, 0.25), c(0, 0.5, 0.5))

For any generation of people, there is a distribution of class memebership, and we can derive the subsequent distribution by cross-multiplying

dsn <- runif(3)

dsn <- dsn/sum(dsn)
overtime <- matrix(0, nrow = 100, ncol = 3)
overtime[1, ] <- dsn
for (i in 2:100) {
    dsn <- dsn %*% transition
    # dsn <- dsn/sum(dsn)
    overtime[i, ] <- dsn
}



matplot(overtime, ylim = c(0, 1), type = "o")

e <- eigen(t(transition))


print(e)

eigen() decomposition
$values
[1]  1.000000e+00  5.000000e-01 -7.019029e-17

$vectors
          [,1]          [,2]       [,3]
[1,] 0.4082483 -7.071068e-01  0.4082483
[2,] 0.8164966  8.112147e-16 -0.8164966
[3,] 0.4082483  7.071068e-01  0.4082483

print(dsn)

     [,1] [,2] [,3]
[1,] 0.25  0.5 0.25

Notice that the eigenvectors have one with all positive values, and it is (.4,.8,.4). If we normalize this, we get (.25, .5, .25), which is the stationary distribution of this stochastic matrix. Sometimes, there are multiple stable distributions depending on the starting conditions, and these would each map onto a distinct eigenvector. in this case, the second eigenvector is essentially -1,0,1; this vector cannot happen for a transition matrix (you can’t have negative probability), so in this case there are no other stable distributions.

Here is another case where there is complete ability to transition between groups. Now, the first eigenvector is all negative, which also predicts the stationary distribution, at least ordinally.

transition <- rbind(c(0.5, 0.4, 0.1), c(0.4, 0.5, 0.1), c(0.1, 0.4, 0.5))
e <- eigen(t(transition))
print("Eigen decomposition:")

[1] "Eigen decomposition:"

print(e)

eigen() decomposition
$values
[1] 1.0 0.4 0.1

$vectors
           [,1]          [,2]          [,3]
[1,] -0.6337502 -7.071068e-01 -7.071068e-01
[2,] -0.7242860  2.266233e-16  7.071068e-01
[3,] -0.2716072  7.071068e-01  3.462330e-17

e$vectors[, 1]/sum(e$vectors[, 1])

[1] 0.3888889 0.4444444 0.1666667

dsn <- runif(3)

dsn <- dsn/sum(dsn)
overtime <- matrix(0, nrow = 100, ncol = 3)
overtime[1, ] <- dsn
for (i in 2:100) {
    dsn <- dsn %*% transition
    # dsn <- dsn/sum(dsn)
    overtime[i, ] <- dsn
}


print("----------------------")

[1] "----------------------"

matplot(overtime, ylim = c(0, 1))

print("Stationary distribution:")

[1] "Stationary distribution:"

print(dsn)

          [,1]      [,2]      [,3]
[1,] 0.3888889 0.4444444 0.1666667

plot(dsn, e$vectors[, 1], main = "Comparing  stable distribution to first eigenvector")

Factors and Dimensions

Oftentimes, the number of things we can measure about a person are greater than the number of theoretically-interesting constructs we think these measureables explain. For example, on a personality questionairre, one asks a bunch of related questions that we believe all involve the same thing, with the hope that we will get robust and reliable measures. A person might by chance agree with one of the statements, but it is unlikely that they will agree with them all just by chance.

Measurements and questions that hang together like this are conceived of as a factor. This is similar to the previous notions of clusters, but there can be differences. A factor is a set of questions that covary together; not one in which the same answers are given by everybody. If you know all but one answer someone gives to the answers that are in a factor, you should have a pretty good idea of what the remaining answer is, regardless of whether their answer is low or high on the factor.

Thus, if you have a 10 measurements (answers) that are all a part of a common factor, antd if you can project them down onto a single combined dimension without losing much information, you have a succeeded. For example, if you were to measure left arm length, right arm length, left thigh length, right thigh length, left forearm length, right forearm length, left calf length, and right calf length, chances are these will all vary fairly closely together. A factor describing these might be a weighted average of these 8 measures (maybe even weighing them all equally). These would project down from 8 dimensions to just 1 or 2. If you started adding other measurements (measurements of shoe size, weight, waist size, etc.), these may not be a closely related, and may in fact form other factors.

A number of related methods have been developed to do this sort of analysis, and they go generally under the name ‘factor analysis’, but some are also called Principal Components Analysis, exploratory factor analysis (EFA), Confirmatory Factor Analysis (CFA), and as you make more complex hypotheses about the relationship between variable, there are other related methods (Structural Equation Modeling, LISREL, Path Analysis, etc.). The naming of these is often confusing and historically has been tied to specific software packages rather than generic analysis routines, so it is important to understand clearly what methods you are using, and what methods other are using, when you do this sort of analysis.

Typically, factor analysis focuses on analyzing the correlational structure of a set of questions or measures. When thinking about correlations, a factor might be a set of questions that are highly positively correlated within the set, but are uncorrelated with questions outside the set. Typically, negatively-correlated items are incorporated into the factor the same way positively-correlated items would be. This is somewhat like looking int the distance measures we examined with previous methods. With N questions, there are N^2 correlations in a matrix, which is sort of like placing each question is an N-dimensional space. Factor analysis attempts to find a number of dimensions smaller than N that explains most of the same data.

Example: four questions with two factors

Let’s suppose we have a set of 4 questions related to two different factors: happiness and intelligence. for example:

data <- rbind(c(5, 5.5, 100, 3), c(4.5, 5, 110, 3.5), c(3, 3.5, 130, 4), c(1, 
    1, 140, 3.9), c(0, 1, 80, 2.5), c(1, 2, 90, 3.1), c(6, 6, 70, 2.3), c(2, 
    2.5, 60, 1.9))

# pairs(data) cor(data)
library(GGally)
ggpairs(as.data.frame(data))

Notice that Q1/Q2 are highly correlated, ase are Q3/Q4, but the there is low correlations between these clusters. A factor analysis should be able to figure it out.

Principal Components Analysis/eigen decomposition as factor analysis

Principal Components Analysis (PCA) is the simplest version of factor analysis. It is typically based on eigen decomposition, which produces an outcome somewhat like an ANOVA model does. That is, for an ANOVA model, you might try to predict responses based on a weighted average of predictor variables. The outcome of an eigen decomposition of the correlation matrix finds a weighted average of predictor variables that can reproduce the correlation matrix…without having the predictor variables to start with.

What PCA does is transforms the data onto a new set of axes that best account for common data.

e <- eigen(cor(data))
plot(e$values)

e$vectors

           [,1]       [,2]       [,3]       [,4]
[1,]  0.4968782 -0.5053846 -0.3846191  0.5914107
[2,]  0.5179202 -0.4791360  0.3608988 -0.6098684
[3,] -0.5020219 -0.4959033 -0.5930808 -0.3876973
[4,] -0.4825399 -0.5187437  0.6083382  0.3577496

barplot(t(e$vectors[, 1:2]), beside = T)

# 

If we look at e\$values, this gives a measure of the importance of each dimension. e$vectors gives us the actual dimensions. Each column of e$vectors is considered a factor–a way of weighting answers from each question to compute a composite. Here, after 2 dimensions, the eigen values go below 1.0, and fall off considerably, indicating most of the variability is accounted for. Let’s look at another slightly more complex data set: one on human abilities. We will use ability.cov–the covariance of these abilitiies.

Example: Human abilities

library(GPArotation)
data(ability.cov)
e2 <- eigen(cov2cor(ability.cov$cov))
plot(e2$values)

evec <- e2$vectors
rownames(evec) <- rownames(ability.cov$cov)
round(evec, 3)

          [,1]   [,2]   [,3]   [,4]   [,5]   [,6]
general -0.471  0.002  0.072  0.863  0.037 -0.164
picture -0.358  0.408  0.593 -0.268  0.531  0.002
blocks  -0.434  0.404  0.064 -0.201 -0.775  0.051
maze    -0.288  0.404 -0.794 -0.097  0.334  0.052
reading -0.440 -0.507 -0.015 -0.100  0.056  0.733
vocab   -0.430 -0.501 -0.090 -0.352  0.021 -0.657

Here, there is a large drop-off between 1 and 2 eigenvalues, but it really flattens out after 3 dimensions. Looking at the eigenvectors, vector values that are different from 0 indicate a high positive or negative involvement of a particular question in a particular factor. Here, the first factor is generally equally involved in each question–maybe general intelligence. This makes sense if you look at the covariance, because there is a lot of positive correlation among the whole set. The second factor is Q234 - Q56; maybe visual versus lexical intelligence. This is one way of carving up the fact that there might be a general intelligence that impacts everyone, and a second factor that isn’t general but might be visual or vocabulary. This is a bit misleading though–if someone is good at both or bad at both, it gets explained as general intelligence; if they are good or bad at one or the other, it gets explained as a specific type of intelligence.

The third vector is essentially picture versus maze, which load together on factor 2. Maybe there are aspects of these two visual tasks that differ. Looking at this, a hierarchical structure appears, which might indicate a hierarchical clustering would be better. I’ll create a distance measure from the covariance converting to a correlation and subtracting from 1.0.

library(cluster)

cor <- ability.cov$cov/sqrt(diag(ability.cov$cov) %*% t(diag(ability.cov$cov)))
dist <- as.dist(1 - cor)
a <- agnes(dist)
plot(a, which = 2)

This sort of maps onto our intuition; maze here looks like a bit of an outlier.

The princomp function does the same thing as eigen(), with a few additional bells and whistles. Notice that the vectors are identical. Another option in the psych::principal function.

p <- princomp(covmat = cov2cor(ability.cov$cov))
plot(p)
summary(p)

Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
Standard deviation     1.7540877 1.0675615 0.9039839 0.64133665 0.59588119
Proportion of Variance 0.5128039 0.1899479 0.1361978 0.06855212 0.05917907
Cumulative Proportion  0.5128039 0.7027518 0.8389497 0.90750178 0.96668085
                           Comp.6
Standard deviation     0.44711845
Proportion of Variance 0.03331915
Cumulative Proportion  1.00000000

loadings(p)

Loadings:
        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
general  0.471                0.863         0.164
picture  0.358  0.408  0.593 -0.268  0.531
blocks   0.434  0.404        -0.201 -0.775
maze     0.288  0.404 -0.794         0.334
reading  0.440 -0.507        -0.100        -0.733
vocab    0.430 -0.501        -0.352         0.657

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
SS loadings     1.000  1.000  1.000  1.000  1.000  1.000
Proportion Var  0.167  0.167  0.167  0.167  0.167  0.167
Cumulative Var  0.167  0.333  0.500  0.667  0.833  1.000

screeplot(p)

The advantage of princomp is that it lets you use the raw data, and will compute scores based on those base data. There is a related function called prcomp, which will do a similar analysis, but handles situations where you have fewer participants than questions.

data <- read.csv("bigfive.csv")

dat.vals <- data[, -1]  ##remove subject code
qtype <- c("E", "A", "C", "N", "O", "E", "A", "C", "N", "O", "E", "A", "C", 
    "N", "O", "E", "A", "C", "N", "O", "E", "A", "C", "N", "O", "E", "A", "C", 
    "N", "O", "E", "A", "C", "N", "O", "E", "A", "C", "N", "O", "O", "A", "C", 
    "O")
valence <- c(1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, 1, 
    1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, -1, 
    1, -1, 1)
add <- c(6, 0, 0)[valence + 2]


tmp <- dat.vals
tmp[is.na(tmp)] <- 3.5
reversed <- t(t(tmp) * valence + add)
## reverse code questions:

bytype <- reversed[, order(qtype)]
key <- sort(qtype)
colnames(bytype) <- paste(key, 1:44, sep = "")

pc <- princomp(bytype)
plot(pc)

cbt <- cor(bytype)
image(cbt, col = grey(100:0/100))

length <- length(key)
noise <- 0.1
f.a <- (key == "A") + runif(length(key)) * noise
f.c <- (key == "C") + runif(length(key)) * noise
f.e <- (key == "E") + runif(length(key)) * noise
f.n <- (key == "N") + runif(length(key)) * noise
f.o <- (key == "O") + runif(length(key)) * noise

cor.a <- f.a %*% t(f.a)
cor.c <- f.c %*% t(f.c)
cor.e <- f.e %*% t(f.e)
cor.n <- f.n %*% t(f.n)
cor.o <- f.o %*% t(f.o)
image(cor.a, col = grey(100:0/100))

image(cor.c, col = grey(100:0/100))

image(cor.e, col = grey(100:0/100))

image(cor.n, col = grey(100:0/100))

image(cor.o, col = grey(100:0/100))

par(mfrow = c(1, 2))
image(cor.a + cor.c + cor.e + cor.n + cor.o, col = grey(100:0/100))
image(cbt, col = grey(100:0/100))

e <- eigen(cbt)
par(mfrow = c(2, 3))
image(cbt, col = grey(100:0/100), main = "Complete correlation")
image((e$vectors[, 1] * t(e$values[1])) %*% t(e$vectors[, 1]), col = grey(100:0/100), 
    main = "One eigenvector")
image((e$vectors[, 1:2] * (e$values[1:2])) %*% t(e$vectors[, 1:2]), col = grey(100:0/100), 
    main = "Two eigenvectors")
image((e$vectors[, 1:3] * (e$values[1:3])) %*% t(e$vectors[, 1:3]), col = grey(100:0/100), 
    main = "Three eigenvectors")
image((e$vectors[, 1:4] * (e$values[1:4])) %*% t(e$vectors[, 1:4]), col = grey(100:0/100), 
    main = "Four eigenvectors")
image((e$vectors[, 1:5] * (e$values[1:5])) %*% t(e$vectors[, 1:5]), col = grey(100:0/100), 
    main = "Five eigenvectors")

par(mfrow = c(1, 2))
image(cbt, col = grey(100:0/100), main = "Complete correlation")
image((e$vectors[, 1:10] * e$values[1:10]) %*% t(e$vectors[, 1:10]), col = grey(100:0/100), 
    main = "Ten eigenvectors")

factors1 <- eigen(cov(bytype))

## a 15-factor solution using SVD
factors2 <- svd(bytype, nu = 15)
summary(factors1)

        Length Class  Mode
values    44   -none- numeric
vectors 1936   -none- numeric

summary(factors2)

  Length Class  Mode
d    44  -none- numeric
u 15255  -none- numeric
v  1936  -none- numeric

par(mfrow = c(1, 3))
plot(factors1$values, main = "Eigen scree plot")
plot(factors2$d, main = "SVD scree plot")
plot(factors2$d[-1], main = "SVD scree plot removing first dimension")

barplot(factors2$v[, 1])

par(mfcol = c(4, 2))
barplot(factors1$vectors[, 1], main = "Eigen D1")
barplot(factors1$vectors[, 2], main = "Eigen D2")
barplot(factors1$vectors[, 3], main = "Eigen D3")
barplot(factors1$vectors[, 4], main = "Eigen D4")


barplot(factors2$v[, 2], main = "SVD D2")
barplot(factors2$v[, 3], main = "SVD D3")
barplot(factors2$v[, 4], main = "SVD D4")
barplot(factors2$v[, 5], main = "SVD D5")

round(cor(factors1$vectors, factors2$v), 3)[1:10, 1:10]

        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,]  0.339  0.985  0.018  0.120 -0.033 -0.050 -0.050 -0.044 -0.034
 [2,] -0.116 -0.312  0.875  0.330 -0.026 -0.061 -0.087 -0.083 -0.044
 [3,]  0.441 -0.345 -0.685  0.593 -0.021 -0.074 -0.124 -0.122 -0.055
 [4,]  0.341  0.188  0.117  0.865 -0.266  0.124  0.209  0.192  0.070
 [5,] -0.408  0.054  0.027  0.203  0.973  0.030  0.059  0.056  0.016
 [6,]  0.511  0.012  0.013  0.027  0.015 -0.998  0.041  0.031  0.013
 [7,]  0.147  0.033  0.019  0.047  0.016 -0.004 -0.950  0.286  0.042
       [,10]
 [1,]  0.032
 [2,]  0.052
 [3,]  0.072
 [4,] -0.099
 [5,] -0.027
 [6,] -0.014
 [7,] -0.063
 [ reached getOption("max.print") -- omitted 3 rows ]