We’re going to be needing some different packages this week (no more lme4!).
Make sure you have these packages installed:
psych
GPArotation
car
Reducing the dimensionality of job performance
Data: jobperfbonus.csv
A company has asked line managers to fill out a questionnaire that asks them to capture ‘performance indicators’ for each of their employees. These 15 indicators are a mix of metrics from the companies work management software, and questions about manager’s perceptions of their employees’ work over the last year.
The company’s leadership team finds the 15-indicator report a bit confusing and overwhelming. They want to create a simplified “Performance Dashboard” that captures only the most important directions of variance in their employees.
We’re probably going to want to use PCA, because it is hard to understand what an underlying latent dimension of ‘performance’ would be — i.e., it makes more sense to think of any overall measure of ‘performance’ as just being the combination of all of those indicators.
One way to reason about this is to task which of these makes more sense:
A. If my ‘job performance’ increases, then my scores on the indicators will go up
B. If my scores on the indicators go up, then my ‘job performance’ will increase
Really, we are just asking about the direction of the arrows. Which picture makes more sense?
Question 2
Explore the relationships between variables in the data.
You’re probably going to want to subset out just the relevant variables (q1 to q15).
Hints
There are lots of variables! There are various things we can do here. Try some of them to see what they do:
cor(data)
heatmap(cor(data))
pairs.panels(data) (from the psych package)
Let’s keep the full data with the names, but make a new object that is just the performance data:
qperf <- jobperf |>select(-name)
It looks like we’ve got some very strong groups of questions there - q1 to q5 are all highly related to one another, as are q6 to q10, and q11 to q15. Furthermore, the relations between those sets are very weak, suggesting 3 groups that are fairly distinct.
heatmap(cor(qperf))
You’ll probably have to zoom in when you do this yourself, as there are a lot of little plots there!
library(psych)pairs.panels(qperf)
Question 3
How much variance in the set of performance indicators will be captured by 15 principal components?
Note: We can figure this out without having to do anything - it’s a theoretical question!
All of it!
The right hand entry of the “Cumulative Var” row is 1 - it explains everything.
Remember, the company is wanting us to simplify things down from having 15 indicators of ‘job performance’ to a smaller set, that captures the most important ways in which employees vary in their performance.
We’re going to use PCA.
How many components should we keep?
According to the Scree plot, I would suggest keeping 3 components
scree(qperf, factors =FALSE)
According to the MAP, we should keep 3 components
VSS(qperf, n =ncol(qperf), rotate ="none", fm ="pc", plot =FALSE)
Very Simple Structure
Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm,
n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
VSS complexity 1 achieves a maximimum of 0.96 with 3 factors
VSS complexity 2 achieves a maximimum of 0.98 with 6 factors
The Velicer MAP achieves a minimum of 0.02 with 3 factors
BIC achieves a minimum of Inf with factors
Sample Size adjusted BIC achieves a minimum of Inf with factors
Statistics by number of factors
vss1 vss2 map dof chisq prob sqresid fit RMSEA BIC SABIC complex eChisq
1 0.41 0.00 0.143 0 NA NA 2.6e+01 0.41 NA NA NA NA NA
2 0.70 0.70 0.071 0 NA NA 1.3e+01 0.70 NA NA NA NA NA
3 0.96 0.97 0.021 0 NA NA 1.5e+00 0.97 NA NA NA NA NA
4 0.96 0.97 0.031 0 NA NA 1.2e+00 0.97 NA NA NA NA NA
5 0.84 0.97 0.047 0 NA NA 9.0e-01 0.98 NA NA NA NA NA
6 0.84 0.98 0.069 0 NA NA 6.4e-01 0.99 NA NA NA NA NA
7 0.76 0.96 0.106 0 NA NA 4.8e-01 0.99 NA NA NA NA NA
8 0.75 0.91 0.125 0 NA NA 3.3e-01 0.99 NA NA NA NA NA
9 0.67 0.92 0.174 0 NA NA 2.2e-01 0.99 NA NA NA NA NA
10 0.67 0.91 0.211 0 NA NA 1.6e-01 1.00 NA NA NA NA NA
11 0.67 0.90 0.264 0 NA NA 1.1e-01 1.00 NA NA NA NA NA
12 0.61 0.90 0.406 0 NA NA 6.3e-02 1.00 NA NA NA NA NA
13 0.59 0.87 0.580 0 NA NA 2.9e-02 1.00 NA NA NA NA NA
14 0.59 0.87 1.000 0 NA NA 8.2e-03 1.00 NA NA NA NA NA
15 0.59 0.87 NA 0 NA NA 1.5e-27 1.00 NA NA NA NA NA
SRMR eCRMS eBIC
1 NA NA NA
2 NA NA NA
3 NA NA NA
4 NA NA NA
5 NA NA NA
6 NA NA NA
7 NA NA NA
8 NA NA NA
9 NA NA NA
10 NA NA NA
11 NA NA NA
12 NA NA NA
13 NA NA NA
14 NA NA NA
15 NA NA NA
According to the parallel analysis, we should keep 3 components
fa.parallel(qperf, fa="pc", n.iter =500)
Parallel analysis suggests that the number of factors = NA and the number of components = 3
In this case, everything agrees that we should retain 3 components:
guides
suggestion
Kaiser
3
Scree
3
MAP
3
Parallel Analysis
3
Question 5
Conduct a principal components analysis extracting the number of components you decided on from the previous question.
Be sure to set rotate = "none" (a conventional PCA does not use rotations - it is simply about data reduction. The line is a bit blurred here, but once we start introducing rotations, we are moving more towards a form of EFA).
Examine the loadings for the components. By thinking in relation to the questions that were asked (refer back to Table 1), what do you think each component is capturing?
It’s rarely as neat as this, but we can see the patterns of higher loadings for the groups of 5 questions on each component.
PC1 has higher loadings for q1 to q5
PC2 has higher loadings for q6 to q10
PC3 has higher loadings for q11 to q15
Looking back at the questions in Table 1, questions 1 to 5 were all about stuff related to ‘project management’. Questions 6 to 10 were all about collaboration with team members, and questions 11 to 15 were all about technical skill.
principal(qperf, nfactors =3, rotate ="none")
Principal Components Analysis
Call: principal(r = qperf, nfactors = 3, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
PC1 PC2 PC3 h2 u2 com
q1 0.90 0.29 -0.15 0.91 0.089 1.3
q2 0.87 0.22 -0.17 0.83 0.169 1.2
q3 0.75 0.20 -0.17 0.64 0.361 1.3
q4 0.85 0.29 -0.16 0.83 0.167 1.3
q5 0.89 0.28 -0.15 0.89 0.106 1.3
q6 -0.33 0.85 -0.08 0.83 0.167 1.3
q7 -0.32 0.83 -0.11 0.80 0.200 1.3
q8 -0.31 0.82 -0.05 0.77 0.233 1.3
q9 -0.24 0.73 -0.01 0.59 0.411 1.2
q10 -0.29 0.77 -0.11 0.69 0.315 1.3
q11 0.17 0.17 0.89 0.84 0.155 1.1
q12 0.11 0.10 0.76 0.60 0.398 1.1
q13 0.15 0.10 0.80 0.67 0.327 1.1
q14 0.22 0.11 0.75 0.62 0.378 1.2
q15 0.07 0.13 0.80 0.67 0.334 1.1
PC1 PC2 PC3
SS loadings 4.20 3.61 3.38
Proportion Var 0.28 0.24 0.23
Cumulative Var 0.28 0.52 0.75
Proportion Explained 0.38 0.32 0.30
Cumulative Proportion 0.38 0.70 1.00
Mean item complexity = 1.2
Test of the hypothesis that 3 components are sufficient.
The root mean square of the residuals (RMSR) is 0.04
with the empirical chi square 134.75 with prob < 3.9e-07
Fit based upon off diagonal values = 0.99
Question 6
Extract the scores on the principal components.
The company wants to reward teamwork. Pick 10 people they should give a bonus to.
Hints
See Chapter 2#PCA-walkthrough - scores for how to extract scores (a score on each component for each row of our original dataset) using the $scores from the object fitted with principal().
This will contain as many sets of scores as there are components. One of these (given the previous question) might be of use here.
You’ll likely want to join them back to the column of names. So we can figure out who gets the bonus. cbind() or bind_cols() might help here.
The second component seems to capture a lot of ‘teamwork’ related questions. So if we extract scores on that second component, we could then pick the 10 people for whom they are highest on that component.
Let’s combine them with the original data that contains the employee names, so we can figure out who to give bonuses to.
We could do this with cbind(), or bind_cols(),
# first we bind the columns of the scores, back to # the original data which contains the namesjobperf <-bind_cols(principal(qperf, nfactors =3, rotate ="none")$scores, jobperf ) # we can then choose just the 10 people # who have the highest scores on PC2 (the teamwork component)jobperf |>select(name, PC2) |>arrange(desc(PC2))
# A tibble: 400 × 2
name PC2
<chr> <dbl>
1 Robinson, Shayla 2.88
2 Hunt, Maliek 2.63
3 Kearns, Sarah 2.40
4 Billings, Justin 2.39
5 Phillips, Nashelle 2.33
6 al-Masih, Sawada 2.26
7 Vargas, Joseph 2.16
8 Reider, Samuel 2.12
9 Dhanoa, Eunice 2.10
10 Brown, Jameela 1.98
# ℹ 390 more rows
Understanding Conduct Problems
Data: Conduct Problems
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:
Breaking curfew
Vandalism
Skipping school
Bullying
Spreading malicious rumours
Fighting
Lying
Using a weapon
Stealing
Threatening others
Our task is to use the dimension reduction techniques we learned about in the lecture to help inform how to organise the items she has developed into subscales.
Read in the dataset.
Create a correlation matrix for the items, and inspect the items to check their suitability for exploratory factor analysis (see Chapter 3#EFA-initial checks).
cpdata <-read.csv("https://uoepsy.github.io/data/conduct_probs_scale.csv")# discard the first columncpdata <- cpdata[,-1]
Here’s a correlation matrix. There’s no obvious blocks of items here, but we can see that there are some fairly high correlations, as well as some weaker ones. All are positive.
heatmap(cor(cpdata))
The Bartlett’s test comes out with a p-value of 0 (which isn’t possible, but it’s been rounded for some reason). This suggests that we reject the null of this test (that our correlation matrix is proportional to the identity matrix). This is good. It basically means “we have some non-zero correlations”!
cortest.bartlett(cor(cpdata), n=450)
$chisq
[1] 2237.533
$p.value
[1] 0
$df
[1] 45
The overall sampling adequacy is 0.87, which is pretty good! (or rather, which is ‘meritorious’!). MSA for all items is >.8
Finally, all the relationships here look fairly linear:
pairs.panels(cpdata)
Question 8
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, fa = "both") to conduct both parallel analysis and view the scree plot!
The scree plot shows a kink at 3, which suggests retaining 2 components.
scree(cpdata)
The MAP suggests retaining 2 factors. I’m just extracting the actuall map values here to save having to show all the other output. We can see that the 2nd entry is the smallest:
Parallel analysis suggests that the number of factors = 2 and the number of components = 2
Again, a quite clear picture that 2 factors is preferred:
guides
suggestion
Scree
2
MAP
2
Parallel Analysis
2
Question 9
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 (rotate = ?) and extraction method (fm = ?).
myfa <-fa(data, nfactors = ?, rotate = ?, fm = ?)
Hints
Would you expect factors to be correlated? If so, you’ll want an oblique rotation.
If there are two+ underlying factors of ‘conduct problems’, if someone is higher on one of them, would we expect them to therefore be high/low on the other? If so, then we the two factors are correlated!
For example, you could choose an oblimin rotation to allow factors to correlate and use ml as the extraction method.
Inspect your solution. Make sure to look at and think about the loadings, the variance accounted for, and the factor correlations (if estimated).
What we’re doing here is essentially evaluating whether our solution looks theoretically coherent. A big part of this is a subjective decision about the groupings of items, but we would also like the numerical parts of our model to meet certain criteria (see Chapter 3: EFA#what-makes-a-good-factor-solution).
Hints
Just printing an fa object:
myfa <-fa(data, ..... )myfa
Will give you lots and lots of information.
You can extract individual parts using:
myfa$loadings for the loadings
myfa$Vaccounted for the variance accounted for by each factor
myfa$Phi for the factor correlation matrix
You can find a quick guide to reading the fa output here: efa_output.pdf.
Things look pretty good here. Each item has a clear primary loading on to one of the factors, and the complexity for all items is 1 (meaning they’re clearly link to just one of the factors). The h2 column is showing that the 2 factor solution is explaining 38%+ of the variance in each item. Both factors are well determined, having at least 3 salient loadings.
The 2 factors together explain 57% of the variance in the data - both factors explain a similar amount (29% for factor 1, 28% for factor 2).
We can also see that there is a moderate correlation between the two factors. Use of an oblique rotation was appropriate - if the correlation had been very weak, then it might not have differed much from if we used an orthogonal rotation.
We can see that, ordered like this, we have five items that have high loadings for one factor and another five items that have high loadings for the other.
The five items for factor 1 all have in common that they are non-aggressive forms of conduct problems. The five items for factor 2 are all more aggressive behaviours. We could, therefore, label our factors: ‘non-aggressive’ and ‘aggressive’ conduct problems.
Question 12
Compare your three different solutions:
your current solution from the previous questions
one where you fit 1 more factor
one where you fit 1 fewer factors
Which one looks best?
Hints
We’re looking here to assess:
how much variance is accounted for by each solution
do all factors load on 3+ items at a salient level?
do all items have at least one loading at a salient level?
are there any “Heywood cases” (communalities or standardised loadings that are >1)?
should we perhaps remove some of the more complex items?
is the factor structure (items that load on to each factor) coherent, and does it make theoretical sense?
The 1-factor model explains 37% of the variance (as opposed to the 57% explained by the 2 factor solution), and all items load fairly high on the factor. The downside here is that we’re not discerning between different types of conduct problems that we did in the 2 factor solution.
Factor Analysis using method = ml
Call: fa(r = cpdata, nfactors = 1, fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
ML1 h2 u2 com
breaking_curfew 0.42 0.18 0.82 1
vandalism 0.42 0.18 0.82 1
skipping_school 0.35 0.13 0.87 1
bullying 0.89 0.79 0.21 1
spreading_rumours 0.90 0.81 0.19 1
fighting 0.64 0.41 0.59 1
lying 0.43 0.18 0.82 1
using_weapon 0.68 0.46 0.54 1
stealing 0.42 0.17 0.83 1
threatening_others 0.61 0.38 0.62 1
ML1
SS loadings 3.69
Proportion Var 0.37
Mean item complexity = 1
Test of the hypothesis that 1 factor is sufficient.
df null model = 45 with the objective function = 5.03 with Chi Square = 2237.53
df of the model are 35 and the objective function was 1.78
The root mean square of the residuals (RMSR) is 0.19
The df corrected root mean square of the residuals is 0.22
The harmonic n.obs is 450 with the empirical chi square 1464.73 with prob < 1.9e-285
The total n.obs was 450 with Likelihood Chi Square = 789.3 with prob < 3.4e-143
Tucker Lewis Index of factoring reliability = 0.557
RMSEA index = 0.219 and the 90 % confidence intervals are 0.206 0.232
BIC = 575.48
Fit based upon off diagonal values = 0.8
Measures of factor score adequacy
ML1
Correlation of (regression) scores with factors 0.96
Multiple R square of scores with factors 0.92
Minimum correlation of possible factor scores 0.84
The 3-factor model explains 60% of the variance (only 3% more than the 2-factor model). Notably, the third factor is not very clearly defined - it only has 1 salient loading (possibly 2 if we consider the 0.31 to be salient, but that item is primarily loaded on the 2nd factor).
Factor Analysis using method = ml
Call: fa(r = cpdata, nfactors = 3, rotate = "oblimin", fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
ML1 ML2 ML3 h2 u2 com
breaking_curfew -0.02 0.61 0.31 0.71 0.29 1.5
vandalism 0.06 0.12 0.74 0.72 0.28 1.1
skipping_school 0.00 0.56 0.14 0.44 0.56 1.1
bullying 0.90 0.09 -0.10 0.82 0.18 1.0
spreading_rumours 0.92 -0.02 0.03 0.85 0.15 1.0
fighting 0.65 -0.13 0.14 0.43 0.57 1.2
lying 0.02 0.85 -0.06 0.67 0.33 1.0
using_weapon 0.63 0.08 0.02 0.45 0.55 1.0
stealing 0.06 0.69 0.02 0.53 0.47 1.0
threatening_others 0.62 -0.08 0.09 0.39 0.61 1.1
ML1 ML2 ML3
SS loadings 2.93 2.14 0.94
Proportion Var 0.29 0.21 0.09
Cumulative Var 0.29 0.51 0.60
Proportion Explained 0.49 0.36 0.16
Cumulative Proportion 0.49 0.84 1.00
With factor correlations of
ML1 ML2 ML3
ML1 1.00 0.39 0.32
ML2 0.39 1.00 0.68
ML3 0.32 0.68 1.00
Mean item complexity = 1.1
Test of the hypothesis that 3 factors are sufficient.
df null model = 45 with the objective function = 5.03 with Chi Square = 2237.53
df of the model are 18 and the objective function was 0.02
The root mean square of the residuals (RMSR) is 0.01
The df corrected root mean square of the residuals is 0.02
The harmonic n.obs is 450 with the empirical chi square 3.98 with prob < 1
The total n.obs was 450 with Likelihood Chi Square = 10.51 with prob < 0.91
Tucker Lewis Index of factoring reliability = 1.009
RMSEA index = 0 and the 90 % confidence intervals are 0 0.016
BIC = -99.46
Fit based upon off diagonal values = 1
Measures of factor score adequacy
ML1 ML2 ML3
Correlation of (regression) scores with factors 0.96 0.93 0.88
Multiple R square of scores with factors 0.93 0.86 0.77
Minimum correlation of possible factor scores 0.85 0.72 0.53
Question 13
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 some bullet points to summarise your methods and the results of your chosen optimal model.
Remember, the main principles governing the reporting of statistical methods are transparency and reproducibility (i.e., someone should be able to reproduce your analysis based on your description).
Solution 1. 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 maximum likelihood estimation 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 21. 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 57% of the variance in the items, suggesting that a two-factor solution effectively summarised the over half of the variation in the items.
Table 2: Factor Loadings
ML1
ML2
spreading_rumours
0.93
bullying
0.90
fighting
0.65
using_weapon
0.63
threatening_others
0.62
vandalism
0.88
stealing
0.77
lying
0.70
skipping_school
0.69
breaking_curfew
0.67
Footnotes
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.↩︎
