EFA and PCA

Thinking about measurement

In the second block of the Multivariate Stats course, we’re going to look at “latent variable modeling.”

Take a moment to think about the various constructs that you are interested in as a researcher. This might be anything from personality traits, to language proficiency, social identity, anxiety etc. How we measure such constructs is a very important consideration for research. The things we’re interested in are very rarely the things we are directly measuring.

The most common way to start thinking about latent variables is to consider how we might assess levels of anxiety or depression. More often than not, we measure these using questionnaire based methods. Can we ever directly measure anxiety?.¹ This inherently leads to the notion of measurement error - two people both scoring 11 on the Generalised Anxiety Disorder 7 (GAD-7) scale will not have identical levels of anxiety.

Before we start to explore how we might use a latent variable modelling approach to accommodate this idea, we’re going to start with the basic concept of reducing the dimensionality of our data.

New packages!

We’re going to be needing some different packages this week, no more lme4!
Make sure you have these packages installed.

• psych
  
• GPArotation
  
• car

Principal component analysis (PCA)

The goal of principal component analysis (PCA) is to find a smaller number of uncorrelated variables which are linear combinations of the original (many) variables and explain most of the variation in the data.

Job Performance: Data Codebook

The file job_performance.csv (available at https://uoepsy.github.io/data/job_performance.csv) contains data on fifty police officers who were rated in six different categories as part of an HR procedure. The rated skills were:

• communication skills: commun
• problem solving: probl_solv
• logical ability: logical
• learning ability: learn
• physical ability: physical
• appearance: appearance

Question A1

Load the job performance data into R and all it job. Check whether or not the data were read correctly into R - do the dimensions correspond to the description of the data above?

Solution

Let’s load the data:

library(tidyverse)

job <- read_csv('https://uoepsy.github.io/data/job_performance.csv')
dim(job)

## [1] 50  6

There are 50 observations on 6 variables.

The top 6 rows in the data are:

head(job)

## # A tibble: 6 x 6
##   commun probl_solv logical learn physical appearance
##    <dbl>      <dbl>   <dbl> <dbl>    <dbl>      <dbl>
## 1     12         52      20    44       48         16
## 2     12         57      25    45       50         16
## 3     12         54      21    45       50         16
## 4     13         52      21    46       51         17
## 5     14         54      24    46       51         17
## 6     14         48      20    47       51         18

Question A2

Provide descriptive statistics for each variable in the dataset.

Solution

We now inspect some descriptive statistics for each variable in the dataset:

# Quick summary
summary(job)

##      commun        probl_solv       logical          learn      
##  Min.   :12.00   Min.   :48.00   Min.   :20.00   Min.   :44.00  
##  1st Qu.:16.00   1st Qu.:52.25   1st Qu.:22.00   1st Qu.:48.00  
##  Median :18.00   Median :54.00   Median :24.00   Median :50.00  
##  Mean   :17.68   Mean   :54.16   Mean   :24.02   Mean   :50.28  
##  3rd Qu.:19.75   3rd Qu.:56.00   3rd Qu.:26.00   3rd Qu.:52.00  
##  Max.   :24.00   Max.   :59.00   Max.   :31.00   Max.   :56.00  
##     physical       appearance   
##  Min.   :48.00   Min.   :16.00  
##  1st Qu.:52.25   1st Qu.:19.00  
##  Median :54.00   Median :21.00  
##  Mean   :54.16   Mean   :21.06  
##  3rd Qu.:56.00   3rd Qu.:23.00  
##  Max.   :59.00   Max.   :28.00

OPTIONAL

If you wish to create a nice looking table for a report, you could try the following code. However, I should warn you: this code is quite difficult to understand so, if you are interested, attend a drop-in!

library(gt)

# Mean and SD of each variable
job %>% 
  summarise(across(everything(), list(M = mean, SD = sd))) %>%
  pivot_longer(everything()) %>% 
  mutate(
    value = round(value, 2),
    name = str_replace(name, '_M', '.M'),
    name = str_replace(name, '_SD', '.SD')
  ) %>%
  separate(name, into = c('variable', 'summary'), sep = '\\.') %>%
  pivot_wider(names_from = summary, values_from = value) %>% 
  gt()

variable	M	SD
commun	17.68	2.74
probl_solv	54.16	2.41
logical	24.02	2.49
learn	50.28	2.84
physical	54.16	2.41
appearance	21.06	2.99

Is PCA needed?

If the original variables are highly correlated, it is possible to reduce the dimensionality of the problem under investigation without losing too much information.

On the other side, when the correlation between the variables under study is weak, a larger number of components is needed in order to explain sufficient variability.

Question A3

Investigate whether or not the recorded variables are highly correlated and explain whether or not you PCA might be useful in this case.

Hint: We only have 6 variables here, but if we had many, how might you visualise cor(job)?

Solution

Let’s start by looking at the correlation matrix of the data:

library(pheatmap)

R <- cor(job)

pheatmap(R, breaks = seq(-1, 1, length.out = 100))

Figure 1: Correlation between the variables in the ``Job’’ dataset

The correlation between the variables seems to be quite large (it doesn’t matter about direction here, only magnitude; if negative correlations were present, we would think in absolute value).

There appears to be a group of highly correlated variables comprising physical ability, appearance, communication skills, and learning ability which are correlated among themselves but uncorrelated with another group of variables. The second group comprises problem solving and logical ability.

This suggests that PCA might be useful in this problem to reduce the dimensionality without a significant loss of information.

Cov vs Cor

Should we perform PCA on the covariance or the correlation matrix?

This depends on the variances of the variables in the dataset. If the variables have large differences in their variances, then the variables with the largest variances will tend to dominate the first few principal components.

A solution to this is to standardise the variables prior to computing the covariance matrix - i.e., compute the correlation matrix!

# show that the correlation matrix and the covariance matrix of the standardized variables are identical
all.equal(cor(job), cov(scale(job)))

## [1] TRUE

Question A4

Look at the variance of the variables in the data set. Do you think that PCA should be carried on the covariance matrix or the correlation matrix?

Solution

Let’s have a look at the standard deviation of each variable:

job %>% 
  summarise(across(everything(), sd))

## # A tibble: 1 x 6
##   commun probl_solv logical learn physical appearance
##    <dbl>      <dbl>   <dbl> <dbl>    <dbl>      <dbl>
## 1   2.74       2.41    2.49  2.84     2.41       2.99

As the standard deviations appear to be fairly similar (and so will the variances) we can perform PCA using the covariance matrix.

Perform PCA

Question A5

Using the principal() function from the psych package, perform a PCA of the job performance data, Call the output job_pca.

job_pca <- principal(job, nfactors = ncol(job), covar = ..., rotate = 'none')
job_pca$loadings

Depending on your answer to the previous question, either set covar = TRUE or covar = FALSE within the principal() function.

Warning: the output of the function will be in terms of standardized variables nevertheless. So you will see output with standard deviation of 1.

Solution

library(psych)

job_pca <- principal(job, nfactors = ncol(job), covar = TRUE, rotate = 'none')

The output

job_pca$loadings

## 
## Loadings:
##            PC1    PC2    PC3    PC4    PC5    PC6   
## commun      0.984 -0.120                0.101       
## probl_solv  0.223  0.810  0.543                     
## logical     0.329  0.747 -0.578                     
## learn       0.987 -0.110                       0.105
## physical    0.988                      -0.110       
## appearance  0.979 -0.125         0.161              
## 
##                  PC1   PC2   PC3   PC4   PC5   PC6
## SS loadings    4.035 1.261 0.631 0.035 0.022 0.016
## Proportion Var 0.673 0.210 0.105 0.006 0.004 0.003
## Cumulative Var 0.673 0.883 0.988 0.994 0.997 1.000

The output is made up of two parts.

First, it shows the loading matrix. In each column of the loading matrix we find how much each of the measured variables contributes to the computed new axis/direction (that is, the principal component). Notice that there are as many principal components as variables.

The second part of the output displays the contribution of each component to the total variance.

Before interpreting it however, let’s focus on the last row of that output called “Cumulative Var.” This displays the cumulative sum of the variances of each principal component. Taken all together, the six principal components taken explain all of the total variance in the original data. In other words, the total variance of the principal components (the sum of their variances) is equal to the total variance in the original data (the sum of the variances of the variables).

However, our goal is to reduce the dimensionality of our data, so it comes natural to wonder which of the six principal components explain most of the variability, and which components instead do not contribute substantially to the total variance.

To that end, the second row “Proportion Var” displays the proportion of the total variance explained by each component, i.e. the variance of the principal component divided by the total variance.

The last row, as we saw, is the cumulative proportion of explained variance: 0.673, 0.673 + 0.210, 0.673 + 0.210 + 0.105, and so on.

We also notice that the first PC alone explains 67.3% of the total variability, while the first two components together explain almost 90% of the total variability. From the third component onwards, we do not see such a sharp increase in the proportion of explained variance, and the cumulative proportion slowly reaches the total ratio of 1 (or 100%).

Optional: (some of) the math behind it

Doing data reduction can feel a bit like magic, and in part that’s just because it’s quite complicated.

The intuition

Consider one way we might construct a correlation matrix - as the product of vector \(\mathbf{f}\) with \(\mathbf{f'}\) (f transposed): \[ \begin{equation*} \mathbf{f} = \begin{bmatrix} 0.9 \\ 0.8 \\ 0.7 \\ 0.6 \\ 0.5 \\ 0.4 \\ \end{bmatrix} \qquad \mathbf{f} \mathbf{f'} = \begin{bmatrix} 0.9 \\ 0.8 \\ 0.7 \\ 0.6 \\ 0.5 \\ 0.4 \\ \end{bmatrix} \begin{bmatrix} 0.9, 0.8, 0.7, 0.6, 0.5, 0.4 \\ \end{bmatrix} \qquad = \qquad \begin{bmatrix} 0.81, 0.72, 0.63, 0.54, 0.45, 0.36 \\ 0.72, 0.64, 0.56, 0.48, 0.40, 0.32 \\ 0.63, 0.56, 0.49, 0.42, 0.35, 0.28 \\ 0.54, 0.48, 0.42, 0.36, 0.30, 0.24 \\ 0.45, 0.40, 0.35, 0.30, 0.25, 0.20 \\ 0.36, 0.32, 0.28, 0.24, 0.20, 0.16 \\ \end{bmatrix} \end{equation*} \]

But we constrain this such that the diagonal has values of 1 (the correlation of a variable with itself is 1), and lets call it R. \[ \begin{equation*} \mathbf{R} = \begin{bmatrix} 1.00, 0.72, 0.63, 0.54, 0.45, 0.36 \\ 0.72, 1.00, 0.56, 0.48, 0.40, 0.32 \\ 0.63, 0.56, 1.00, 0.42, 0.35, 0.28 \\ 0.54, 0.48, 0.42, 1.00, 0.30, 0.24 \\ 0.45, 0.40, 0.35, 0.30, 1.00, 0.20 \\ 0.36, 0.32, 0.28, 0.24, 0.20, 1.00 \\ \end{bmatrix} \end{equation*} \]

PCA is about trying to determine a vector f which generates the correlation matrix R. a bit like unscrambling eggs!

in PCA, we express \(\mathbf{R = CC'}\), where \(\mathbf{C}\) are our principal components.
If \(n\) is number of variables in \(R\), then \(i^{th}\) component \(C_i\) is the linear sum of each variable multiplied by some weighting:
\[ C_i = \sum_{j=1}^{n}w_{ij}x_{j} \]

How do we find \(C\)?

This is where “eigen decomposition” comes in.
For the \(n \times n\) correlation matrix \(\mathbf{R}\), there is an eigenvector \(x_i\) that solves the equation \[ \mathbf{x_i R} = \lambda_i \mathbf{x_i} \] Where the vector multiplied by the correlation matrix is equal to some eigenvalue \(\lambda_i\) multiplied by that vector.
We can write this without subscript \(i\) as: \[ \begin{align} & \mathbf{R X} = \mathbf{X \lambda} \\ & \text{where:} \\ & \mathbf{R} = \text{correlation matrix} \\ & \mathbf{X} = \text{matrix of eigenvectors} \\ & \mathbf{\lambda} = \text{vector of eigenvalues} \end{align} \] the vectors which make up \(\mathbf{X}\) must be orthogonal (\(\mathbf{XX' = I}\)), which means that \(\mathbf{R = X \lambda X'}\)

We can actually do this in R manually. Creating a correlation matrix

# lets create a correlation matrix, as the produce of ff'
f <- seq(.9,.4,-.1)
R <- f %*% t(f)
#give rownames and colnames
rownames(R)<-colnames(R)<-paste0("V",seq(1:6))
#constrain diagonals to equal 1
diag(R)<-1
R

##      V1   V2   V3   V4   V5   V6
## V1 1.00 0.72 0.63 0.54 0.45 0.36
## V2 0.72 1.00 0.56 0.48 0.40 0.32
## V3 0.63 0.56 1.00 0.42 0.35 0.28
## V4 0.54 0.48 0.42 1.00 0.30 0.24
## V5 0.45 0.40 0.35 0.30 1.00 0.20
## V6 0.36 0.32 0.28 0.24 0.20 1.00

Eigen Decomposition

# do eigen decomposition
e <- eigen(R)
print(e, digits=2)

## eigen() decomposition
## $values
## [1] 3.16 0.82 0.72 0.59 0.44 0.26
## 
## $vectors
##       [,1]   [,2]   [,3]  [,4]   [,5]   [,6]
## [1,] -0.50 -0.061  0.092  0.14  0.238  0.816
## [2,] -0.47 -0.074  0.121  0.21  0.657 -0.533
## [3,] -0.43 -0.096  0.182  0.53 -0.675 -0.184
## [4,] -0.39 -0.142  0.414 -0.78 -0.201 -0.104
## [5,] -0.34 -0.299 -0.860 -0.20 -0.108 -0.067
## [6,] -0.28  0.934 -0.178 -0.10 -0.067 -0.045

The eigenvectors are orthogonal (\(\mathbf{XX' = I}\)):

round(e$vectors %*% t(e$vectors),2)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    0    0    0    0
## [2,]    0    1    0    0    0    0
## [3,]    0    0    1    0    0    0
## [4,]    0    0    0    1    0    0
## [5,]    0    0    0    0    1    0
## [6,]    0    0    0    0    0    1

The Principal Components \(\mathbf{C}\) are the eigenvectors scaled by the square root of the eigenvalues:

#eigenvectors
e$vectors

##            [,1]        [,2]        [,3]       [,4]        [,5]        [,6]
## [1,] -0.4964630 -0.06108095  0.09230163  0.1390722  0.23845173  0.81552052
## [2,] -0.4680485 -0.07428828  0.12095171  0.2141704  0.65661667 -0.53269894
## [3,] -0.4326827 -0.09628978  0.18197078  0.5298520 -0.67508185 -0.18417914
## [4,] -0.3899672 -0.14160502  0.41427181 -0.7780083 -0.20056923 -0.10357339
## [5,] -0.3397764 -0.29920322 -0.86039187 -0.1967243 -0.10763437 -0.06685542
## [6,] -0.2823444  0.93375805 -0.17844447 -0.1002466 -0.06668308 -0.04515620

#scaled by sqrt of eigenvalues
diag(sqrt(e$values))

##         [,1]      [,2]      [,3]      [,4]      [,5]     [,6]
## [1,] 1.77897 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## [2,] 0.00000 0.9064417 0.0000000 0.0000000 0.0000000 0.000000
## [3,] 0.00000 0.0000000 0.8476448 0.0000000 0.0000000 0.000000
## [4,] 0.00000 0.0000000 0.0000000 0.7694699 0.0000000 0.000000
## [5,] 0.00000 0.0000000 0.0000000 0.0000000 0.6641042 0.000000
## [6,] 0.00000 0.0000000 0.0000000 0.0000000 0.0000000 0.511868

C <- e$vectors %*% diag(sqrt(e$values))
C

##            [,1]        [,2]       [,3]        [,4]        [,5]        [,6]
## [1,] -0.8831928 -0.05536631  0.0782390  0.10701187  0.15835679  0.41743889
## [2,] -0.8326442 -0.06733799  0.1025241  0.16479768  0.43606188 -0.27267156
## [3,] -0.7697296 -0.08728107  0.1542466  0.40770515 -0.44832468 -0.09427542
## [4,] -0.6937401 -0.12835669  0.3511553 -0.59865399 -0.13319887 -0.05301591
## [5,] -0.6044520 -0.27121026 -0.7293067 -0.15137345 -0.07148043 -0.03422115
## [6,] -0.5022823  0.84639719 -0.1512575 -0.07713676 -0.04428451 -0.02311401

And we can reproduce our correlation matrix, because \(\mathbf{R = CC'}\).

C %*% t(C)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.00 0.72 0.63 0.54 0.45 0.36
## [2,] 0.72 1.00 0.56 0.48 0.40 0.32
## [3,] 0.63 0.56 1.00 0.42 0.35 0.28
## [4,] 0.54 0.48 0.42 1.00 0.30 0.24
## [5,] 0.45 0.40 0.35 0.30 1.00 0.20
## [6,] 0.36 0.32 0.28 0.24 0.20 1.00

Now lets imagine we only consider 1 principal component.
We can do this with the principal() function:

library(psych)
pc1<-principal(R, nfactors = 1, covar = FALSE, rotate = 'none')
pc1

## Principal Components Analysis
## Call: principal(r = R, nfactors = 1, rotate = "none", covar = FALSE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##     PC1   h2   u2 com
## V1 0.88 0.78 0.22   1
## V2 0.83 0.69 0.31   1
## V3 0.77 0.59 0.41   1
## V4 0.69 0.48 0.52   1
## V5 0.60 0.37 0.63   1
## V6 0.50 0.25 0.75   1
## 
##                 PC1
## SS loadings    3.16
## Proportion Var 0.53
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 component is sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.09 
## 
## Fit based upon off diagonal values = 0.95

Look familiar? It looks like the first component we computed manually. The first column of \(\mathbf{C}\):

cbind(pc1$loadings, C=C[,1])

##          PC1          C
## V1 0.8831928 -0.8831928
## V2 0.8326442 -0.8326442
## V3 0.7697296 -0.7697296
## V4 0.6937401 -0.6937401
## V5 0.6044520 -0.6044520
## V6 0.5022823 -0.5022823

We can now ask “how well does this component (on its own) recreate our correlation matrix?”

C[,1] %*% t(C[,1])

##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.7800296 0.7353854 0.6798197 0.6127063 0.5338477 0.4436121
## [2,] 0.7353854 0.6932964 0.6409109 0.5776386 0.5032935 0.4182224
## [3,] 0.6798197 0.6409109 0.5924836 0.5339923 0.4652646 0.3866215
## [4,] 0.6127063 0.5776386 0.5339923 0.4812753 0.4193326 0.3484533
## [5,] 0.5338477 0.5032935 0.4652646 0.4193326 0.3653623 0.3036056
## [6,] 0.4436121 0.4182224 0.3866215 0.3484533 0.3036056 0.2522875

It looks close, but not quite. How much not quite? Measurably so!

R - (C[,1] %*% t(C[,1]))

##             V1          V2          V3          V4          V5          V6
## V1  0.21997041 -0.01538541 -0.04981966 -0.07270625 -0.08384772 -0.08361212
## V2 -0.01538541  0.30670361 -0.08091089 -0.09763864 -0.10329350 -0.09822244
## V3 -0.04981966 -0.08091089  0.40751635 -0.11399225 -0.11526463 -0.10662154
## V4 -0.07270625 -0.09763864 -0.11399225  0.51872473 -0.11933260 -0.10845334
## V5 -0.08384772 -0.10329350 -0.11526463 -0.11933260  0.63463772 -0.10360556
## V6 -0.08361212 -0.09822244 -0.10662154 -0.10845334 -0.10360556  0.74771250

Notice the values on the diagonals?

diag(C[,1] %*% t(C[,1]))

## [1] 0.7800296 0.6932964 0.5924836 0.4812753 0.3653623 0.2522875

diag(R) - diag(C[,1] %*% t(C[,1]))

##        V1        V2        V3        V4        V5        V6 
## 0.2199704 0.3067036 0.4075164 0.5187247 0.6346377 0.7477125

pc1$communality

##        V1        V2        V3        V4        V5        V6 
## 0.7800296 0.6932964 0.5924836 0.4812753 0.3653623 0.2522875

pc1$uniquenesses

##        V1        V2        V3        V4        V5        V6 
## 0.2199704 0.3067036 0.4075164 0.5187247 0.6346377 0.7477125

How many components to keep?

There is no single best method to select the optimal number of components to keep, while discarding the remaining ones (which are then considered as noise components).

The following three heuristic rules are commonly used in the literature:

• The cumulative proportion of explained variance criterion
• Kaiser’s rule
• The scree plot
• Velicer’s Minimum Average Partial method
• Parallel analysis

In the next sections we will analyse each of them in turn.

The cumulative proportion of explained variance criterion

The rule suggests to keep as many principal components as needed in order to explain approximately 80-90% of the total variance.

Question A6

Looking again at the PCA output, how many principal components would you keep if you were following the cumulative proportion of explained variance criterion?

Solution

Let’s look again at the PCA summary:

job_pca$loadings

## 
## Loadings:
##            PC1    PC2    PC3    PC4    PC5    PC6   
## commun      0.984 -0.120                0.101       
## probl_solv  0.223  0.810  0.543                     
## logical     0.329  0.747 -0.578                     
## learn       0.987 -0.110                       0.105
## physical    0.988                      -0.110       
## appearance  0.979 -0.125         0.161              
## 
##                  PC1   PC2   PC3   PC4   PC5   PC6
## SS loadings    4.035 1.261 0.631 0.035 0.022 0.016
## Proportion Var 0.673 0.210 0.105 0.006 0.004 0.003
## Cumulative Var 0.673 0.883 0.988 0.994 0.997 1.000

The following part of the output tells us that the first two components explain 88.3% of the total variance.

Cumulative Var 0.673 0.883 0.988 0.994 0.997 1.000

According to this criterion, we should keep 2 principal components.

Kaiser’s rule

According to Kaiser’s rule, we should keep the principal components having variance larger than 1. Standardized variables have a variance equal 1. Because we have 6 variables in the data set, and the total variance is 6, the value 1 represents the average variance in the data: \[ \frac{1 + 1 + 1 + 1 + 1 + 1}{6} = 1 \]

Hint:

The variances of each PC are shown in the row of the output named SS loadings and also in job_pca$values. The average variance is:

mean(job_pca$values)

## [1] 1

Question A7

Looking again at the PCA output, how many principal components would you keep if you were following Kaiser’s criterion?

Solution

job_pca$loadings

## 
## Loadings:
##            PC1    PC2    PC3    PC4    PC5    PC6   
## commun      0.984 -0.120                0.101       
## probl_solv  0.223  0.810  0.543                     
## logical     0.329  0.747 -0.578                     
## learn       0.987 -0.110                       0.105
## physical    0.988                      -0.110       
## appearance  0.979 -0.125         0.161              
## 
##                  PC1   PC2   PC3   PC4   PC5   PC6
## SS loadings    4.035 1.261 0.631 0.035 0.022 0.016
## Proportion Var 0.673 0.210 0.105 0.006 0.004 0.003
## Cumulative Var 0.673 0.883 0.988 0.994 0.997 1.000

The variances are shown in the row

SS loadings    4.035 1.261 0.631 0.035 0.022 0.016

From the result we see that only the first two principal components have variance greater than 1, so this rule suggests to keep 2 PCs only.

The scree plot

The scree plot is a graphical criterion which involves plotting the variance for each principal component. This can be easily done by calling plot on the variances, which are stored in job_pca$values

plot(x = 1:length(job_pca$values), y = job_pca$values, 
     type = 'b', xlab = '', ylab = 'Variance', 
     main = 'Police officers: scree plot', frame.plot = FALSE)

where the argument type = 'b' tells R that the plot should have both points and lines.

A typical scree plot features higher variances for the initial components and quickly drops to small variances where the curve is almost flat. The flat part of the curve represents the noise components, which are not able to capture the main sources of variability in the system.

According to the scree plot criterion, we should keep as many principal components as where the “elbow” in the plot occurs. By elbow we mean the variance before the curve looks almost flat.

Alternatively, some people prefer to use the function scree() from the psych package:

scree(job, factors = FALSE)

This also draws a horizontal line at y = 1. So, if you are making a decision about how many PCs to keep by looking at where the plot falls below the y = 1 line, you are basically following Kaiser’s rule. In fact, Kaiser’s criterion tells you to keep as many PCs as are those with a variance (= eigenvalue) greater than 1.

Question A8

According to the scree plot, how many principal components would you retain?

Solution

This criterion then suggests to keep three principal components.

Velicer’s Minimum Average Partial method

The Minimum Average Partial (MAP) test computes the partial correlation matrix (removing and adjusting for a component from the correlation matrix), sequentially partialling out each component. At each step, the partial correlations are squared and their average is computed.
At first, the components which are removed will be those that are most representative of the shared variance between 2+ variables, meaning that the “average squared partial correlation” will decrease. At some point in the process, the components being removed will begin represent variance that is specific to individual variables, meaning that the average squared partial correlation will increase.
The MAP method is to keep the number of components for which the average squared partial correlation is at the minimum.

We can conduct MAP in R using:

VSS(data, plot = FALSE, method="pc", n = ncol(data))

(be aware there is a lot of other information in this output too! For now just focus on the map column)

Question A9

How many components should we keep according to the MAP method?

Solution

job_map <- VSS(job, plot=FALSE, method="pc", n = ncol(job))$map
paste("MAP is lowest for", which.min(job_map), "components")

## [1] "MAP is lowest for 2 components"

According to the MAP criterion we should keep 2 principal components.

Parallel analysis

Parallel analysis involves simulating lots of datasets of the same dimension but in which the variables are uncorrelated. For each of these simulations, a PCA is conducted on its correlation matrix, and the eigenvalues are extracted. We can then compare our eigenvalues from the PCA on our actual data to the average eigenvalues across these simulations. In theory, for uncorrelated variables, no components should explain more variance than any others, and eigenvalues should be equal to 1. In reality, variables are rarely truly uncorrelated, and so there will be slight variation in the magnitude of eigenvalues simply due to chance. The parallel analysis method suggests keeping those components for which the eigenvalues are greater than those from the simulations.

It can be conducted in R using:

fa.parallel(job, fa="pc", quant=.95)

Question A10

How many components should we keep according to parallel analysis?

Solution

fa.parallel(job, fa="pc", quant=.95)

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  1

Parallel analysis suggests to keep 1 principal component only as there is only one PC with an eigenvalue higher than the simulated random ones in red.

Interpretation

Because three out of the five selection criteria introduced above suggest to keep 2 principal components, in the following we will work with the first two PCs only.

Let’s have a look at the selected principal components:

job_pca$loadings[, 1:2]

##                  PC1         PC2
## commun     0.9840708 -0.11972298
## probl_solv 0.2233011  0.80950107
## logical    0.3287266  0.74661139
## learn      0.9869336 -0.10970488
## physical   0.9883609 -0.07835277
## appearance 0.9787497 -0.12532863

and at their corresponding proportion of total variance explained:

job_pca$values / sum(job_pca$values)

## [1] 0.672527675 0.210155954 0.105099073 0.005767644 0.003721180 0.002728475

We see that the first PC accounts for 67.3% of the total variability. All loadings seem to have the same magnitude apart from probl_solv and logical which are closer to zero. The first component looks like a sort of average of the officers performance scores excluding problem solving and logical ability.

The second principal component, which explains only 21% of the total variance, has two loadings clearly distant from zero: the ones associated to problem solving and logical ability. It distinguishes police officers with strong logical and problem solving skills and a low score on the test (note the negative magnitude) from the other officers.

We have just seen how to interpret the first components by looking at the magnitude and sign of the coefficients for each measured variable.

For interpretation purposes, it might help hiding very small scores in the loadings. This can be done by specifying the cutoff value in the print() function. However, this only works when you pass the loadings for all the PCs:

print(job_pca$loadings, cutoff = 0.3)

## 
## Loadings:
##            PC1    PC2    PC3    PC4    PC5    PC6   
## commun      0.984                                   
## probl_solv         0.810  0.543                     
## logical     0.329  0.747 -0.578                     
## learn       0.987                                   
## physical    0.988                                   
## appearance  0.979                                   
## 
##                  PC1   PC2   PC3   PC4   PC5   PC6
## SS loadings    4.035 1.261 0.631 0.035 0.022 0.016
## Proportion Var 0.673 0.210 0.105 0.006 0.004 0.003
## Cumulative Var 0.673 0.883 0.988 0.994 0.997 1.000

In the literature, some authors also suggest to look at the correlation between each principal component and the measured variables:

# First PC
cor(job_pca$scores[,1], job)

##         commun probl_solv   logical     learn  physical appearance
## [1,] 0.9854237  0.2137116 0.3187957 0.9882455 0.9885528  0.9809939

The first PC is strongly correlated with all the measured variables except probl_solv and logical. As we mentioned above, all variables seem to contributed to the first PC.

# Second PC
cor(job_pca$scores[,2], job)

##          commun probl_solv   logical      learn   physical appearance
## [1,] -0.1634967  0.7924169 0.7378895 -0.1536895 -0.1220727 -0.1693193

The second PC is strongly correlated with probl_solv and logical, and slightly negatively correlated with the remaining variables. This separates police offices with clear logical and problem solving skills and a small score on the test (negative sign) from the others.

Plotting the retained principal components

We can now visualise the statistical units (police officers) in the reduced space given by the retained principal components.

tibble(pc1 = job_pca$scores[, 1],
       pc2 = job_pca$scores[, 2]) %>%
  ggplot(.,aes(x=pc1,y=pc2))+
  geom_point()

Optional: How well are the units represented in the reduced space?

We now focus our attention on the following question: Are all the statistical units (police officers) well represented in the 2D plot?

The 2D representation of the original data, which comprise 6 measured variables, is an approximation and henceforth it may happen that not all units are well represented in this new space.

Typically, it is good to assess the approximation for each statistical unit by inspecting the scores on the discarded principal components. If a unit has a high score on those components, then this is a sign that the unit might be highly misplaced in the new space and misrepresented.

Consider the 3D example below. There are three cases (= units or individuals). In the original space they are all very different from each other. For example, cases 1 and 2 are very different in their x and y values, but very similar in their z value. Cases 2 and 3 are very similar in their x and y values but very different in their z value. Cases 1 and 3 have very different values for all three variables x, y, and z.

However, when represented in the 2D space given by the two principal components, units 2 and 3 seems like they are very similar when, in fact, they were very different in the original space which also accounted for the z variable.

We typically measure how badly a unit is represented in the new coordinate system by considering the sum of squared scores on the discarded principal components:

scores_discarded <- job_pca$scores[, -(1:2)]
sum_sq <- rowSums(scores_discarded^2)
sum_sq

##  [1]  28.509764  46.889599  63.690141  64.235774  36.577460  17.393796
##  [7]  49.237998  35.095353  18.564496  19.272346  18.564496  24.435681
## [13]  12.389191  59.101079  24.432716  33.184844  13.401839  12.685414
## [19]  11.222085  78.868819  14.155536  34.175498  95.567877  18.399495
## [25]  16.447519  14.408862  31.966172  33.521642  40.122681  32.480075
## [31]  16.851942  24.845993  30.836630  16.001226  29.587409  11.013435
## [37]   8.065003  18.181302  14.604670  23.732693  29.823514  41.371227
## [43]   9.298762  65.418475  21.976763  63.973763  36.092667  84.980382
## [49] 129.650749  88.002871

Units with a high score should be considered for further inspection as it may happen that they are represented as close to another unit when, in fact, they might be very different.

boxplot(sum_sq)

There seem to be only five outliers, and they are not too high compared to the rest of the scores. For this reason, we will consider the 2D representation of the data to be satisfactory.

Exploratory Factor Analysis (EFA)

Where PCA aims to summarise a set of measured variables into a set of orthogonal (uncorrelated) components as linear combinations (a weighted average) of the measured variables, Factor Analysis (FA) assumes that the relationships between a set of measured variables can be explained by a number of underlying latent factors.

Note how the directions of the arrows in Figure 2 are different between PCA and FA - in PCA, each component \(C_i\) is the weighted combination of the observed variables \(y_1, ...,y_n\), whereas in FA, each measured variable \(y_i\) is seen as generated by some latent factor(s) \(F_i\) plus some unexplained variance \(u_i\).

It might help to read the \(\lambda\)s as beta-weights (\(b\), or \(\beta\)), because that’s all they really are. The equation \(y_i = \lambda_{1i}F_1 + \lambda_{2i}F_2 + u_i\) is just our way of saying that the variable \(y_i\) is the manifestation of some amount (\(\lambda_{1i}\)) of an underlying factor \(F_1\), some amount (\(\lambda_{2i}\)) of some other underlying factor \(F_2\), and some error (\(u_i\)).

Figure 2: Path diagrams for PCA and FA

In Exploratory Factor Analysis (EFA), we are starting with no hypothesis about either the number of latent factors or about the specific relationships between latent factors and measured variables (known as the factor structure). Typically, all variables will load on all factors, and a transformation method such as a rotation (we’ll cover this in more detail below) is used to help make the results more easily interpretable.

A brief look to next week

When we have some clear hypothesis about relationships between measured variables and latent factors, we might want to impose a specific factor structure on the data (e.g., items 1 to 10 all measure social anxiety, items 11 to 15 measure health anxiety, and so on). When we impose a specific factor structure, we are doing Confirmatory Factor Analysis (CFA).

This will be the focus of next week, but it’s important to note that in practice EFA is not wholly “exploratory” (your theory will influence the decisions you make) nor is CFA wholly “confirmatory” (in which you will inevitably get tempted to explore how changing your factor structure might improve fit).

Conduct Problems: Data codebook

A researcher is developing a new brief measure of Conduct Problems. She has collected data from n=450 adolescents on 10 items, which cover the following behaviours:

1. Stealing
2. Lying
3. Skipping school
4. Vandalism
5. Breaking curfew
6. Threatening others
7. Bullying
8. Spreading malicious rumours
9. Using a weapon
10. Fighting

Your task is to use the dimension reduction techniques you learned about in the lecture to help inform how to organise the items she has developed into subscales.

The data can be found at https://uoepsy.github.io/data/conduct_probs.csv

Question B1

Read in the dataset from https://uoepsy.github.io/data/conduct_probs.csv.
The first column is clearly an ID column, and it is easiest just to discard this when we are doing factor analysis.

Create a correlation matrix for the items.
Inspect the items to check their suitability for exploratory factor analysis.

• You can use a function such as cor or corr.test(data) (from the psych package) to create the correlation matrix.
  
• The function cortest.bartlett(cor(data), n = nrow(data)) conducts Bartlett’s test that the correlation matrix is proportional to the identity matrix (a matrix of all 0s except for 1s on the diagonal).
  
• You can check linearity of relations using pairs.panels(data) (also from psych), and you can view the histograms on the diagonals, allowing you to check univariate normality (which is usually a good enough proxy for multivariate normality).
• You can check the “factorability” of the correlation matrix using KMO(data) (also from psych!).
  • Rules of thumb:
    • \(0.8 < MSA < 1\): the sampling is adequate
    • \(MSA <0.6\): sampling is not adequate
    • \(MSA \sim 0\): large partial correlations compared to the sum of correlations. Not good for FA

Optional Kaiser’s suggested cuts

• 0.00 to 0.49 unacceptable
  
• 0.50 to 0.59 miserable
  
• 0.60 to 0.69 mediocre
  
• 0.70 to 0.79 middling
  
• 0.80 to 0.89 meritorious
  
• 0.90 to 1.00 marvelous

Solution

df <- read.csv("https://edin.ac/2Vk1BVU")
# discard the first column
df <- df[,-1]

corr.test(df)  

## Call:corr.test(x = df)
## Correlation matrix 
##        item1 item2 item3 item4 item5 item6 item7 item8 item9 item10
## item1   1.00  0.59  0.49  0.48  0.60  0.17  0.30  0.32  0.26   0.20
## item2   0.59  1.00  0.53  0.51  0.66  0.20  0.33  0.30  0.29   0.19
## item3   0.49  0.53  1.00  0.49  0.55  0.15  0.25  0.24  0.25   0.15
## item4   0.48  0.51  0.49  1.00  0.65  0.23  0.29  0.32  0.28   0.25
## item5   0.60  0.66  0.55  0.65  1.00  0.21  0.30  0.29  0.27   0.21
## item6   0.17  0.20  0.15  0.23  0.21  1.00  0.54  0.57  0.41   0.44
## item7   0.30  0.33  0.25  0.29  0.30  0.54  1.00  0.83  0.61   0.58
## item8   0.32  0.30  0.24  0.32  0.29  0.57  0.83  1.00  0.61   0.59
## item9   0.26  0.29  0.25  0.28  0.27  0.41  0.61  0.61  1.00   0.44
## item10  0.20  0.19  0.15  0.25  0.21  0.44  0.58  0.59  0.44   1.00
## Sample Size 
## [1] 450
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##        item1 item2 item3 item4 item5 item6 item7 item8 item9 item10
## item1      0     0     0     0     0     0     0     0     0      0
## item2      0     0     0     0     0     0     0     0     0      0
## item3      0     0     0     0     0     0     0     0     0      0
## item4      0     0     0     0     0     0     0     0     0      0
## item5      0     0     0     0     0     0     0     0     0      0
## item6      0     0     0     0     0     0     0     0     0      0
## item7      0     0     0     0     0     0     0     0     0      0
## item8      0     0     0     0     0     0     0     0     0      0
## item9      0     0     0     0     0     0     0     0     0      0
## item10     0     0     0     0     0     0     0     0     0      0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

cortest.bartlett(cor(df), n=450)

## $chisq
## [1] 2237.533
## 
## $p.value
## [1] 0
## 
## $df
## [1] 45

KMO(df)  

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = df)
## Overall MSA =  0.87
## MSA for each item = 
##  item1  item2  item3  item4  item5  item6  item7  item8  item9 item10 
##   0.90   0.88   0.92   0.88   0.84   0.94   0.82   0.81   0.95   0.94

pairs.panels(df)

or alternatively, if you want a ggplot based approach:

library(GGally)
ggpairs(data=df, diag=list(continuous="density"), axisLabels="show")

How many factors?

Question B2

How many dimensions should be retained? This question can be answered in the same way as we did above for PCA.

Use a scree plot, parallel analysis, and MAP test to guide you.
You can use fa.parallel(data, fm = "fa") to conduct both parallel analysis and view the scree plot!

Solution

fa.parallel(df, fa = "fa")

## Parallel analysis suggests that the number of factors =  2  and the number of components =  NA

In this case the scree plot has a kink at the third factor, so we probably want to retain 2 factors.

We can conduct the MAP test using VSS(data).

VSS(df, plot = FALSE, n = ncol(df))$map

##  [1] 0.10575332 0.03377130 0.05762202 0.10352974 0.14935664 0.25195168
##  [7] 0.39743503 0.45518643 1.00000000         NA

The MAP test suggests retaining 2 factors.

Perform EFA

Now we need to perform the factor analysis. But there are two further things we need to consider, and they a) whether we want to apply a rotation to our factor loadings, in order to make them easier to interpret, and b) how do we want to extract our factors (it turns out there are loads of different approaches!).

Rotations?

Rotations are so called because they transform our loadings matrix in such a way that it can make it more easy to interpret. You can think of it as a transformation applied to our loadings in order to optimise interpretability, by maximising the loading of each item onto one factor, while minimising its loadings to others. We can do this by simple rotations, but maintaining our axes (the factors) as perpendicular (i.e., uncorrelated) as in Figure 4, or we can allow them to be transformed beyond a rotation to allow the factors to correlate (Figure 5).

Figure 3: No rotation

Figure 4: Orthogonal rotation

Figure 5: Oblique rotation

In our path diagram of the model (Figure 6), all the factor loadings remain present, but some of them become negligible. We can also introduce the possible correlation between our factors, as indicated by the curved arrow between \(F_1\) and \(F_2\).

Figure 6: Path diagrams for EFA with rotation

Factor Extraction

PCA (using eigendecomposition) is itself a method of extracting the different dimensions from our data. However, there are lots more available for factor analysis.

You can find a lot of discussion about different methods both in the help documentation for the fa() function from the psych package:

Factoring method fm=“minres” will do a minimum residual as will fm=“uls.” Both of these use a first derivative. fm=“ols” differs very slightly from “minres” in that it minimizes the entire residual matrix using an OLS procedure but uses the empirical first derivative. This will be slower. fm=“wls” will do a weighted least squares (WLS) solution, fm=“gls” does a generalized weighted least squares (GLS), fm=“pa” will do the principal factor solution, fm=“ml” will do a maximum likelihood factor analysis. fm=“minchi” will minimize the sample size weighted chi square when treating pairwise correlations with different number of subjects per pair. fm =“minrank” will do a minimum rank factor analysis. “old.min” will do minimal residual the way it was done prior to April, 2017 (see discussion below). fm=“alpha” will do alpha factor analysis as described in Kaiser and Coffey (1965)

And there are lots of discussions both in papers and on forums.

As you can see, this is a complicated issue, but when you have a large sample size, a large number of variables, for which you have similar communalities, then the extraction methods tend to agree. For now, don’t fret too much about the factor extraction method.²

Question B3

Use the function fa() from the psych package to conduct and EFA to extract 2 factors (this is what we suggest based on the various tests above, but you might feel differently - the ideal number of factors is subjective!). Use a suitable rotation and extraction method (fm).

conduct_efa <- fa(data, nfactors = ?, rotate = ?, fm = ?)

Solution

For example, you could choose an oblimin rotation to allow factors to correlate and use minres as the extraction method.

conduct_efa <- fa(df, nfactors=2, rotate='oblimin', fm="minres")

Inspect

Question B4

Inspect the loadings (conduct_efa$loadings) and give the factors you extracted labels based on the patterns of loadings.

Look back to the description of the items, and suggest a name for your factors

Solution

You can inspect the loadings using:

print(conduct_efa$loadings, sort=T)

## 
## Loadings:
##        MR1    MR2   
## item1   0.613  0.375
## item2   0.654  0.420
## item3   0.554  0.384
## item4   0.619  0.339
## item5   0.693  0.502
## item6   0.527 -0.337
## item7   0.766 -0.458
## item8   0.782 -0.484
## item9   0.613 -0.276
## item10  0.552 -0.360
## 
##                  MR1   MR2
## SS loadings    4.128 1.593
## Proportion Var 0.413 0.159
## Cumulative Var 0.413 0.572

We can see that the first five items have high loadings for one factor and the second five items have high loadings for the other.

The first five items all have in common that they are non-aggressive forms of conduct problems, while the last five items are all aggressive behaviours. We could, therefore, label our factors: ‘non-aggressive’ and ‘aggressive’ conduct problems.

Question B5

How correlated are your factors?

We can inspect the factor correlations (if we used an oblique rotation) using:

conduct_efa$Phi

Solution

conduct_efa$Phi

## NULL

We can see here that there is a moderate correlation between the two factors. An oblique rotation would be appropriate here.

Write-up

Question B6

Drawing on your previous answers and conducting any additional analyses you believe would be necessary to identify an optimal factor structure for the 10 conduct problems, write a brief text that summarises your method and the results from your chosen optimal model.

Solution

The main principles governing the reporting of statistical results are transparency and reproducibility (i.e., someone should be able to reproduce your analysis based on your description).

An example summary would be:

First, the data were checked for their suitability for factor analysis. Normality was checked using visual inspection of histograms, linearity was checked through the inspection of the linear and lowess lines for the pairwise relations of the variables, and factorability was confirmed using a KMO test, which yielded an overall KMO of \(.87\) with no variable KMOs \(<.50\). An exploratory factor analysis was conducted to inform the structure of a new conduct problems test. Inspection of a scree plot alongside parallel analysis (using principal components analysis; PA-PCA) and the MAP test were used to guide the number of factors to retain. All three methods suggested retaining two factors; however, a one-factor and three-factor solution were inspected to confirm that the two-factor solution was optimal from a substantive and practical perspective, e.g., that it neither blurred important factor distinctions nor included a minor factor that would be better combined with the other in a one-factor solution. These factor analyses were conducted using minres extraction and (for the two- and three-factor solutions) an oblimin rotation, because it was expected that the factors would correlate. Inspection of the factor loadings and correlations reinforced that the two-factor solution was optimal: both factors were well-determined, including 5 loadings \(>|0.3|\) and the one-factor model blurred the distinction between different forms of conduct problems. The factor loadings are provided in Table 1³. Based on the pattern of factor loadings, the two factors were labelled ‘aggressive conduct problems’ and ‘non-aggressive conduct problems.’ These factors had a correlation of \(r=.43\). Overall, they accounted for 52% of the variance in the items, suggesting that a two-factor solution effectively summarised the variation in the items.

Table 1: Factor loadings.

	MR1	MR2
item1	0.61	0.38
item2	0.65	0.42
item3	0.55	0.38
item4	0.62	0.34
item5	0.69	0.50
item6	0.53	-0.34
item7	0.77	-0.46
item8	0.78	-0.48
item9	0.61
item10	0.55	-0.36

Comparison exercise

Question B7

Using the same data, conduct a PCA using the principal() function.

What differences do you notice compared to your EFA?

Do you think a PCA or an EFA is more appropriate in this particular case?

Solution

We can use:

principal(df, nfactors=2)

## Principal Components Analysis
## Call: principal(r = df, nfactors = 2)
## Standardized loadings (pattern matrix) based upon correlation matrix
##         RC1  RC2   h2   u2 com
## item1  0.17 0.77 0.62 0.38 1.1
## item2  0.17 0.81 0.68 0.32 1.1
## item3  0.11 0.75 0.58 0.42 1.0
## item4  0.21 0.74 0.60 0.40 1.2
## item5  0.16 0.85 0.75 0.25 1.1
## item6  0.73 0.08 0.53 0.47 1.0
## item7  0.87 0.20 0.80 0.20 1.1
## item8  0.88 0.19 0.82 0.18 1.1
## item9  0.72 0.21 0.56 0.44 1.2
## item10 0.75 0.09 0.57 0.43 1.0
## 
##                        RC1  RC2
## SS loadings           3.29 3.22
## Proportion Var        0.33 0.32
## Cumulative Var        0.33 0.65
## Proportion Explained  0.51 0.49
## Cumulative Proportion 0.51 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.06 
##  with the empirical chi square  166.43  with prob <  1.9e-22 
## 
## Fit based upon off diagonal values = 0.98

We can see that while the loadings differ somewhat between the EFA and the PCA, the overall pattern is quite similar. This is not always the case, especially when the item communalities are low.

In terms of which method is more appropriate, arguably EFA would be more appropriate in this case because our researcher wishes to measure a theoretical construct (conduct problems), rather than simply reduce the dimensions of her data.

---

1. Even if we cut open someone’s brain, it’s unclear what we would be looking for in order to ‘measure’ it. It is unclear whether anxiety even exists as a physical thing, or rather if it is simply the overarching concept we apply to a set of behaviours and feelings↩︎

2. (It’s a bit like the optimiser issue in the multi-level model block)↩︎

3. You should provide the table of factor loadings. It is conventional to omit factor loadings \(<|0.3|\); however, be sure to ensure that you mention this in a table note.↩︎