How can I create a parallel analysis scree in R?

Question

Dataset Dput

structure(list(V1 = structure(c(4, 4, 2, 2, 2, 2, 2, 2, 4, 4, 
2, 3, 2, 3, 4, 2, 2, 2, 3, 3, 2, 3, 1, 3, 3, 3, 3, 4, 1, 2, 4, 
1, 2, 3, 2, 3, 1, 1, 2, 2, 4, 3, 2, 1, 2, 3, 3, 4, 3, 3, 2, 3, 
1, 4, 3, 2, 3, 4, 1, 3, 3, 3, 2, 2, 1, 2, 3, 4, 4, 2, 4, 3, 2, 
3, 3, 3, 3, 2, 4, 3, 3, 3, 2, 2, 3, 4, 2, 4, 4, 2, 2, 3, 3), format.spss = "F8.0"), 
    V2 = structure(c(4, 4, 3, 4, 3, 4, 3, 2, 4, 1, 3, 3, 3, 4, 
    3, 3, 2, 3, 4, 3, 1, 4, 2, 3, 4, 2, 4, 3, 3, 2, 3, 2, 3, 
    3, 4, 3, 3, 3, 3, 3, 3, 2, 4, 2, 2, 2, 4, 3, 4, 4, 2, 4, 
    2, 3, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 3, 3, 4, 4, 4, 4, 4, 
    3, 4, 3, 3, 3, 4, 2, 4, 3, 4, 3, 3, 2, 3, 3, 4, 3, 4, 3, 
    4, 4, 3), format.spss = "F8.0"), V3 = structure(c(4, 4, 4, 
    4, 4, 4, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
    3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4), format.spss = "F8.0"), 
    V4 = structure(c(4, 4, 3, 4, 3, 4, 2, 1, 3, 2, 3, 1, 4, 4, 
    2, 3, 2, 2, 2, 4, 1, 2, 2, 2, 3, 2, 3, 2, 2, 1, 3, 1, 1, 
    2, 4, 1, 1, 2, 3, 2, 2, 1, 1, 1, 3, 2, 4, 3, 3, 3, 3, 3, 
    3, 4, 3, 1, 4, 3, 4, 3, 2, 3, 2, 1, 4, 1, 4, 1, 2, 4, 4, 
    4, 3, 3, 3, 2, 2, 1, 4, 3, 2, 3, 2, 1, 3, 4, 1, 2, 4, 3, 
    4, 2, 2), format.spss = "F8.0"), V5 = structure(c(3, 3, 3, 
    4, 3, 4, 3, 1, 1, 1, 1, 2, 1, 2, 2, 2, 1, 2, 2, 2, 3, 2, 
    2, 2, 2, 4, 2, 3, 2, 3, 4, 1, 4, 2, 3, 3, 2, 2, 3, 2, 2, 
    3, 3, 2, 3, 3, 3, 2, 2, 2, 3, 2, 3, 3, 2, 2, 3, 3, 2, 3, 
    2, 2, 3, 3, 3, 2, 3, 3, 3, 4, 3, 2, 3, 3, 3, 3, 3, 3, 4, 
    3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 4, 3, 3), format.spss = "F8.0"), 
    V6 = structure(c(4, 4, 3, 4, 3, 4, 4, 1, 3, 3, 3, 3, 2, 3, 
    4, 2, 4, 3, 3, 3, 3, 4, 4, 3, 3, 3, 4, 4, 4, 3, 4, 4, 3, 
    3, 3, 4, 2, 2, 3, 3, 3, 4, 2, 4, 3, 4, 4, 4, 3, 4, 2, 4, 
    3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 3, 1, 4, 4, 4, 4, 4, 4, 
    4, 3, 4, 4, 4, 4, 2, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 3, 
    4, 4, 4), format.spss = "F8.0"), V7 = structure(c(4, 4, 2, 
    4, 2, 4, 4, 3, 3, 3, 2, 2, 4, 4, 3, 3, 1, 4, 3, 3, 1, 2, 
    4, 3, 4, 2, 4, 4, 3, 3, 2, 2, 3, 2, 4, 3, 3, 3, 3, 3, 3, 
    1, 4, 3, 2, 2, 4, 3, 4, 4, 2, 4, 2, 3, 4, 3, 3, 3, 4, 3, 
    4, 4, 3, 4, 4, 3, 4, 4, 4, 4, 3, 4, 4, 4, 3, 3, 4, 3, 4, 
    3, 3, 3, 3, 2, 2, 4, 4, 4, 4, 2, 4, 4, 3), format.spss = "F8.0"), 
    V8 = structure(c(4, 4, 2, 1, 2, 1, 1, 1, 3, 3, 2, 3, 2, 3, 
    4, 2, 2, 2, 3, 3, 2, 3, 1, 3, 3, 3, 3, 4, 1, 2, 4, 1, 2, 
    3, 2, 3, 1, 1, 2, 2, 3, 1, 1, 1, 2, 3, 3, 4, 3, 3, 2, 3, 
    1, 3, 4, 2, 3, 4, 1, 3, 3, 3, 2, 2, 1, 2, 3, 4, 4, 2, 4, 
    3, 4, 4, 4, 4, 3, 2, 4, 3, 3, 3, 2, 2, 3, 4, 2, 4, 4, 2, 
    1, 3, 4), format.spss = "F8.0"), V9 = structure(c(4, 4, 4, 
    4, 4, 4, 4, 4, 3, 3, 2, 3, 3, 3, 3, 2, 3, 3, 2, 3, 4, 4, 
    4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 4, 3, 2, 4, 3, 4, 
    4, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 4, 3, 4, 3, 2, 4, 
    3, 3, 4, 4, 4, 3, 4, 4, 4, 4, 4, 3, 4, 3, 4, 3, 4, 4, 4, 
    4, 3, 4, 4, 4, 4, 4, 3, 2, 4, 4, 4, 4, 4), format.spss = "F8.0"), 
    V10 = structure(c(4, 4, 2, 4, 2, 4, 3, 2, 3, 3, 3, 2, 4, 
    4, 2, 2, 1, 3, 4, 4, 1, 4, 2, 3, 3, 2, 4, 3, 2, 3, 3, 1, 
    3, 2, 4, 3, 2, 3, 3, 3, 3, 1, 2, 4, 2, 3, 4, 4, 3, 3, 2, 
    4, 2, 4, 3, 3, 4, 3, 4, 3, 4, 4, 4, 1, 4, 3, 3, 4, 3, 4, 
    4, 3, 3, 3, 3, 3, 4, 1, 4, 3, 3, 3, 3, 2, 3, 4, 4, 2, 4, 
    2, 4, 4, 3), format.spss = "F8.0"), V11 = structure(c(3, 
    3, 1, 4, 1, 4, 1, 1, 1, 1, 2, 1, 1, 1, 3, 2, 2, 2, 2, 1, 
    2, 3, 1, 2, 3, 3, 2, 1, 2, 2, 2, 3, 2, 2, 3, 2, 1, 2, 2, 
    1, 1, 4, 3, 1, 3, 2, 3, 1, 2, 1, 2, 1, 2, 2, 1, 2, 2, 3, 
    2, 2, 2, 2, 2, 2, 1, 1, 1, 3, 3, 4, 2, 1, 2, 2, 3, 3, 3, 
    3, 4, 3, 2, 3, 3, 2, 2, 2, 2, 1, 3, 1, 4, 1, 3), format.spss = "F8.0"), 
    V12 = structure(c(4, 4, 3, 2, 3, 2, 3, 1, 3, 3, 3, 3, 2, 
    3, 3, 2, 4, 3, 3, 4, 4, 3, 3, 4, 4, 3, 3, 3, 4, 3, 4, 4, 
    3, 3, 3, 4, 2, 2, 3, 3, 3, 4, 2, 4, 3, 4, 4, 4, 3, 4, 2, 
    4, 3, 3, 3, 3, 4, 3, 3, 2, 2, 1, 1, 3, 1, 4, 4, 4, 4, 4, 
    4, 4, 3, 3, 2, 2, 2, 2, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 
    3, 2, 3, 4), format.spss = "F8.0")), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -93L))

EFA Before Scree Plot

I have done the work of performing exploratory factor analysis on the data:

fa3 <- fa(hwk2,
   nfactors = 3,
   n.obs = 93,
   rotate = "oblimin",
   max.iter = 100)
fa3

Which gives me this:

      MR1   MR3   MR2   h2   u2 com
V1   0.03  0.87 -0.05 0.77 0.23 1.0
V2   0.75  0.05  0.09 0.63 0.37 1.0
V3   0.13  0.06  0.67 0.53 0.47 1.1
V4   0.50  0.07  0.08 0.31 0.69 1.1
V5   0.03 -0.06  0.88 0.77 0.23 1.0
V6   0.00  0.47  0.32 0.37 0.63 1.8
V7   0.80 -0.08 -0.04 0.60 0.40 1.0
V8   0.05  0.88 -0.03 0.80 0.20 1.0
V9  -0.22  0.02  0.58 0.34 0.66 1.3
V10  0.75  0.10  0.01 0.63 0.37 1.0
V11  0.03  0.00  0.53 0.29 0.71 1.0
V12 -0.24  0.52  0.14 0.28 0.72 1.6

                       MR1  MR3  MR2
SS loadings           2.18 2.09 2.03
Proportion Var        0.18 0.17 0.17
Cumulative Var        0.18 0.36 0.53
Proportion Explained  0.35 0.33 0.32
Cumulative Proportion 0.35 0.68 1.00

 With factor correlations of 
     MR1  MR3  MR2
MR1 1.00 0.32 0.19
MR3 0.32 1.00 0.15
MR2 0.19 0.15 1.00

Basic Scree

Making a normal scree plot from there is quite simple. I just add this to my script:

scree(hwk2,
      pc=T,
      factors = F,
      main = "Scree Plot of Eigenvalues")

Which creates this:

What I Want

However, I want to graph simulated parallel analysis with it. In Jamovi this is super easy to accomplish:

However, I don't see an option for this so far. There is another version of scree I have tried fa.parallel but the legend comes out really strange:

fa.parallel(
  hwk2,
  n.obs = 93,
  fm = "minres",
  nfactors = 3,
  main = "Parallel Analysis Scree Plots",
  n.iter = 100,
  error.bars = FALSE,
  se.bars = FALSE,
  SMC = FALSE,
  ylabel = NULL,
  show.legend = F,
  sim = TRUE,
  quant = .95,
  use = "pairwise",
  plot = TRUE,
  correct = .5
)

I get either this if I remove the legend:

Or I get this annoying one with the legend:

Basically, I just need factor analysis and don't need principal components in the plot, but I can't figure out how to remove it.

I cannot reproduce your results from the sample data (e.g. V3 and V9 have 0 variance and there are only 6 observations). I think you can use the `fa` argument inside `fa.parallel` with `fa="pa"`, it's default is `fa="both"` which does both pc and pa. — Claudiu Papasteri, Feb 12 '22 at 11:55
I originally just used the head of my dput. Just added the full dput in my question. — Shawn Hemelstrand, Feb 12 '22 at 12:49

score 2 · Accepted Answer · answered Feb 12 '22 at 13:17

2

The only problem is that there are Heywood cases, so the fa analysis isn't trustworthy.

library(psych)

fa.parallel(
  hwk2,
  n.obs = 93,
  fa = "fa",                # you want only "fa", not "pc"
  show.legend = TRUE,       # show legend
  fm = "minres",
  nfactors = 3,
  main = "Parallel Analysis Scree Plots",
  n.iter = 100,
  error.bars = FALSE,
  se.bars = FALSE,
  SMC = FALSE,
  ylabel = NULL,
  sim = TRUE,
  quant = .95,
  use = "pairwise",
  plot = TRUE,
  correct = .5
)

answered Feb 12 '22 at 13:17

Claudiu Papasteri

2,469
1
17
30

Awesome. Thanks for the swift answer. Haven't been informed about what Heywood cases are but I'll educate myself with the link. – Shawn Hemelstrand Feb 12 '22 at 13:46
1

You are welcome. I didn't add to the Heywood case information because it is not relevant to the code, but being mindful of different estimation problems is important. Most of the time model specification (in fa, number of factors to extract and estimation method) and exclusion of problematic observed variables (e.g. very highly correlated) is enough to resolve such cases. – Claudiu Papasteri Feb 12 '22 at 14:00
Interesting. Thanks for the insight! – Shawn Hemelstrand Feb 12 '22 at 14:02

How can I create a parallel analysis scree in R?

Dataset Dput

EFA Before Scree Plot

Basic Scree

What I Want

1 Answers1