Chapter 15:
Phil Wood
2025-08-10
Reasons for Drinking (Simulated) data
Load RFD Data
<a https://docs.google.com/document/d/17C7eM7Mc3pqBqe4eChuFqIIKSD7WVbhVoBfeUtA8fqM/edit?usp=drive_link
Plot shows regression lines, confidence interval for regression line and 50% confidence ellipse to communicate variation and slope
Read in Reasons for Drinking Data
# Load necessary libraries
#Read in data
AllRFDData <- read.table(
"RFDSim.txt",
sep="\t", header=TRUE)
myvars <- c(
"rridaw0", "rridbw0", "rridcw0", "rriddw0", "rridew0", "rridfw0", "rridgw0", "rridhw0", "rridiw0", "rridjw0", "rridkw0", "rridlw0", "rridmw0", "rridnw0", "rridow0", "rridpw0", "rridqw0", "rridrw0", "rridsw0", "rridtw0", "rriduw0", "rridvw0", "rridww0")
#subset variables of interest
RFD=AllRFDData[myvars]
#This is legacy code- original data had missing observations. Simulated data does not
CompleteRFD <- na.omit(RFD)
#Descriptive Statistics
library(psych)## Warning: package 'psych' was built under R version 4.5.1#descriptive statistics
describe(CompleteRFD)## vars n mean sd median trimmed mad min max range skew kurtosis
## rridaw0 1 2422 3.02 0.92 3.03 3.02 0.93 -0.24 6.36 6.60 -0.05 -0.03
## rridbw0 2 2422 2.94 0.92 2.95 2.94 0.93 -0.42 5.73 6.16 -0.06 -0.06
## rridcw0 3 2422 3.07 0.93 3.08 3.07 0.96 0.09 6.09 6.00 -0.04 -0.18
## rriddw0 4 2422 2.99 0.96 2.99 2.99 0.97 -0.23 6.20 6.43 0.01 -0.13
## rridew0 5 2422 3.12 0.88 3.14 3.13 0.85 -0.05 5.92 5.97 -0.09 0.11
## rridfw0 6 2422 1.90 1.04 1.93 1.91 1.03 -1.97 5.66 7.63 -0.08 0.00
## rridgw0 7 2422 1.76 0.98 1.76 1.76 0.97 -1.88 4.84 6.72 -0.04 -0.02
## rridhw0 8 2422 1.88 1.02 1.88 1.88 0.99 -1.90 5.14 7.04 -0.03 0.09
## rridiw0 9 2422 1.93 1.04 1.95 1.94 1.00 -1.79 5.42 7.21 -0.06 0.12
## rridjw0 10 2422 1.80 1.00 1.81 1.81 1.01 -1.99 5.12 7.10 -0.08 0.04
## rridkw0 11 2422 2.88 1.03 2.90 2.89 1.04 -0.63 6.31 6.95 -0.10 0.01
## rridlw0 12 2422 2.69 1.03 2.71 2.70 1.04 -0.67 5.73 6.39 -0.10 -0.15
## rridmw0 13 2422 1.39 0.67 1.42 1.39 0.69 -0.87 3.69 4.56 -0.08 -0.17
## rridnw0 14 2422 1.84 0.98 1.83 1.84 0.99 -2.26 5.44 7.70 -0.03 0.07
## rridow0 15 2422 2.57 1.03 2.58 2.56 1.08 -0.86 5.79 6.64 0.00 -0.19
## rridpw0 16 2422 2.19 1.12 2.21 2.20 1.12 -1.48 6.20 7.68 -0.09 0.01
## rridqw0 17 2422 2.77 1.01 2.76 2.77 1.02 -1.02 6.59 7.61 0.02 -0.03
## rridrw0 18 2422 3.10 0.94 3.10 3.10 0.95 -0.82 5.90 6.72 -0.06 -0.01
## rridsw0 19 2422 1.44 0.75 1.44 1.44 0.77 -1.10 3.76 4.86 0.01 -0.03
## rridtw0 20 2422 1.32 0.65 1.32 1.32 0.65 -1.32 3.35 4.67 0.02 0.01
## rriduw0 21 2422 1.30 0.64 1.30 1.30 0.64 -0.76 3.53 4.29 -0.07 -0.16
## rridvw0 22 2422 1.31 0.66 1.31 1.31 0.64 -0.88 3.32 4.20 -0.06 -0.03
## rridww0 23 2422 1.47 0.79 1.47 1.48 0.82 -1.55 3.98 5.53 -0.07 -0.04
## se
## rridaw0 0.02
## rridbw0 0.02
## rridcw0 0.02
## rriddw0 0.02
## rridew0 0.02
## rridfw0 0.02
## rridgw0 0.02
## rridhw0 0.02
## rridiw0 0.02
## rridjw0 0.02
## rridkw0 0.02
## rridlw0 0.02
## rridmw0 0.01
## rridnw0 0.02
## rridow0 0.02
## rridpw0 0.02
## rridqw0 0.02
## rridrw0 0.02
## rridsw0 0.02
## rridtw0 0.01
## rriduw0 0.01
## rridvw0 0.01
## rridww0 0.02Example of an Exploratory Factor analysis Setup
These two approaches yield the same values although in practice minor differences may emerge between the two approaches having to do with when the rotation is specified relative to extraction. In this analysis F1, F2, F3, and F4 from fa are represented as ML1, ML3, ML4, and ML2 respectively.
Using fa()
# Packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(psych) # fa()
library(GPArotation) # geominQ / geominT rotations##
## Attaching package: 'GPArotation'## The following objects are masked from 'package:psych':
##
## equamax, varimin# EFA: 4 factors, ML extraction, Geomin (oblique)
fit <- fa(CompleteRFD, nfactors = 4, rotate = "geominQ", fm = "ml", n.obs = N)
# 4) Inspect
print(fit$loadings, cutoff = .0, sort = TRUE) # pattern matrix (sorted)##
## Loadings:
## ML2 ML3 ML1 ML4
## rridsw0 0.720 -0.007 0.001 0.037
## rridtw0 0.818 -0.013 0.004 -0.048
## rriduw0 0.901 -0.003 0.005 0.006
## rridvw0 0.906 0.026 -0.007 0.017
## rridww0 0.783 -0.004 0.064 -0.030
## rridfw0 -0.024 0.884 0.036 -0.028
## rridgw0 -0.004 0.917 -0.021 -0.027
## rridhw0 -0.018 0.824 0.015 0.069
## rridjw0 0.020 0.894 -0.005 -0.010
## rridaw0 0.019 0.027 0.766 0.074
## rridbw0 0.118 0.014 0.681 -0.018
## rridcw0 -0.003 -0.002 0.952 -0.012
## rriddw0 -0.001 0.009 0.893 0.027
## rridew0 -0.042 0.066 0.543 0.130
## rridkw0 -0.057 0.043 0.084 0.772
## rridow0 0.079 0.013 0.118 0.659
## rridpw0 0.070 0.072 -0.035 0.690
## rridqw0 -0.006 -0.018 0.007 0.887
## rridrw0 -0.023 -0.060 0.288 0.645
## rridiw0 0.125 0.474 0.140 0.135
## rridlw0 -0.033 0.012 -0.022 0.402
## rridmw0 0.301 0.208 -0.091 0.159
## rridnw0 0.145 0.067 -0.103 0.252
##
## ML2 ML3 ML1 ML4
## SS loadings 3.595 3.391 3.199 3.013
## Proportion Var 0.156 0.147 0.139 0.131
## Cumulative Var 0.156 0.304 0.443 0.574fit$Phi # factor correlations## ML2 ML3 ML1 ML4
## ML2 1.0000000 0.3245435 0.1735281 0.1780090
## ML3 0.3245435 1.0000000 0.3994886 0.4786475
## ML1 0.1735281 0.3994886 1.0000000 0.7094108
## ML4 0.1780090 0.4786475 0.7094108 1.0000000fit$communality # communalities## rridaw0 rridbw0 rridcw0 rriddw0 rridew0 rridfw0 rridgw0
## 0.69815462 0.49683285 0.88711963 0.83909097 0.44340620 0.77016990 0.80139371
## rridhw0 rridiw0 rridjw0 rridkw0 rridlw0 rridmw0 rridnw0
## 0.73908529 0.46940354 0.79996012 0.71569284 0.15038265 0.21161744 0.08743975
## rridow0 rridpw0 rridqw0 rridrw0 rridsw0 rridtw0 rriduw0
## 0.59725654 0.51762742 0.77948583 0.70828875 0.52612973 0.65253282 0.81432887
## rridvw0 rridww0
## 0.84077930 0.62211295fit$uniquenesses # uniquenesses## rridaw0 rridbw0 rridcw0 rriddw0 rridew0 rridfw0 rridgw0 rridhw0
## 0.3018454 0.5031672 0.1128804 0.1609090 0.5565938 0.2298301 0.1986063 0.2609147
## rridiw0 rridjw0 rridkw0 rridlw0 rridmw0 rridnw0 rridow0 rridpw0
## 0.5305965 0.2000399 0.2843072 0.8496174 0.7883826 0.9125602 0.4027435 0.4823726
## rridqw0 rridrw0 rridsw0 rridtw0 rriduw0 rridvw0 rridww0
## 0.2205142 0.2917113 0.4738703 0.3474672 0.1856711 0.1592207 0.3778871fit$STATISTIC; fit$dof; fit$PVAL # chi-square test (ML)## [1] 2332.49## [1] 167## [1] 0Using lavaan
library(lavaan)## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covvars <- c(
"rridaw0","rridbw0","rridcw0","rriddw0","rridew0","rridfw0","rridgw0","rridhw0",
"rridiw0","rridjw0","rridkw0","rridlw0","rridmw0","rridnw0","rridow0","rridpw0",
"rridqw0","rridrw0","rridsw0","rridtw0","rriduw0","rridvw0","rridww0"
)
model <- paste0('
# 4 exploratory factors in the same block ("efa4")
efa("efa4")*F1 + efa("efa4")*F2 + efa("efa4")*F3 + efa("efa4")*F4 =~
', paste(vars, collapse = " + "), '
')
fit_lav <- cfa(
model,
data = CompleteRFD,
std.lv = TRUE, # set factor variances = 1
missing = "fiml",
estimator = "ml",
rotation = "geomin", # <-- rotation name goes here
rotation.args = list(
geomin.epsilon = 0.01, # typical epsilon
orthogonal = FALSE # oblique rotation (factors can correlate)
)
)
summary(fit_lav, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 144
## Row rank of the constraints matrix 12
##
## Rotation method GEOMIN OBLIQUE
## Geomin epsilon 0.01
## Rotation algorithm (rstarts) GPA (30)
## Standardized metric TRUE
## Row weights None
##
## Number of observations 2422
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 2344.266
## Degrees of freedom 167
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 40456.892
## Degrees of freedom 253
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.946
## Tucker-Lewis Index (TLI) 0.918
##
## Robust Comparative Fit Index (CFI) 0.946
## Robust Tucker-Lewis Index (TLI) 0.918
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -54129.245
## Loglikelihood unrestricted model (H1) -52957.112
##
## Akaike (AIC) 108522.490
## Bayesian (BIC) 109287.080
## Sample-size adjusted Bayesian (SABIC) 108867.686
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.073
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.076
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 0.000
##
## Robust RMSEA 0.076
## 90 Percent confidence interval - lower 0.073
## 90 Percent confidence interval - upper 0.079
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 0.013
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.034
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~ efa4
## rridaw0 0.708 0.020 36.007 0.000 0.708 0.766
## rridbw0 0.625 0.022 28.931 0.000 0.625 0.681
## rridcw0 0.883 0.017 52.488 0.000 0.883 0.952
## rriddw0 0.858 0.018 47.849 0.000 0.858 0.893
## rridew0 0.479 0.023 20.560 0.000 0.479 0.543
## rridfw0 0.037 0.015 2.548 0.011 0.037 0.036
## rridgw0 -0.021 0.013 -1.652 0.098 -0.021 -0.021
## rridhw0 0.015 0.014 1.076 0.282 0.015 0.015
## rridiw0 0.145 0.026 5.570 0.000 0.145 0.140
## rridjw0 -0.005 0.013 -0.414 0.679 -0.005 -0.005
## rridkw0 0.086 0.021 4.036 0.000 0.086 0.084
## rridlw0 -0.023 0.030 -0.777 0.437 -0.023 -0.022
## rridmw0 -0.061 0.019 -3.214 0.001 -0.061 -0.091
## rridnw0 -0.100 0.030 -3.318 0.001 -0.100 -0.103
## rridow0 0.122 0.025 4.920 0.000 0.122 0.118
## rridpw0 -0.040 0.021 -1.901 0.057 -0.040 -0.035
## rridqw0 0.007 0.014 0.527 0.598 0.007 0.007
## rridrw0 0.270 0.021 12.955 0.000 0.270 0.288
## rridsw0 0.001 0.015 0.054 0.957 0.001 0.001
## rridtw0 0.002 0.011 0.218 0.827 0.002 0.004
## rriduw0 0.003 0.008 0.378 0.706 0.003 0.005
## rridvw0 -0.005 0.008 -0.575 0.566 -0.005 -0.007
## rridww0 0.050 0.015 3.257 0.001 0.050 0.064
## F2 =~ efa4
## rridaw0 0.025 0.012 2.187 0.029 0.025 0.027
## rridbw0 0.012 0.014 0.866 0.386 0.012 0.014
## rridcw0 -0.002 0.008 -0.184 0.854 -0.002 -0.002
## rriddw0 0.008 0.009 0.905 0.365 0.008 0.009
## rridew0 0.058 0.017 3.491 0.000 0.058 0.066
## rridfw0 0.917 0.018 50.076 0.000 0.917 0.884
## rridgw0 0.899 0.017 52.080 0.000 0.899 0.917
## rridhw0 0.837 0.018 46.142 0.000 0.837 0.824
## rridiw0 0.491 0.021 23.270 0.000 0.491 0.474
## rridjw0 0.893 0.017 51.715 0.000 0.893 0.894
## rridkw0 0.044 0.013 3.345 0.001 0.044 0.043
## rridlw0 0.013 0.024 0.541 0.588 0.013 0.012
## rridmw0 0.140 0.016 8.873 0.000 0.140 0.208
## rridnw0 0.065 0.024 2.709 0.007 0.065 0.067
## rridow0 0.014 0.014 1.009 0.313 0.014 0.013
## rridpw0 0.080 0.020 3.986 0.000 0.080 0.072
## rridqw0 -0.018 0.012 -1.565 0.118 -0.018 -0.018
## rridrw0 -0.056 0.012 -4.539 0.000 -0.056 -0.060
## rridsw0 -0.005 0.012 -0.418 0.676 -0.005 -0.007
## rridtw0 -0.008 0.009 -0.959 0.337 -0.008 -0.013
## rriduw0 -0.002 0.007 -0.259 0.796 -0.002 -0.003
## rridvw0 0.017 0.007 2.489 0.013 0.017 0.026
## rridww0 -0.003 0.011 -0.295 0.768 -0.003 -0.004
## F3 =~ efa4
## rridaw0 0.069 0.018 3.784 0.000 0.069 0.074
## rridbw0 -0.016 0.019 -0.870 0.384 -0.016 -0.018
## rridcw0 -0.012 0.011 -1.015 0.310 -0.012 -0.012
## rriddw0 0.026 0.013 1.947 0.051 0.026 0.027
## rridew0 0.115 0.024 4.725 0.000 0.115 0.130
## rridfw0 -0.029 0.014 -2.064 0.039 -0.029 -0.028
## rridgw0 -0.026 0.013 -1.955 0.051 -0.026 -0.027
## rridhw0 0.070 0.018 3.833 0.000 0.070 0.069
## rridiw0 0.140 0.028 5.005 0.000 0.140 0.135
## rridjw0 -0.010 0.013 -0.802 0.423 -0.010 -0.010
## rridkw0 0.795 0.024 33.162 0.000 0.795 0.772
## rridlw0 0.416 0.033 12.502 0.000 0.416 0.402
## rridmw0 0.107 0.021 5.033 0.000 0.107 0.159
## rridnw0 0.246 0.033 7.406 0.000 0.246 0.252
## rridow0 0.680 0.027 25.663 0.000 0.680 0.659
## rridpw0 0.772 0.028 27.684 0.000 0.772 0.690
## rridqw0 0.898 0.021 43.150 0.000 0.898 0.887
## rridrw0 0.606 0.022 27.269 0.000 0.606 0.645
## rridsw0 0.028 0.017 1.678 0.093 0.028 0.037
## rridtw0 -0.032 0.013 -2.512 0.012 -0.032 -0.048
## rriduw0 0.004 0.009 0.454 0.650 0.004 0.006
## rridvw0 0.011 0.008 1.295 0.195 0.011 0.017
## rridww0 -0.023 0.015 -1.590 0.112 -0.023 -0.030
## F4 =~ efa4
## rridaw0 0.018 0.010 1.813 0.070 0.018 0.019
## rridbw0 0.108 0.015 7.121 0.000 0.108 0.118
## rridcw0 -0.003 0.007 -0.400 0.689 -0.003 -0.003
## rriddw0 -0.001 0.008 -0.130 0.896 -0.001 -0.001
## rridew0 -0.037 0.014 -2.723 0.006 -0.037 -0.042
## rridfw0 -0.025 0.010 -2.500 0.012 -0.025 -0.024
## rridgw0 -0.004 0.009 -0.488 0.626 -0.004 -0.004
## rridhw0 -0.018 0.010 -1.860 0.063 -0.018 -0.018
## rridiw0 0.130 0.018 7.356 0.000 0.130 0.125
## rridjw0 0.020 0.009 2.233 0.026 0.020 0.020
## rridkw0 -0.059 0.012 -4.849 0.000 -0.059 -0.057
## rridlw0 -0.035 0.021 -1.656 0.098 -0.035 -0.033
## rridmw0 0.202 0.014 14.541 0.000 0.202 0.301
## rridnw0 0.142 0.021 6.672 0.000 0.142 0.145
## rridow0 0.082 0.015 5.424 0.000 0.082 0.079
## rridpw0 0.078 0.018 4.438 0.000 0.078 0.070
## rridqw0 -0.006 0.010 -0.652 0.514 -0.006 -0.006
## rridrw0 -0.022 0.008 -2.610 0.009 -0.022 -0.023
## rridsw0 0.540 0.014 38.593 0.000 0.540 0.720
## rridtw0 0.535 0.012 46.053 0.000 0.535 0.818
## rriduw0 0.575 0.010 55.294 0.000 0.575 0.901
## rridvw0 0.596 0.011 56.462 0.000 0.596 0.906
## rridww0 0.619 0.014 44.002 0.000 0.619 0.783
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.399 0.018 22.319 0.000 0.399 0.399
## F3 0.709 0.013 53.564 0.000 0.709 0.709
## F4 0.174 0.020 8.609 0.000 0.174 0.174
## F2 ~~
## F3 0.479 0.017 28.321 0.000 0.479 0.479
## F4 0.325 0.019 17.306 0.000 0.325 0.325
## F3 ~~
## F4 0.178 0.020 8.791 0.000 0.178 0.178
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw0 3.023 0.019 160.998 0.000 3.023 3.271
## .rridbw0 2.935 0.019 157.441 0.000 2.935 3.199
## .rridcw0 3.067 0.019 162.812 0.000 3.067 3.308
## .rriddw0 2.995 0.020 153.421 0.000 2.995 3.117
## .rridew0 3.122 0.018 174.332 0.000 3.122 3.542
## .rridfw0 1.904 0.021 90.315 0.000 1.904 1.835
## .rridgw0 1.761 0.020 88.422 0.000 1.761 1.797
## .rridhw0 1.878 0.021 90.961 0.000 1.878 1.848
## .rridiw0 1.935 0.021 91.941 0.000 1.935 1.868
## .rridjw0 1.801 0.020 88.784 0.000 1.801 1.804
## .rridkw0 2.881 0.021 137.706 0.000 2.881 2.798
## .rridlw0 2.692 0.021 128.059 0.000 2.692 2.602
## .rridmw0 1.389 0.014 101.651 0.000 1.389 2.065
## .rridnw0 1.839 0.020 92.487 0.000 1.839 1.879
## .rridow0 2.566 0.021 122.358 0.000 2.566 2.486
## .rridpw0 2.187 0.023 96.154 0.000 2.187 1.954
## .rridqw0 2.767 0.021 134.527 0.000 2.767 2.734
## .rridrw0 3.098 0.019 162.275 0.000 3.098 3.297
## .rridsw0 1.441 0.015 94.437 0.000 1.441 1.919
## .rridtw0 1.320 0.013 99.373 0.000 1.320 2.019
## .rriduw0 1.296 0.013 99.917 0.000 1.296 2.030
## .rridvw0 1.310 0.013 98.102 0.000 1.310 1.993
## .rridww0 1.469 0.016 91.442 0.000 1.469 1.858
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridaw0 0.258 0.009 30.128 0.000 0.258 0.302
## .rridbw0 0.424 0.013 32.506 0.000 0.424 0.503
## .rridcw0 0.097 0.006 16.854 0.000 0.097 0.113
## .rriddw0 0.148 0.006 23.283 0.000 0.148 0.161
## .rridew0 0.432 0.013 33.516 0.000 0.432 0.557
## .rridfw0 0.247 0.009 26.184 0.000 0.247 0.230
## .rridgw0 0.191 0.008 23.653 0.000 0.191 0.199
## .rridhw0 0.269 0.010 28.098 0.000 0.269 0.261
## .rridiw0 0.569 0.017 33.553 0.000 0.569 0.531
## .rridjw0 0.199 0.008 24.248 0.000 0.199 0.200
## .rridkw0 0.301 0.011 26.797 0.000 0.301 0.284
## .rridlw0 0.909 0.027 34.003 0.000 0.909 0.850
## .rridmw0 0.357 0.010 34.212 0.000 0.357 0.788
## .rridnw0 0.874 0.025 34.375 0.000 0.874 0.913
## .rridow0 0.429 0.014 30.457 0.000 0.429 0.403
## .rridpw0 0.604 0.020 30.938 0.000 0.604 0.482
## .rridqw0 0.226 0.011 21.107 0.000 0.226 0.221
## .rridrw0 0.257 0.009 29.042 0.000 0.257 0.292
## .rridsw0 0.267 0.009 31.418 0.000 0.267 0.474
## .rridtw0 0.148 0.005 28.844 0.000 0.148 0.347
## .rriduw0 0.076 0.003 23.728 0.000 0.076 0.186
## .rridvw0 0.069 0.003 21.173 0.000 0.069 0.159
## .rridww0 0.236 0.008 30.423 0.000 0.236 0.378
## F1 1.000 1.000 1.000
## F2 1.000 1.000 1.000
## F3 1.000 1.000 1.000
## F4 1.000 1.000 1.000lavInspect(fit_lav, "cor.lv") # factor correlations## F1 F2 F3 F4
## F1 1.000
## F2 0.399 1.000
## F3 0.709 0.479 1.000
## F4 0.174 0.325 0.178 1.000library(dplyr)
library(tidyr)
# Extract standardized factor loadings
loadings_df <- standardizedSolution(fit_lav) %>%
filter(op == "=~") %>% # keep only loadings
select(Factor = lhs, Variable = rhs, Loading = est.std)
# Sort within each factor by absolute value of loading
loadings_sorted <- loadings_df %>%
arrange(Factor, desc(abs(Loading)))
# Make wide table: variables × factors
loadings_wide <- loadings_sorted %>%
pivot_wider(names_from = Factor, values_from = Loading)
library(knitr)
loadings_wide %>%
mutate(across(-Variable, ~ round(.x, 3))) %>%
kable()| Variable | F1 | F2 | F3 | F4 |
|---|---|---|---|---|
| rridcw0 | 0.952 | -0.002 | -0.012 | -0.003 |
| rriddw0 | 0.893 | 0.009 | 0.027 | -0.001 |
| rridaw0 | 0.766 | 0.027 | 0.074 | 0.019 |
| rridbw0 | 0.681 | 0.014 | -0.018 | 0.118 |
| rridew0 | 0.543 | 0.066 | 0.130 | -0.042 |
| rridrw0 | 0.288 | -0.060 | 0.645 | -0.023 |
| rridiw0 | 0.140 | 0.474 | 0.135 | 0.125 |
| rridow0 | 0.118 | 0.013 | 0.659 | 0.079 |
| rridnw0 | -0.103 | 0.067 | 0.252 | 0.145 |
| rridmw0 | -0.091 | 0.208 | 0.159 | 0.301 |
| rridkw0 | 0.084 | 0.043 | 0.772 | -0.057 |
| rridww0 | 0.064 | -0.004 | -0.030 | 0.783 |
| rridfw0 | 0.036 | 0.884 | -0.028 | -0.024 |
| rridpw0 | -0.035 | 0.072 | 0.690 | 0.070 |
| rridlw0 | -0.022 | 0.012 | 0.402 | -0.033 |
| rridgw0 | -0.021 | 0.917 | -0.027 | -0.004 |
| rridhw0 | 0.015 | 0.824 | 0.069 | -0.018 |
| rridqw0 | 0.007 | -0.018 | 0.887 | -0.006 |
| rridvw0 | -0.007 | 0.026 | 0.017 | 0.906 |
| rridjw0 | -0.005 | 0.894 | -0.010 | 0.020 |
| rriduw0 | 0.005 | -0.003 | 0.006 | 0.901 |
| rridtw0 | 0.004 | -0.013 | -0.048 | 0.818 |
| rridsw0 | 0.001 | -0.007 | 0.037 | 0.720 |
write.table(loadings_wide, "clipboard", sep="\t", row.names = FALSE)Heatmaps of Correlation Matrices (Figures 15.1 and 15.3)
Simple Heatmaps using corrplot
library(corrplot)## corrplot 0.95 loaded# Simple Correlation heatmap plot using corrplot
library(corrplot)
m <- cor(CompleteRFD, use = "pairwise.complete.obs")
corrplot(m, method = "color", type = "upper",
col = colorRampPalette(c("navy","white","firebrick"))(200),
addCoef.col = "black", # show correlation values
tl.col = "black", tl.srt = 90,tl.cex=.75,number.cex=.5)
Clustered Heatmap using Corrplot
# Two methods commonly used are complete and Ward's (modified) clustering
# These show slightly different patterns for these data (not discussed in book)
corrplot(m, method = "color", type = "upper",
col = colorRampPalette(c("navy","white","firebrick"))(200),
addCoef.col = "black", # show correlation values
tl.col = "black", tl.srt = 90,tl.cex=.75,number.cex=.5,
order = "hclust",
hclust.method = "complete"
)
corrplot(m, method = "color", type = "upper",
col = colorRampPalette(c("navy","white","firebrick"))(200),
addCoef.col = "black", # show correlation values
tl.col = "black", tl.srt = 90,tl.cex=.75,number.cex=.5,
order = "hclust",
hclust.method = "ward.D2"
)
Clustering information in addition to correlations using pheatmap
Useful for thinking about distinct clusters within item group (Not discussed in book)
library(pheatmap)## Warning: package 'pheatmap' was built under R version 4.5.1pheatmap(m,
color = colorRampPalette(c("navy","white","firebrick"))(200),
border_color = NA,
cluster_rows = TRUE,
cluster_cols = TRUE)
Heatmaps using ggplot
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(dplyr)
library(tidyr)
# m: your square, named correlation matrix
vars <- colnames(m)
# Build long data + keep upper triangle
corr_df <- expand.grid(var1 = vars, var2 = vars, stringsAsFactors = FALSE) %>%
mutate(
i = match(var1, vars),
j = match(var2, vars),
correlation = m[cbind(i, j)]
) %>%
filter(i >= j) %>% # keep upper triangle (incl. diagonal)
mutate(label = sprintf("%.2f", correlation))
ggplot(corr_df, aes(var1, var2, fill = correlation)) +
geom_tile(color = "white") +
geom_text(aes(label = label), size = 2) +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
scale_x_discrete(limits = vars) + # left -> right in original order
scale_y_discrete(limits = rev(vars)) + # TOP -> BOTTOM reversed so diag is TL->BR
coord_equal() +
labs(fill = "r", x = NULL, y = NULL) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
)
Hierarchical Clustering with ggplot
library(ggplot2)
library(dplyr)
library(tidyr)
# m: square correlation matrix with row/col names
vars <- colnames(m)
# 1) Cluster variables (distance = 1 - |r|; change if you want sign-sensitive)
dist_mat <- as.dist(1 - abs(m))
hc <- hclust(dist_mat, method = "complete")
ord <- vars[hc$order] # clustered order
# 2) Build long data in clustered order and keep the UPPER triangle
corr_df <- expand.grid(var1 = ord, var2 = ord, stringsAsFactors = FALSE) %>%
mutate(
i = match(var1, ord),
j = match(var2, ord),
# index into the ORIGINAL matrix by positions in 'vars'
ii = match(var1, vars),
jj = match(var2, vars),
correlation = m[cbind(ii, jj)]
) %>%
filter(j >= i) %>% # keep upper triangle (incl. diagonal)
mutate(label = sprintf("%.2f", correlation))
# 3) Plot (diag runs top-left -> bottom-right)
ggplot(corr_df, aes(var1, var2, fill = correlation)) +
geom_tile(color = "white") +
geom_text(aes(label = label), size = 2) +
scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0) +
scale_x_discrete(limits = ord) +
scale_y_discrete(limits = rev(ord)) + # flip y so it looks like corrplot(type="upper")
coord_equal() +
labs(fill = "r", x = NULL, y = NULL) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
)
Example Heatmap for Longitudinal Data (Figure 15.2)
Read in Reasons for Drinking Scale Dataset (zipped)
<a https://drive.google.com/file/d/1UjlcmACcAxOt4GRU91hqJdMKswNPiR54/view?usp=sharing
Simple Corrplot
vars <- c(
"enhw0", "enhw1", "enhw3", "enhw5", "enhw7",
"socw0", "socw1", "socw3", "socw5", "socw7",
"cnfw0", "cnfw1", "cnfw3", "cnfw5", "cnfw7",
"copw0", "copw1", "copw3", "copw5", "copw7",
"tstw0", "tstw1", "tstw3", "tstw5", "tstw7"
)
RFDScales <- read.table(
"RFDScales.txt",
sep="\t", header=TRUE)
RFDSubset=RFDScales[vars]
m <- cor(RFDSubset, use = "pairwise.complete.obs")
corrplot(m, method = "color", type = "upper",
col = colorRampPalette(c("navy","white","firebrick"))(200),
addCoef.col = "black", # show correlation values
tl.col = "black", tl.srt = 90,tl.cex=.75,number.cex=.5)
Pairplots of individual scales (Figures 15.4 and 15.5)
This plot will not resemble Figures 15.4 and 15.5 because the data are simulated. The regression lines will therefore look straight and the likert-format of the questions is not evident. Reversing the order of the questions (from rridew0 to rridaw0) will produce Figure 15.5
# packages
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(ggplot2)
library(dplyr)
# choose your data frame here (replace `CompleteRFD` if needed)
df <- CompleteRFD %>%
select(all_of(c("rridaw0","rridbw0","rridcw0","rriddw0","rridew0"))) %>%
mutate(across(everything(), as.numeric)) # ensure numeric
# custom continuous panel: jittered points + loess + 75% ellipse
cont_panel <- function(data, mapping, ...) {
ggplot(data = data, mapping = mapping) +
geom_point(
position = position_jitter(width = 0.15, height = 0.15),
alpha = 0.6, size = 1.3, na.rm = TRUE
) +
stat_smooth(method = "loess", se = FALSE, linewidth = 0.7, na.rm = TRUE) +
stat_ellipse(level = 0.75, type = "norm", linewidth = 0.7, na.rm = TRUE) +
theme_minimal(base_size = 10)
}
# build the pairplot
GGally::ggpairs(
df,
lower = list(continuous = cont_panel),
upper = list(continuous = GGally::wrap(GGally::ggally_cor, method = "pearson", size = 4)),
diag = list(continuous = "densityDiag")
) +
theme_minimal(base_size = 11) +
theme(
panel.grid = element_blank(),
strip.text = element_text(face = "bold")
)## `geom_smooth()` using formula = 'y ~ x'## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Correlations of Items for Two Scales of Reasons for Drinking (Tables 15.2 & 15.6)
library(psych)
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.5.1##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsvars <- c("rridcw0", "rriddw0", "rridgw0", "rridjw0",
"rridkw0", "rridqw0", "rriduw0", "rridvw0")
RFD=AllRFDData[vars]
apa_cor_table <- function(data, vars = NULL, digits = 2,
show_n = TRUE, star_levels = c(.001, .01, .05)) {
if (!is.null(vars)) data <- dplyr::select(data, all_of(vars))
data <- mutate(across(where(is.numeric), as.numeric), .data = data)
ct <- psych::corr.test(data) # r, p, n
r <- ct$r
p <- ct$p
n <- ct$n
# Significance stars
star_fun <- function(pv) {
ifelse(is.na(pv), "",
ifelse(pv < star_levels[1], "***",
ifelse(pv < star_levels[2], "**",
ifelse(pv < star_levels[3], "*", ""))))
}
stars <- matrix(star_fun(p), nrow = nrow(p), dimnames = dimnames(p))
# Format r to two decimals (exactly)
r_fmt <- formatC(r, format = "f", digits = digits)
# Keep only lower triangle (APA style); blank upper triangle
lower_mask <- row(r) > col(r)
upper_mask <- row(r) < col(r)
diag(r_fmt) <- formatC(1, format = "f", digits = digits) # APA often shows 1.00
r_fmt[upper_mask] <- "" # blank upper triangle
# append stars to lower triangle cells
r_fmt[lower_mask] <- paste0(r_fmt[lower_mask], stars[lower_mask])
# Build table
out <- as.data.frame(r_fmt, stringsAsFactors = FALSE)
out <- tibble::rownames_to_column(out, "Variable")
# Optional N row at bottom
if (show_n) {
n_row <- c("N", rep(as.character(n), ncol(r)))
out <- rbind(out, setNames(as.list(n_row), names(out)))
}
out
}
tab <- apa_cor_table(RFD, vars = vars, digits = 2, show_n = TRUE)
kable(tab, align = "c", caption = "APA-Style Correlation Matrix (lower triangle)") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = FALSE) %>%
add_footnote(label = c(
"* p < .05, ** p < .01, *** p < .001"
), notation = "symbol")| Variable | rridcw0 | rriddw0 | rridgw0 | rridjw0 | rridkw0 | rridqw0 | rriduw0 | rridvw0 |
|---|---|---|---|---|---|---|---|---|
| rridcw0 | 1.00 | |||||||
| rriddw0 | 0.87*** | 1.00 | ||||||
| rridgw0 | 0.30*** | 0.31*** | 1.00 | |||||
| rridjw0 | 0.33*** | 0.33*** | 0.78*** | 1.00 | ||||
| rridkw0 | 0.60*** | 0.60*** | 0.35*** | 0.38*** | 1.00 | |||
| rridqw0 | 0.59*** | 0.59*** | 0.34*** | 0.35*** | 0.75*** | 1.00 | ||
| rriduw0 | 0.15*** | 0.16*** | 0.26*** | 0.28*** | 0.11*** | 0.14*** | 1.00 | |
| rridvw0 | 0.16*** | 0.16*** | 0.29*** | 0.31*** | 0.12*** | 0.16*** | 0.85*** | 1.00 |
| N | 2422 | 2422 | 2422 | 2422 | 2422 | 2422 | 2422 | 2422 |
| p < .05, ** p < .01, *** p < .001 |
Eigenvectors (Table 15.3)
# install.packages(c("knitr","kableExtra","dplyr")) # if not installed
library(dplyr)
library(knitr)
library(kableExtra)
# --- 1. Calculate correlation matrix and eigen decomposition ---
RFD_cor <- cor(RFD, use = "pairwise.complete.obs") # correlation matrix
eig <- eigen(RFD_cor)
# --- 2. Extract eigenvectors ---
eigenvectors <- eig$vectors
colnames(eigenvectors) <- paste0("Factor ", 1:ncol(eigenvectors))
rownames(eigenvectors) <- colnames(RFD)
# --- 3. Format to two decimals ---
eig_df <- as.data.frame(eigenvectors) %>%
mutate(across(everything(), ~formatC(.x, digits = 2, format = "f"))) %>%
tibble::rownames_to_column("Variable")
# --- 4. Create APA-style table ---
kable(eig_df, align = c("l", rep("c", ncol(eig_df)-1)),
caption = "Table 15.3 Eigenvectors of Correlation Matrix in Figure 15.2
Eigenvectors from Principal Components Analysis of the RFD Dataset") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
column_spec(1, bold = TRUE) %>%
add_header_above(c(" " = 1, "Eigenvectors" = ncol(eig_df)-1))| Variable | Factor 1 | Factor 2 | Factor 3 | Factor 4 | Factor 5 | Factor 6 | Factor 7 | Factor 8 |
|---|---|---|---|---|---|---|---|---|
| rridcw0 | -0.41 | 0.25 | -0.24 | -0.46 | 0.01 | -0.02 | -0.10 | 0.70 |
| rriddw0 | -0.41 | 0.25 | -0.24 | -0.46 | 0.02 | 0.02 | 0.09 | -0.70 |
| rridgw0 | -0.34 | -0.19 | 0.59 | -0.10 | 0.24 | 0.65 | 0.05 | 0.03 |
| rridjw0 | -0.35 | -0.19 | 0.56 | -0.09 | -0.21 | -0.69 | -0.03 | -0.02 |
| rridkw0 | -0.40 | 0.25 | -0.05 | 0.50 | -0.68 | 0.24 | -0.01 | 0.00 |
| rridqw0 | -0.40 | 0.23 | -0.10 | 0.55 | 0.66 | -0.21 | 0.03 | 0.00 |
| rriduw0 | -0.22 | -0.59 | -0.33 | 0.04 | -0.05 | -0.02 | 0.69 | 0.10 |
| rridvw0 | -0.23 | -0.59 | -0.31 | 0.06 | 0.02 | 0.05 | -0.70 | -0.09 |
eig <- eigen(RFD_cor)Eigenvalues (Table 15.4)
# 1. Eigenvalues
eigen_df <- data.frame(
Factor = paste0("Factor ", 1:length(eig$values)),
Eigenvalue = eig$values,
PercentVariance = eig$values / sum(eig$values) * 100,
CumulativeVariance = cumsum(eig$values) / sum(eig$values) * 100
) %>%
mutate(across(-Factor, ~round(.x, 2))) # APA: two decimals
# 2. Format as APA-style table
kable(eigen_df,
caption = "Table 15.4
Eigenvalues and Variance Explained from Principal Components Analysis of the RFD Dataset",
align = c("l", rep("c", 3))) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
add_header_above(c(" " = 1, "Eigenvalues" = 3))| Factor | Eigenvalue | PercentVariance | CumulativeVariance |
|---|---|---|---|
| Factor 1 | 3.73 | 46.67 | 46.67 |
| Factor 2 | 1.79 | 22.38 | 69.05 |
| Factor 3 | 1.12 | 14.01 | 83.06 |
| Factor 4 | 0.61 | 7.64 | 90.70 |
| Factor 5 | 0.25 | 3.11 | 93.81 |
| Factor 6 | 0.21 | 2.67 | 96.48 |
| Factor 7 | 0.15 | 1.90 | 98.38 |
| Factor 8 | 0.13 | 1.62 | 100.00 |
Principal Components (Table 15.5)
# install.packages(c("psych","knitr","kableExtra","dplyr","tibble")) # if needed
library(psych)
library(dplyr)
library(tibble)
library(knitr)
library(kableExtra)
apa_pca_table <- function(data, k = NULL, digits = 2, cutoff = 0) {
# 1) Correlation matrix (pairwise) and n for reference
R <- cor(data, use = "pairwise.complete.obs")
n <- nrow(data)
# 2) Eigenvalues and default number of components (Kaiser, >= 1.0) if k not given
ev <- eigen(R, symmetric = TRUE, only.values = TRUE)$values
if (is.null(k)) k <- max(1, sum(ev >= 1)) # at least 1 component
# 3) Unrotated PCA via psych::principal on the correlation matrix
pc <- principal(r = R, nfactors = k, rotate = "none", n.obs = n)
# 4) Loadings table (variables x components), APA formatting
L <- as.matrix(pc$loadings)
colnames(L) <- paste0("PC", seq_len(ncol(L)))
rownames(L) <- colnames(R)
# optional: blank very small loadings (cutoff default 0 = show all)
if (cutoff > 0) L[abs(L) < cutoff] <- NA
L_fmt <- apply(L, 2, function(col) formatC(col, digits = digits, format = "f"))
L_df <- as.data.frame(L_fmt, stringsAsFactors = FALSE) %>%
rownames_to_column("Variable")
# 5) Eigenvalues table for the same k components
eig_df <- data.frame(
Component = paste0("PC", seq_len(k)),
Eigenvalue = ev[seq_len(k)],
PercentVariance = ev[seq_len(k)] / sum(ev) * 100,
CumulativeVariance = cumsum(ev[seq_len(k)]) / sum(ev) * 100
) %>%
mutate(across(-Component, ~ round(.x, digits)))
list(loadings = L_df, eigen = eig_df)
}
# ---- Run on your data frame RFD ----
out <- apa_pca_table(RFD, k = 8, digits = 2, cutoff = 0) # Extracted 8 components
# APA-style Loadings table (Unrotated)
kable(out$loadings,
caption = "Table 15.5 Unrotated Principal Component Loadings (Correlation Matrix, APA Format)",
align = c("l", rep("c", ncol(out$loadings) - 1))) %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
column_spec(1, bold = TRUE)| Variable | PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 |
|---|---|---|---|---|---|---|---|---|
| rridcw0 | 0.79 | -0.34 | 0.25 | -0.36 | -0.00 | -0.01 | -0.04 | 0.25 |
| rriddw0 | 0.79 | -0.33 | 0.25 | -0.36 | -0.01 | 0.01 | 0.03 | -0.25 |
| rridgw0 | 0.65 | 0.25 | -0.63 | -0.08 | -0.12 | 0.30 | 0.02 | 0.01 |
| rridjw0 | 0.68 | 0.26 | -0.60 | -0.07 | 0.10 | -0.32 | -0.01 | -0.01 |
| rridkw0 | 0.78 | -0.34 | 0.05 | 0.39 | 0.34 | 0.11 | -0.00 | 0.00 |
| rridqw0 | 0.77 | -0.31 | 0.10 | 0.43 | -0.33 | -0.10 | 0.01 | 0.00 |
| rriduw0 | 0.43 | 0.78 | 0.35 | 0.03 | 0.03 | -0.01 | 0.27 | 0.03 |
| rridvw0 | 0.45 | 0.78 | 0.32 | 0.05 | -0.01 | 0.02 | -0.27 | -0.03 |
First Principal Component (Figure 15.6)
Path Diagram
(Numbers do not exactly match those in the book because data are simulated)
First Principal Component
R Code
library(lavaan)
RFD_std <- as.data.frame(scale(RFD))
# 2. Build lavaan syntax: one factor loading on all variables
FirstPCmodel <- paste0("PC1 =~ ", paste(vars, collapse = " + "))
# 3. Fit the model (std.lv = TRUE makes the latent variable variance = 1)
FirstPCfit <- cfa(FirstPCmodel, data = RFD_std, std.lv = TRUE, estimator = "ULS")
# 4. Summary of loadings
summary(FirstPCfit, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 22 iterations
##
## Estimator ULS
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 2422
##
## Model Test User Model:
##
## Test statistic 2611.338
## Degrees of freedom 20
## P-value (Unknown) NA
##
## Model Test Baseline Model:
##
## Test statistic 13221.019
## Degrees of freedom 28
## P-value NA
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.804
## Tucker-Lewis Index (TLI) 0.725
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.231
## 90 Percent confidence interval - lower 0.224
## 90 Percent confidence interval - upper 0.239
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.173
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## PC1 =~
## rridcw0 0.777 0.013 57.843 0.000 0.777 0.777
## rriddw0 0.780 0.013 57.983 0.000 0.780 0.780
## rridgw0 0.563 0.012 45.678 0.000 0.563 0.563
## rridjw0 0.592 0.012 47.615 0.000 0.592 0.592
## rridkw0 0.749 0.013 56.500 0.000 0.749 0.749
## rridqw0 0.741 0.013 56.105 0.000 0.741 0.741
## rriduw0 0.331 0.012 28.452 0.000 0.331 0.331
## rridvw0 0.349 0.012 29.846 0.000 0.349 0.349
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridcw0 0.397 0.029 13.624 0.000 0.397 0.397
## .rriddw0 0.392 0.029 13.424 0.000 0.392 0.392
## .rridgw0 0.683 0.025 27.768 0.000 0.683 0.683
## .rridjw0 0.649 0.025 25.848 0.000 0.649 0.649
## .rridkw0 0.440 0.028 15.481 0.000 0.440 0.440
## .rridqw0 0.452 0.028 16.010 0.000 0.452 0.452
## .rriduw0 0.890 0.022 40.946 0.000 0.890 0.890
## .rridvw0 0.878 0.022 40.117 0.000 0.878 0.878
## PC1 1.000 1.000 1.000require(lavaangui)## Loading required package: lavaangui## Warning: package 'lavaangui' was built under R version 4.5.1## This is lavaangui 0.2.5
## lavaangui is BETA software! Please report any bugs at https://github.com/karchjd/lavaangui/issues#plot_lavaan(FirstPCfit)Implied and Residual Correlation Matrices for First Component
# Factor loadings
lambda <- c(
rridcw0 = 0.777,
rriddw0 = 0.780,
rridgw0 = 0.563,
rridjw0 = 0.592,
rridkw0 = 0.749,
rridqw0 = 0.741,
rriduw0 = 0.331,
rridvw0 = 0.349
)
# Convert to matrix (8×1)
Lambda <- matrix(lambda, ncol = 1)
# Assume factor variance = 1
Phi <- matrix(1, nrow = 1, ncol = 1)
# Unique variances: here we assume standardized observed variables
# so unique variance = 1 - loading^2
Theta <- diag(1 - lambda^2)
# Model-implied covariance
Sigma_hat <- Lambda %*% Phi %*% t(Lambda)
# Add row/col names
dimnames(Sigma_hat) <- list(names(lambda), names(lambda))
round(Sigma_hat,2)## rridcw0 rriddw0 rridgw0 rridjw0 rridkw0 rridqw0 rriduw0 rridvw0
## rridcw0 0.60 0.61 0.44 0.46 0.58 0.58 0.26 0.27
## rriddw0 0.61 0.61 0.44 0.46 0.58 0.58 0.26 0.27
## rridgw0 0.44 0.44 0.32 0.33 0.42 0.42 0.19 0.20
## rridjw0 0.46 0.46 0.33 0.35 0.44 0.44 0.20 0.21
## rridkw0 0.58 0.58 0.42 0.44 0.56 0.56 0.25 0.26
## rridqw0 0.58 0.58 0.42 0.44 0.56 0.55 0.25 0.26
## rriduw0 0.26 0.26 0.19 0.20 0.25 0.25 0.11 0.12
## rridvw0 0.27 0.27 0.20 0.21 0.26 0.26 0.12 0.12# Observed Correlations
RFD_cor <- cor(RFD, use = "pairwise.complete.obs")
#Residual matrix
Resid <- RFD_cor- Sigma_hat
round(Resid,2)## rridcw0 rriddw0 rridgw0 rridjw0 rridkw0 rridqw0 rriduw0 rridvw0
## rridcw0 0.40 0.26 -0.14 -0.13 0.02 0.01 -0.11 -0.11
## rriddw0 0.26 0.39 -0.13 -0.13 0.01 0.01 -0.10 -0.11
## rridgw0 -0.14 -0.13 0.68 0.45 -0.07 -0.08 0.07 0.09
## rridjw0 -0.13 -0.13 0.45 0.65 -0.06 -0.09 0.09 0.10
## rridkw0 0.02 0.01 -0.07 -0.06 0.44 0.20 -0.14 -0.14
## rridqw0 0.01 0.01 -0.08 -0.09 0.20 0.45 -0.11 -0.10
## rriduw0 -0.11 -0.10 0.07 0.09 -0.14 -0.11 0.89 0.73
## rridvw0 -0.11 -0.11 0.09 0.10 -0.14 -0.10 0.73 0.88Extract First Two Components (Figure 15.7)
Path Diagram
First Two Principal Components
R Program
library(lavaan)
library(dplyr)
# --- Prep: standardize so we're on the correlation scale (PCA scale) ---
RFD_std <- as.data.frame(scale(RFD))
vars <- colnames(RFD_std) # or your own character vector of variable names
# ---------- 1) First component via 1-factor CFA ----------
mod_pc1 <- paste0("PC1 =~ ", paste(vars, collapse = " + "))
fit_pc1 <- cfa(
mod_pc1,
data = RFD_std,
std.lv = TRUE, # Var(PC1)=1
estimator = "ULS" # or "ML"
)
summary(fit_pc1, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 22 iterations
##
## Estimator ULS
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 2422
##
## Model Test User Model:
##
## Test statistic 2611.338
## Degrees of freedom 20
## P-value (Unknown) NA
##
## Model Test Baseline Model:
##
## Test statistic 13221.019
## Degrees of freedom 28
## P-value NA
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.804
## Tucker-Lewis Index (TLI) 0.725
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.231
## 90 Percent confidence interval - lower 0.224
## 90 Percent confidence interval - upper 0.239
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.173
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## PC1 =~
## rridcw0 0.777 0.013 57.843 0.000 0.777 0.777
## rriddw0 0.780 0.013 57.983 0.000 0.780 0.780
## rridgw0 0.563 0.012 45.678 0.000 0.563 0.563
## rridjw0 0.592 0.012 47.615 0.000 0.592 0.592
## rridkw0 0.749 0.013 56.500 0.000 0.749 0.749
## rridqw0 0.741 0.013 56.105 0.000 0.741 0.741
## rriduw0 0.331 0.012 28.452 0.000 0.331 0.331
## rridvw0 0.349 0.012 29.846 0.000 0.349 0.349
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridcw0 0.397 0.029 13.624 0.000 0.397 0.397
## .rriddw0 0.392 0.029 13.424 0.000 0.392 0.392
## .rridgw0 0.683 0.025 27.768 0.000 0.683 0.683
## .rridjw0 0.649 0.025 25.848 0.000 0.649 0.649
## .rridkw0 0.440 0.028 15.481 0.000 0.440 0.440
## .rridqw0 0.452 0.028 16.010 0.000 0.452 0.452
## .rriduw0 0.890 0.022 40.946 0.000 0.890 0.890
## .rridvw0 0.878 0.022 40.117 0.000 0.878 0.878
## PC1 1.000 1.000 1.000# extract standardized loadings (lambda) from the first fit
lam1 <- standardizedSolution(fit_pc1) %>%
filter(op == "=~", lhs == "PC1", rhs %in% vars) %>%
select(rhs, est.std) %>%
tibble::deframe() # named numeric vector: names=vars, values=loadings
# ---------- 2) Fix PC1 loadings; add a second orthogonal component (PC2) ----------
# Build the fixed-loading part for PC1
pc1_fixed <- paste0(
"PC1 =~ ",
paste(sprintf("%.10f*%s", lam1[vars], vars), collapse = " + ")
)
# PC2 loads freely on all variables; PC1 ⟂ PC2; both factors var=1 (std.lv=TRUE)
pc2_free <- paste0("PC2 =~ ", paste(vars, collapse = " + "))
mod_pc2 <- paste(
pc1_fixed,
pc2_free,
"PC1 ~~ 0*PC2", # orthogonal, like PCA components
sep = "\n"
)
fit_pc2 <- cfa(
mod_pc2,
data = RFD_std,
std.lv = TRUE, # Var(PC1)=Var(PC2)=1
estimator = "ULS"
)
summary(fit_pc2, standardized = TRUE, fit.measures = TRUE)## lavaan 0.6-19 ended normally after 27 iterations
##
## Estimator ULS
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 2422
##
## Model Test User Model:
##
## Test statistic 801.574
## Degrees of freedom 20
## P-value (Unknown) NA
##
## Model Test Baseline Model:
##
## Test statistic 13221.019
## Degrees of freedom 28
## P-value NA
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.941
## Tucker-Lewis Index (TLI) 0.917
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.127
## 90 Percent confidence interval - lower 0.120
## 90 Percent confidence interval - upper 0.135
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.096
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## PC1 =~
## rridcw0 0.777 0.777 0.777
## rriddw0 0.780 0.780 0.780
## rridgw0 0.563 0.563 0.563
## rridjw0 0.592 0.592 0.592
## rridkw0 0.749 0.749 0.749
## rridqw0 0.741 0.741 0.741
## rriduw0 0.331 0.331 0.331
## rridvw0 0.349 0.349 0.349
## PC2 =~
## rridcw0 0.191 0.017 11.385 0.000 0.191 0.191
## rriddw0 0.186 0.017 11.129 0.000 0.186 0.186
## rridgw0 -0.203 0.017 -12.064 0.000 -0.203 -0.203
## rridjw0 -0.212 0.017 -12.611 0.000 -0.212 -0.212
## rridkw0 0.196 0.017 11.669 0.000 0.196 0.196
## rridqw0 0.164 0.017 9.863 0.000 0.164 0.164
## rriduw0 -0.797 0.032 -24.668 0.000 -0.797 -0.797
## rridvw0 -0.826 0.033 -24.721 0.000 -0.826 -0.826
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## PC1 ~~
## PC2 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rridcw0 0.360 0.021 16.918 0.000 0.360 0.360
## .rriddw0 0.357 0.021 16.812 0.000 0.357 0.357
## .rridgw0 0.642 0.021 29.962 0.000 0.642 0.642
## .rridjw0 0.604 0.022 28.028 0.000 0.604 0.604
## .rridkw0 0.401 0.021 18.794 0.000 0.401 0.401
## .rridqw0 0.425 0.021 20.172 0.000 0.425 0.425
## .rriduw0 0.254 0.055 4.591 0.000 0.254 0.254
## .rridvw0 0.196 0.059 3.333 0.001 0.196 0.196
## PC1 1.000 1.000 1.000
## PC2 1.000 1.000 1.000# Optional: extract tidy loadings for PC2
loadings_pc2 <- standardizedSolution(fit_pc2) %>%
filter(op == "=~", lhs == "PC2", rhs %in% vars) %>%
select(variable = rhs, loading = est.std) %>%
arrange(desc(abs(loading)))
loadings_pc2## variable loading
## 1 rridvw0 -0.826
## 2 rriduw0 -0.797
## 3 rridjw0 -0.212
## 4 rridgw0 -0.203
## 5 rridkw0 0.196
## 6 rridcw0 0.191
## 7 rriddw0 0.186
## 8 rridqw0 0.164#plot_lavaan(fit_pc2)First Four Principal Components
Path Diagram
First Four Principal Components
R Program
library(lavaan)
library(dplyr)
library(semPlot) # for semPaths
# --- Prep ---
RFD_std <- as.data.frame(scale(RFD))
vars <- colnames(RFD_std)
# --- Function to get standardized loadings from a lavaan fit ---
get_loadings <- function(fit, factor_name, vars) {
standardizedSolution(fit) %>%
filter(op == "=~", lhs == factor_name, rhs %in% vars) %>%
select(rhs, est.std) %>%
tibble::deframe()
}
# 1) PC1
fit_pc1 <- cfa(paste0("PC1 =~ ", paste(vars, collapse = " + ")),
data = RFD_std, std.lv = TRUE, estimator = "ULS")
lam1 <- get_loadings(fit_pc1, "PC1", vars)
# 2) PC2 with PC1 fixed
pc1_fixed <- paste0("PC1 =~ ", paste(sprintf("%.10f*%s", lam1, vars), collapse = " + "))
pc2_free <- paste0("PC2 =~ ", paste(vars, collapse = " + "))
fit_pc2 <- cfa(paste(pc1_fixed, pc2_free, "PC1 ~~ 0*PC2", sep = "\n"),
data = RFD_std, std.lv = TRUE, estimator = "ULS")
lam2 <- get_loadings(fit_pc2, "PC2", vars)
# 3) PC3 with PC1 & PC2 fixed
pc2_fixed <- paste0("PC2 =~ ", paste(sprintf("%.10f*%s", lam2, vars), collapse = " + "))
pc3_free <- paste0("PC3 =~ ", paste(vars, collapse = " + "))
fit_pc3 <- cfa(paste(pc1_fixed, pc2_fixed, pc3_free,
"PC1 ~~ 0*PC2", "PC1 ~~ 0*PC3", "PC2 ~~ 0*PC3", sep = "\n"),
data = RFD_std, std.lv = TRUE, estimator = "ULS")
lam3 <- get_loadings(fit_pc3, "PC3", vars)
# 4) PC4 with PC1–PC3 fixed
pc3_fixed <- paste0("PC3 =~ ", paste(sprintf("%.10f*%s", lam3, vars), collapse = " + "))
pc4_free <- paste0("PC4 =~ ", paste(vars, collapse = " + "))
final_model <- paste(pc1_fixed, pc2_fixed, pc3_fixed, pc4_free,
"PC1 ~~ 0*PC2", "PC1 ~~ 0*PC3", "PC1 ~~ 0*PC4",
"PC2 ~~ 0*PC3", "PC2 ~~ 0*PC4",
"PC3 ~~ 0*PC4",
sep = "\n")
fit_pc4 <- cfa(final_model, data = RFD_std, std.lv = TRUE, estimator = "ULS")
# --- Plot the 4-component model ---
#plot_lavaan(fit_pc4)Table of Reduced Correlation Matrix (i.e., SMC’s on diagonal), Eigenvectors and Unrotated PA loadings
# install.packages(c("dplyr","tibble","knitr","kableExtra")) # if needed
library(dplyr)
library(tibble)
library(knitr)
library(kableExtra)
# If you don't already have RFD_std:
# RFD_std <- as.data.frame(scale(RFD))
# 1) Correlation matrix and SMCs on the diagonal
R <- cor(RFD_std, use = "pairwise.complete.obs")
# Squared Multiple Correlations (SMC) for each variable
# SMC = 1 - 1 / diag(R^{-1}) (for standardized variables)
invR <- solve(R)
smc <- 1 - 1 / diag(invR)
R_smc <- R
diag(R_smc) <- smc
# 2) Eigenvectors of the (standard) correlation matrix R
eig <- eigen(R, symmetric = TRUE)
eigenvectors <- eig$vectors
colnames(eigenvectors) <- paste0("EV", seq_len(ncol(eigenvectors)))
rownames(eigenvectors) <- colnames(R)
# 3) First four unrotated, un-iterated principal-axis loadings
# Use the initial PA step: replace diag(R) with SMCs, then eigen-decompose.
eig_pa <- eigen(R_smc, symmetric = TRUE)
k <- 4
# guard against negative/near-zero eigenvalues in the first k
lam <- pmax(eig_pa$values[1:k], 0)
V <- eig_pa$vectors[, 1:k, drop = FALSE]
# Unrotated loadings (initial PA): V * sqrt(lambda)
Loadings_PA <- V %*% diag(sqrt(lam), nrow = k, ncol = k)
rownames(Loadings_PA) <- colnames(R)
colnames(Loadings_PA) <- paste0("PA", 1:k)
# ---- APA-style tables (two decimals) ----
# Helper to format to exactly two decimals
fmt2 <- function(x) formatC(x, digits = 2, format = "f")
# (A) Correlation matrix with SMCs on the diagonal
R_smc_fmt <- apply(R_smc, 2, fmt2) %>% as.data.frame() %>%
rownames_to_column("Variable")
kable(R_smc_fmt,
caption = "Table 15.8(1)\nCorrelation Matrix of RFD_std with SMCs on the Diagonal",
align = c("l", rep("c", ncol(R_smc_fmt) - 1))) %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
column_spec(1, bold = TRUE)| Variable | rridcw0 | rriddw0 | rridgw0 | rridjw0 | rridkw0 | rridqw0 | rriduw0 | rridvw0 |
|---|---|---|---|---|---|---|---|---|
| rridcw0 | 0.77 | 0.87 | 0.30 | 0.33 | 0.60 | 0.59 | 0.15 | 0.16 |
| rriddw0 | 0.87 | 0.77 | 0.31 | 0.33 | 0.60 | 0.59 | 0.16 | 0.16 |
| rridgw0 | 0.30 | 0.31 | 0.62 | 0.78 | 0.35 | 0.34 | 0.26 | 0.29 |
| rridjw0 | 0.33 | 0.33 | 0.78 | 0.63 | 0.38 | 0.35 | 0.28 | 0.31 |
| rridkw0 | 0.60 | 0.60 | 0.35 | 0.38 | 0.62 | 0.75 | 0.11 | 0.12 |
| rridqw0 | 0.59 | 0.59 | 0.34 | 0.35 | 0.75 | 0.60 | 0.14 | 0.16 |
| rriduw0 | 0.15 | 0.16 | 0.26 | 0.28 | 0.11 | 0.14 | 0.72 | 0.85 |
| rridvw0 | 0.16 | 0.16 | 0.29 | 0.31 | 0.12 | 0.16 | 0.85 | 0.72 |
# (B) Eigenvectors of R
EV_fmt <- apply(eigenvectors, 2, fmt2) %>% as.data.frame() %>%
rownames_to_column("Variable")
kable(EV_fmt,
caption = "Table 15.8(2)\nEigenvectors of the Correlation Matrix (RFD_std)",
align = c("l", rep("c", ncol(EV_fmt) - 1))) %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
column_spec(1, bold = TRUE) %>%
add_header_above(c(" " = 1, "Eigenvectors" = ncol(EV_fmt) - 1))| Variable | EV1 | EV2 | EV3 | EV4 | EV5 | EV6 | EV7 | EV8 |
|---|---|---|---|---|---|---|---|---|
| rridcw0 | -0.41 | 0.25 | -0.24 | -0.46 | 0.01 | -0.02 | -0.10 | 0.70 |
| rriddw0 | -0.41 | 0.25 | -0.24 | -0.46 | 0.02 | 0.02 | 0.09 | -0.70 |
| rridgw0 | -0.34 | -0.19 | 0.59 | -0.10 | 0.24 | 0.65 | 0.05 | 0.03 |
| rridjw0 | -0.35 | -0.19 | 0.56 | -0.09 | -0.21 | -0.69 | -0.03 | -0.02 |
| rridkw0 | -0.40 | 0.25 | -0.05 | 0.50 | -0.68 | 0.24 | -0.01 | 0.00 |
| rridqw0 | -0.40 | 0.23 | -0.10 | 0.55 | 0.66 | -0.21 | 0.03 | 0.00 |
| rriduw0 | -0.22 | -0.59 | -0.33 | 0.04 | -0.05 | -0.02 | 0.69 | 0.10 |
| rridvw0 | -0.23 | -0.59 | -0.31 | 0.06 | 0.02 | 0.05 | -0.70 | -0.09 |
# (C) First four unrotated, un-iterated principal-axis loadings
L_fmt <- apply(Loadings_PA, 2, fmt2) %>% as.data.frame() %>%
rownames_to_column("Variable")
kable(L_fmt,
caption = "Table 15.8(3)\nUnrotated, Un-iterated Principal-Axis Loadings (First Four Factors)",
align = c("l", rep("c", ncol(L_fmt) - 1))) %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
column_spec(1, bold = TRUE) %>%
add_header_above(c(" " = 1, "Initial PA Loadings" = ncol(L_fmt) - 1))| Variable | PA1 | PA2 | PA3 | PA4 |
|---|---|---|---|---|
| rridcw0 | -0.78 | 0.32 | -0.23 | -0.23 |
| rriddw0 | -0.78 | 0.32 | -0.23 | -0.23 |
| rridgw0 | -0.61 | -0.21 | 0.52 | -0.07 |
| rridjw0 | -0.64 | -0.21 | 0.51 | -0.07 |
| rridkw0 | -0.73 | 0.28 | -0.02 | 0.28 |
| rridqw0 | -0.72 | 0.26 | -0.06 | 0.30 |
| rriduw0 | -0.42 | -0.73 | -0.28 | 0.02 |
| rridvw0 | -0.44 | -0.73 | -0.25 | 0.03 |
Eigenvalues by Dimensions Plot
R Program
library(paran)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectparan(AllRFDData, iterations = 5000, centile = 0, quietly = FALSE,
status = TRUE, all = TRUE, cfa = TRUE, graph = TRUE, color = TRUE,
col = c("black", "red", "blue"), lty = c(1, 2, 3), lwd = 1, legend = TRUE,seed = 0)##
## Using eigendecomposition of correlation matrix.
## Computing: 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
##
##
## Results of Horn's Parallel Analysis for factor retention
## 5000 iterations, using the mean estimate
##
## --------------------------------------------------
## Factor Adjusted Unadjusted Estimated
## Eigenvalue Eigenvalue Bias
## --------------------------------------------------
## No components passed.
## --------------------------------------------------
## 1 8.035200 8.222706 0.187505
## 2 3.125955 3.285825 0.159870
## 3 1.712588 1.851410 0.138821
## 4 0.751993 0.872420 0.120427
## 5 0.481457 0.585173 0.103716
## 6 0.087695 0.175715 0.088019
## 7 0.075378 0.148489 0.073111
## 8 0.017672 0.076538 0.058865
## 9 -0.016800 0.028235 0.045036
## 10 -0.021868 0.009770 0.031638
## 11 -0.025985 -0.00742 0.018556
## 12 -0.031575 -0.02582 0.005748
## 13 -0.021373 -0.02852 -0.00715
## 14 -0.038524 -0.05859 -0.02006
## 15 -0.030565 -0.06348 -0.03291
## 16 -0.030626 -0.07613 -0.04550
## 17 -0.035716 -0.09434 -0.05862
## 18 -0.029785 -0.10169 -0.07190
## 19 -0.020696 -0.10618 -0.08548
## 20 -0.021428 -0.12109 -0.09966
## 21 -0.013841 -0.12883 -0.11499
## 22 -0.009619 -0.14178 -0.13216
## 23 -0.050226 -0.20424 -0.15401
## --------------------------------------------------
##
## Adjusted eigenvalues > 0 indicate dimensions to retain.
## (8 factors retained)
#Another approach: fa.parallel
library(nFactors)
parallel = fa.parallel(AllRFDData,
fm = 'ml',
fa = 'fa',
n.iter = 50,
SMC = TRUE,
quant = .95)
## Parallel analysis suggests that the number of factors = 8 and the number of components = NAVelicer’s MAP test (Table 15.1)
# Load the psych package
library(psych)
# Run Velicer's MAP test on raw data
map_result <- VSS(AllRFDData, rotate = "none", fm = "minres", n = ncol(AllRFDData))
# The MAP test results
map_result$map # Minimum Average Partial values for each factor## [1] 0.06672034 0.04420930 0.02111018 0.02118145 0.01792197 0.02355937
## [7] 0.02903513 0.03609009 0.04581622 0.06009931 0.08144871 0.09087631
## [13] 0.08399545 0.09473116 0.10452696 0.12588086 0.16923530 0.24193623
## [19] 0.26628233 0.32988296 0.48494337 1.00000000 NAmap_result$n.obs # Number of observations## NULL# The number of factors suggested by the MAP test is the one with the lowest MAP value
# (See the R Console)
suggested_factors <- which.min(map_result$map)
cat("Velicer's MAP test suggests", suggested_factors, "factors.\n")## Velicer's MAP test suggests 5 factors.Cattell’s Scree test
R Program
library(nFactors)
ev <- eigen(cor(AllRFDData, use = "pairwise.complete.obs"), symmetric = TRUE)$values
# "Optimal Coordinates" are the Scree test, along with the parallel results
ns <- nScree(eig = ev) # computes CNG, Kaiser, PA, etc.
plotnScree(ns) # has a dedicated plotting function