GLM EFA
Description
This shows the output of Exploratory Factor Analysis (EFA) functions from the package rwf. EFA uncovers the latent structure underlying a set of observed variables: it identifies groups of variables that share common variance, which are interpreted as factors. The functions include loading plots, scree plots, loading tables, residual diagnostics, and a full EFA report.
Installation instructions for rwf can be found here
The code can be found here
Plot Loadings
plot_loadings produces a heatmap of factor loadings,
making it easy to see which variables belong to which factor and how
strongly they load. Two matrix types are available:
- Structure matrix: correlations between each variable and each factor — shows total association including shared variance between factors.
- Pattern matrix: unique contribution of each factor to each variable after accounting for factor intercorrelations — use this to interpret factor meaning in oblique rotation.
Four datasets are shown to illustrate how the plot looks under different factor structures:
mtcars— a real dataset with a moderate two-factor structure.- A simulated dataset with two clear, well-separated factors (high within-factor correlations ~0.8, low cross-factor correlations ~0.1).
- A simulated dataset where only one factor has strong internal correlations — tests asymmetric structure detection.
- A near-identity correlation matrix (all correlations ~0.01) — no real factor structure; loadings should be near zero.
model<-psych::fa(mtcars,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
plot_loadings(model=model,matrix_type="structure")
cm<-matrix(c(1,.8,.8,.1,.1,.1,
.8,1,.8,.1,.1,.1,
.8,.8,1,.1,.1,.1,
.1,.1,.1,1,.8,.8,
.1,.1,.1,.8,1,.8,
.1,.1,.1,.8,.8,1),
ncol=6,nrow=6)
df1<-generate_correlation_matrix(cm,nrows=10000)
model1<-psych::fa(df1,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
plot_loadings(model=model1,matrix_type="pattern",base_size=10)
cm<-matrix(c(1,.1,.1,.1,.1,.1,
.1,1,.1,.1,.1,.1,
.1,.1,1,.1,.1,.1,
.1,.1,.1,1,.8,.8,
.1,.1,.1,.8,1,.8,
.1,.1,.1,.8,.8,1),
ncol=6,nrow=6)
df1<-generate_correlation_matrix(cm,nrows=10000)
model2<-psych::fa(df1,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
plot_loadings(model=model2,matrix_type="pattern",base_size=10)
cm<-matrix(c(1,.01,.01,.01,.01,.01,
.01,1,.01,.01,.01,.01,
.01,.01,1,.01,.01,.01,
.01,.01,.01,1,.01,.01,
.01,.01,.01,.01,1,.01,
.01,.01,.01,.01,.01,1),
ncol=6,nrow=6)
df1<-generate_correlation_matrix(cm,nrows=10000)
model3<-psych::fa(df1,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
plot_loadings(model=model3,matrix_type="pattern",base_size=10)
Plot Scree
A scree plot shows eigenvalues (the amount of variance each factor explains) in descending order. It is used to decide how many factors to retain. The conventional rule is to keep factors whose eigenvalue exceeds 1 (Kaiser criterion) or to look for the “elbow” — the point where the curve flattens out.
plot_scree plots eigenvalues alongside parallel analysis
results, which compares observed eigenvalues against those from random
data of the same size. Factors with eigenvalues above the parallel
analysis line are unlikely to be noise and should be retained.
Model Loadings
model_loadings returns the factor loading matrix as a
formatted table. The cut argument suppresses loadings below
a threshold so only meaningful loadings are shown — this makes factor
interpretation cleaner. The sort argument reorders
variables so that items loading on the same factor appear together.
Three calls demonstrate the options:
cut = NULL: all loadings shown, pattern matrix.cut = 0.4: only loadings ≥ 0.4 shown, structure matrix — a common threshold for “meaningful” loadings in applied research.cut = 0.4, matrix_type = "all": both pattern and structure matrices side by side, unsorted.
model<-psych::fa(mtcars,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
model_loadings(model=model,cut=NULL,matrix_type="pattern")## Matrix variable PA1 PA2 ## qsec Pattern qsec -0.92 ## hp Pattern hp 0.89 ## carb Pattern carb 0.85 ## vs Pattern vs -0.8 ## cyl Pattern cyl ## mpg Pattern mpg ## am Pattern am 0.93 ## gear Pattern gear 0.93 ## drat Pattern drat 0.78 ## wt Pattern wt ## disp Pattern disp
## Matrix variable PA1 PA2 ## hp Structure hp 0.93 ## cyl Structure cyl 0.85 -0.68 ## vs Structure vs -0.83 ## qsec Structure qsec -0.82 ## carb Structure carb 0.78 ## mpg Structure mpg -0.76 0.72 ## am Structure am 0.89 ## gear Structure gear 0.87 ## drat Structure drat 0.83 ## wt Structure wt 0.61 -0.81 ## disp Structure disp 0.75 -0.77
## Matrix variable PA1 PA2 ## 1 Pattern mpg -0.61 0.55 ## 2 Pattern cyl 0.71 -0.48 ## 3 Pattern disp 0.58 -0.6 ## 4 Pattern hp 0.89 ## 5 Pattern drat 0.78 ## 6 Pattern wt 0.42 -0.7 ## 7 Pattern qsec -0.92 ## 8 Pattern vs -0.8 ## 9 Pattern am 0.93 ## 10 Pattern gear 0.93 ## 11 Pattern carb 0.85 ## 12 Structure mpg -0.76 0.72 ## 13 Structure cyl 0.85 -0.68 ## 14 Structure disp 0.75 -0.77 ## 15 Structure hp 0.93 ## 16 Structure drat 0.83 ## 17 Structure wt 0.61 -0.81 ## 18 Structure qsec -0.82 ## 19 Structure vs -0.83 ## 20 Structure am 0.89 ## 21 Structure gear 0.87 ## 22 Structure carb 0.78
Compute Residual
Residuals in EFA are the differences between the observed correlation matrix and the correlation matrix reproduced by the factor model. A good-fitting model produces small residuals.
compute_residual_stats returns summary statistics of
these residuals (mean, RMSR — root mean square of residuals, and the
proportion of residuals above a given threshold). A low RMSR (typically
< 0.05) and few large residuals indicate that the factor model
adequately reproduces the observed correlations.
model<-psych::fa(mtcars,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
compute_residual_stats(model)## residual_statistics value critical formula ## 1 Root Mean Squared Residual 0.04419293 NA sqrt(mean(residuals^2)) ## 2 Number of absolute residuals > 0.05 13.00000000 NA abs(residuals)>0.05 ## 3 Proportion of absolute residuals > 0.05 0.23636364 0.5 numberLargeResiduals/nrow(residuals)
Report EFA
report_efa bundles all the above steps into a single
call — it fits the EFA model, extracts loadings, computes residual
statistics, and produces plots. It returns a structured list suitable
for inclusion in a report.
The example uses a simulated dataset with a known two-factor
structure (the same 6×6 correlation matrix used in the Plot Loadings
section), so the output should cleanly recover the two factors. The
computation chunk uses results='hide' to suppress
intermediate output; the formatted result is printed in the next
chunk.
cm<-matrix(c(1,.8,.8,.1,.1,.1,
.8,1,.8,.1,.1,.1,
.8,.8,1,.1,.1,.1,
.1,.1,.1,1,.8,.8,
.1,.1,.1,.8,1,.8,
.1,.1,.1,.8,.8,1),
ncol=6,nrow=6)
df1<-generate_correlation_matrix(cm,nrows=10000)
model<-psych::fa(df1,nfactors=2,rotate="oblimin",fm="pa",oblique.scores=TRUE)
result<-report_efa(model=model,df=df1)
## $correlations ## type X1 X2 X3 X4 X5 X6 ## X1 reproduced correlations 7.996617e-01 7.993266e-01 8.002435e-01 1.107260e-01 1.028915e-01 1.089642e-01 ## X2 reproduced correlations 7.993266e-01 7.989942e-01 7.999094e-01 1.092358e-01 1.014080e-01 1.074894e-01 ## X3 reproduced correlations 8.002435e-01 7.999094e-01 8.008264e-01 1.101502e-01 1.023116e-01 1.083938e-01 ## X4 reproduced correlations 1.107260e-01 1.092358e-01 1.101502e-01 8.101146e-01 8.071253e-01 8.018398e-01 ## X5 reproduced correlations 1.028915e-01 1.014080e-01 1.023116e-01 8.071253e-01 8.042173e-01 7.988870e-01 ## X6 reproduced correlations 1.089642e-01 1.074894e-01 1.083938e-01 8.018398e-01 7.988870e-01 7.936501e-01 ## X11 observed correlations 1.000000e+00 7.993694e-01 8.002989e-01 1.082151e-01 1.052364e-01 1.091376e-01 ## X21 observed correlations 7.993694e-01 1.000000e+00 7.999329e-01 1.104121e-01 1.007859e-01 1.069279e-01 ## X31 observed correlations 8.002989e-01 7.999329e-01 1.000000e+00 1.114853e-01 1.005890e-01 1.087810e-01 ## X41 observed correlations 1.082151e-01 1.104121e-01 1.114853e-01 1.000000e+00 8.072830e-01 8.018595e-01 ## X51 observed correlations 1.052364e-01 1.007859e-01 1.005890e-01 8.072830e-01 1.000000e+00 7.988307e-01 ## X61 observed correlations 1.091376e-01 1.069279e-01 1.087810e-01 8.018595e-01 7.988307e-01 1.000000e+00 ## X12 residual correlations 2.003383e-01 4.280987e-05 5.533587e-05 -2.510837e-03 2.344942e-03 1.733888e-04 ## X22 residual correlations 4.280987e-05 2.010058e-01 2.353943e-05 1.176290e-03 -6.221017e-04 -5.614650e-04 ## X32 residual correlations 5.533587e-05 2.353943e-05 1.991736e-01 1.335128e-03 -1.722532e-03 3.871522e-04 ## X42 residual correlations -2.510837e-03 1.176290e-03 1.335128e-03 1.898854e-01 1.577385e-04 1.971173e-05 ## X52 residual correlations 2.344942e-03 -6.221017e-04 -1.722532e-03 1.577385e-04 1.957827e-01 -5.631082e-05 ## X62 residual correlations 1.733888e-04 -5.614650e-04 3.871522e-04 1.971173e-05 -5.631082e-05 2.063499e-01 ## ## $npobs ## X1 X2 X3 X4 X5 X6 ## X1 10000 10000 10000 10000 10000 10000 ## X2 10000 10000 10000 10000 10000 10000 ## X3 10000 10000 10000 10000 10000 10000 ## X4 10000 10000 10000 10000 10000 10000 ## X5 10000 10000 10000 10000 10000 10000 ## X6 10000 10000 10000 10000 10000 10000 ## ## $residual_stats ## residual_statistics value critical formula ## 1 Root Mean Squared Residual 0.001121085 NA sqrt(mean(residuals^2)) ## 2 Number of absolute residuals > 0.05 0.000000000 NA abs(residuals)>0.05 ## 3 Proportion of absolute residuals > 0.05 0.000000000 0.5 numberLargeResiduals/nrow(residuals) ## ## $determinant_test ## determinant above_critical ## 1 0.01039058 TRUE ## ## $bartlett_test ## x_squared[bartlett] df[bartlett] p[bartlett] ## 1 45651.05 15 0 ## ## $kmo_test ## Overall_MSA MSA Kaiser_1974 ## X1 0.7694984 0.7693286 NA ## X2 0.7694984 0.7699709 NA ## X3 0.7694984 0.7686747 NA ## X4 0.7694984 0.7646410 NA ## X5 0.7694984 0.7678678 NA ## X6 0.7694984 0.7766697 NA ## ## $loadings ## Matrix variable PA1 PA2 type row.names.model.Vaccounted. ## 1 Pattern X4 0.9000000 0.000000## 2 Pattern X5 0.9000000 -0.010000 ## 3 Pattern X6 0.8900000 0.000000 ## 4 Pattern X1 0.0000000 0.890000 ## 5 Pattern X2 0.0000000 0.890000 ## 6 Pattern X3 0.0000000 0.890000 ## 7 Structure X4 0.9000000 0.120000 ## 8 Structure X5 0.9000000 0.110000 ## 9 Structure X6 0.8900000 0.120000 ## 10 Structure X1 0.1200000 0.890000 ## 11 Structure X2 0.1200000 0.890000 ## 12 Structure X3 0.1200000 0.890000 ## 13 2.4079385 2.399526 variance accounted SS loadings ## 14 0.4013231 0.399921 variance accounted Proportion Var ## 15 0.4013231 0.801244 variance accounted Cumulative Var ## 16 0.5008750 0.499125 variance accounted Proportion Explained ## 17 0.5008750 1.000000 variance accounted Cumulative Proportion ## ## $instruction_loading_critical_values ## sample critical_loading ## 1 50 0.75 ## 2 60 0.70 ## 3 70 0.65 ## 4 85 0.60 ## 5 100 0.55 ## 6 120 0.50 ## 7 150 0.45 ## 8 200 0.40 ## 9 250 0.35 ## 10 350 0.30 ## ## $weights ## PA1 PA2 ## X1 -0.04254933 0.34888255 ## X2 -0.04328924 0.34777669 ## X3 -0.04330714 0.35124134 ## X4 0.36358068 -0.04350210 ## X5 0.35229804 -0.04546112 ## X6 0.33157292 -0.03994150