Exploratory and Confirmatory Factor Analysis on SF-36 Scales

Last updated: July 15 2023

Description of Example

The SF-36 questionnaire measures self-reported functional status. Its 36 items are combined to form 8 subscales yielding two summary measures: physical and mental health.

The physical health measure includes four subscales: physical functioning (10 items), role-physical (4 items), bodily pain (2 items), and general health (5 items).
The mental health measure includes the other four subscales: vitality (4 items), social functioning (2 items), role-emotional (3 items), and mental health (5 items).
This example uses data from 249 women (age 31 to 87, median 64) who completed the SF-36 following heart surgery. The subscale variable names are: MENHLTH, SOCFUNC, ROLEEMO, VITAL, PHSFUNC, ROLPHYS, BODPAIN, and GENHLTH.
First, an Exploratory Factor Analysis with 2 uncorrelated factors is run using the psych package in which each subscale loads on both factors.
Second, a Confirmatory Factor Analysis with 2 correlated factors is run using the lavaan package in which each subscale loads on only one hypothesized factor.

Setup

Required Packages

knitr::opts_chunk$set(error = TRUE, message = FALSE, warning = TRUE)

library(rio)
library(readr)
library(here)
library(dplyr)
library(lavaan)
library(psych)
library(semoutput)
library(semPlot)
library(sjPlot)
library(GPArotation)
library(semTools)

Import Data

## Import Data
hrtsrg <- rio::import("C:/Users/jbeckste/Desktop/WOMEN_HRTSRG_noLabels.SAV")
attach(hrtsrg)
data <- data.frame(MENHLTH, SOCFUNC, ROLEEMO, VITAL,   PHSFUNC, ROLPHYS, BODPAIN, GENHLTH)

---

Descriptives

# Prints basic descriptive statistics
sem_descriptives(data)

Descriptive Statistics

Variable	n	Mean	SD	min	max	Skewness	Kurtosis	% Missing
MENHLTH	249	18.09	4.43	6	25	-0.68	-0.18	0
SOCFUNC	249	7.15	2.17	2	10	-0.34	-0.74	0
ROLEEMO	249	10.59	3.44	3	15	-0.31	-0.87	0
VITAL	249	11.30	3.34	4	20	-0.16	-0.48	0
PHSFUNC	249	19.69	4.93	10	30	0.17	-0.81	0
ROLPHYS	249	10.83	3.88	4	20	0.23	-0.44	0
BODPAIN	249	7.21	2.38	2	12	0.20	-0.45	0
GENHLTH	249	16.44	4.03	7	25	-0.13	-0.68	0
Total N = 249

---

Correlation Matrix

Bivariate Summary Display

Histograms are shown on the diagonal, scatterplots showing bivariate regression line and centroid are below the diagnonal, correlations are above the diagonal.

psych::pairs.panels(data, lm=T, cex.cor=.7, jiggle=T, factor=3,
                    show.points = T, ellipses=T, pch=".")

Heatmap

Larger correlations are more darkly shaded. Ideally, subscales within each measure should show stronger correlations than subscales from different measures.

#makes heatmap of corr mtx; SORT df first based on FA
psych::cor.plot(data, pval=F,diag=T, scale=F, show.legend=F, upper=F, n=51,
                main="", cex=.8, xlas=2)

Nice Table

# Uses sjPlot to print a nice looking correlation table
tab_corr(data, na.deletion = "listwise", digits = 2, triangle = "lower",
         title="Correlations among SF-36 Scales",
         string.diag=c('1','1','1','1','1','1','1','1'))

Correlations among SF-36 Scales

	MENHLTH	SOCFUNC	ROLEEMO	VITAL	PHSFUNC	ROLPHYS	BODPAIN	GENHLTH
MENHLTH	1
SOCFUNC	0.60***	1
ROLEEMO	0.63***	0.55***	1
VITAL	0.58***	0.58***	0.42***	1
PHSFUNC	0.20**	0.40***	0.31***	0.43***	1
ROLPHYS	0.35***	0.51***	0.53***	0.48***	0.55***	1
BODPAIN	0.34***	0.45***	0.33***	0.36***	0.46***	0.45***	1
GENHLTH	0.51***	0.48***	0.39***	0.54***	0.36***	0.34***	0.34***	1
Computed correlation used pearson-method with listwise-deletion.

---

EFA Results

Based on the scree plot below, a 2-factor solution seems reasonable. This plot was generated by psych::VSS.scree

efa_data <- data
## Determine the number of factors to extract by inspection of scree plot
psych::VSS.scree(efa_data)

## Conduct EFA analysis with 2 factors determined from scree plot
efa_fit <- psych::fa(efa_data, fm = "ml", nfactors = 2, rotate = "varimax")

Summary Output

efa_method(efa_fit)

Extraction Method

Sample.Size	Method	Factors.Extracted	Rotation
249	ml	2	Varimax

efa_var(efa_fit)

Total Variance Explained

Factor	Eigenvalue	Proportion Var	Cumulative Var
ML1	2.391	0.299	0.299
ML2	2.030	0.254	0.553

efa_loadings(efa_fit)

Factor Loadings

	Factors
Variable	ML1	ML2	h2
MENHLTH	0.956	0.094	0.923
SOCFUNC	0.580	0.487	0.574
ROLEEMO	0.622	0.351	0.510
VITAL	0.558	0.461	0.523
PHSFUNC	0.133	0.735	0.558
ROLPHYS	0.297	0.702	0.581
BODPAIN	0.293	0.532	0.369
GENHLTH	0.501	0.365	0.384

efa_rotmatrix(efa_fit)

Factor Rotation Matrix

Factor	1	2
1	0.946	0.324
2	-0.324	0.946

Diagram Output

psych::fa.diagram(efa_fit, cut=.01, simple=F, adj=2, digits=2,
                  main="Exploratory Factor Analysis using ML")

Full Output

efa_fit

## Factor Analysis using method =  ml
## Call: psych::fa(r = efa_data, nfactors = 2, rotate = "varimax", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##          ML1  ML2   h2    u2 com
## MENHLTH 0.96 0.09 0.92 0.077 1.0
## SOCFUNC 0.58 0.49 0.57 0.426 1.9
## ROLEEMO 0.62 0.35 0.51 0.490 1.6
## VITAL   0.56 0.46 0.52 0.477 1.9
## PHSFUNC 0.13 0.73 0.56 0.442 1.1
## ROLPHYS 0.30 0.70 0.58 0.419 1.3
## BODPAIN 0.29 0.53 0.37 0.631 1.6
## GENHLTH 0.50 0.36 0.38 0.616 1.8
## 
##                        ML1  ML2
## SS loadings           2.39 2.03
## Proportion Var        0.30 0.25
## Cumulative Var        0.30 0.55
## Proportion Explained  0.54 0.46
## Cumulative Proportion 0.54 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  28  with the objective function =  3.43 with Chi Square =  838.42
## df of  the model are 13  and the objective function was  0.17 
## 
## The root mean square of the residuals (RMSR) is  0.04 
## The df corrected root mean square of the residuals is  0.06 
## 
## The harmonic n.obs is  249 with the empirical chi square  20.96  with prob <  0.074 
## The total n.obs was  249  with Likelihood Chi Square =  41.84  with prob <  7e-05 
## 
## Tucker Lewis Index of factoring reliability =  0.923
## RMSEA index =  0.094  and the 90 % confidence intervals are  0.063 0.127
## BIC =  -29.89
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    ML1  ML2
## Correlation of (regression) scores with factors   0.96 0.88
## Multiple R square of scores with factors          0.92 0.77
## Minimum correlation of possible factor scores     0.84 0.54

---

CFA in LISREL Notation

The goal of CFA is to reproduce the covariance matrix, \(\mathbf\Sigma\), among the observed variables. This is done using a hypothesized model that consists of three types of parameters contained in three matrices:

\[ \mathbf{\Sigma^* = \Lambda_x \Phi \Lambda'_x + \Theta_\delta} \] where \(\mathbf\Sigma^*\) is the estimated covariance matrix among the eight subscales, \(\mathbf\Lambda_x\) is a matrix of factor loadings, \(\mathbf\Phi\) is the covariance matrix of the latent variables (factors), and \(\mathbf\Theta_\delta\) is a matrix of measurement error variances.

Here is the \(\mathbf\Lambda_x\) for our model of the SF-36 subscales:

\[\mathbf\Lambda_x = \begin{bmatrix} \lambda_{1,1} & 0\\ \lambda_{2,1} & 0\\ \lambda_{3,1} & 0\\ \lambda_{4,1} & 0\\ 0 & \lambda_{5,2}\\ 0 & \lambda_{6,2}\\ 0 & \lambda_{7,2}\\ 0 & \lambda_{8,2}\\ \end{bmatrix} \] The rows correspond to the eight subscales and the columns to the two latent variables (MH and GH). Note the first four subscales load on the first factor (MH) and have zero loadings on the second factor (PH). The reverse pattern is true for the last four subscales.

Here is the \(\mathbf\Phi\) matrix for our model:

\[\mathbf\Phi = \begin{bmatrix} \phi_{1,1} & \phi_{1,2}\\ \phi_{2,1} & \phi_{2,2}\\ \end{bmatrix} \] with the factor variances on the diagonal and their covariance in the off-diagonal. Note that this is a symmetric matrix.

Finally, here is the \(\mathbf\Theta_\delta\) matrix for our model:

\[\mathbf\Theta_\delta = \begin{bmatrix} \theta_{1,1} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \theta_{2,2} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \theta_{3,3} & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \theta_{4,4} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \theta_{5,5} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \theta_{6,6} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \theta_{7,7} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \theta_{8,8}\\ \end{bmatrix} \] Note that \(\mathbf\Theta_\delta\) is a diagonal matrix. This reflects our hypothesis that the measurement error terms for the subscales are uncorreleated (i.e., all off-diagonal elements are hypothesized to be zero).

CFA Results

Based on the results of the exploratory analysis, we decide to run a confirmatory model specifying that the first 4 subscales load on one factor and the last 4 subscales load on a second factor.

# specify the 2-factor CFA model
model <- '
# latent factors
MH =~  MENHLTH + SOCFUNC + ROLEEMO + VITAL
PH =~  PHSFUNC + ROLPHYS + BODPAIN + GENHLTH'
# fit the model
fit <- lavaan::cfa(model = model, data = data, mimic = "lavaan", 
           estimator = "ML", missing = "listwise", check.post=T,
           std.lv = TRUE, std.ov = FALSE, test = "standard", 
           se = "standard", bootstrap = 1000)

Summary Output

sem_tables(fit)

Model Significance

Sample.Size	Chi.Square	df	p.value
249	125.511	19	0

Model Fit Measures

CFI	RMSEA	RMSEA.Lower	RMSEA.Upper	SRMR	AIC	BIC
0.871	0.15	0.126	0.176	0.07	9891.606	9951.403

Factor Loadings

		Standardized
Latent Factor	Indicator	Loadings	sig	p	Lower.CI	Upper.CI	SE	z
MH	MENHLTH	0.756	***	0	0.690	0.822	0.034	22.416
MH	SOCFUNC	0.798	***	0	0.738	0.857	0.030	26.228
MH	ROLEEMO	0.711	***	0	0.637	0.784	0.037	18.985
MH	VITAL	0.730	***	0	0.659	0.800	0.036	20.338
PH	PHSFUNC	0.643	***	0	0.553	0.732	0.046	14.080
PH	ROLPHYS	0.728	***	0	0.650	0.806	0.040	18.279
PH	BODPAIN	0.618	***	0	0.525	0.711	0.047	13.046
PH	GENHLTH	0.607	***	0	0.513	0.701	0.048	12.629

Latent Factor Correlations

Factor 1	Factor 2	r	sig	p	Lower.CI	Upper.CI	SE
MH	PH	0.846	***	0	0.771	0.921	0.038

Diagram Output

factors <- fit@pta$vnames$lv[[1]]
size <- .65
semPaths(fit, latents = factors, whatLabels = "std", layout = "tree2", 
         rotation = 4, style = "lisrel", optimizeLatRes = TRUE, 
         structural = FALSE, layoutSplit = FALSE,
         intercepts = FALSE, residuals = T, 
         curve = 1, curvature = 3, nCharNodes = 8, 
         sizeLat = 11 * size, sizeMan = 11 * size, sizeMan2 = 4 * size, 
         edge.label.cex = 1.2 * size, 
         edge.color = "#000000", edge.label.position = .40)

Residual Correlation Matrix

The values below are the residual correlations (i.e., observed correlation - model estimated correlation). Note the largest residual is |-.21| for the correlation between MENHLTH & PHSFUNC. The average of these residuals is provided by the standarized root mean square residual (SRMR).

sem_residuals(fit)

	MENHLTH	SOCFUNC	ROLEEMO	VITAL	PHSFUNC	ROLPHYS	BODPAIN	GENHLTH
MENHLTH
SOCFUNC	-0.01
ROLEEMO	0.10	-0.02
VITAL	0.03	0.00	-0.10
PHSFUNC	-0.21	-0.03	-0.07	0.03
ROLPHYS	-0.12	0.02	0.09	0.03	0.09
BODPAIN	-0.06	0.03	-0.04	-0.02	0.06	0.0
GENHLTH	0.12	0.07	0.03	0.17	-0.03	-0.1	-0.04

Full Output

Summary

summary(fit, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6.15 ended normally after 22 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        17
## 
##   Number of observations                           249
## 
## Model Test User Model:
##                                                       
##   Test statistic                               125.511
##   Degrees of freedom                                19
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               853.850
##   Degrees of freedom                                28
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.871
##   Tucker-Lewis Index (TLI)                       0.810
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -4928.803
##   Loglikelihood unrestricted model (H1)      -4866.048
##                                                       
##   Akaike (AIC)                                9891.606
##   Bayesian (BIC)                              9951.403
##   Sample-size adjusted Bayesian (SABIC)       9897.512
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.150
##   90 Percent confidence interval - lower         0.126
##   90 Percent confidence interval - upper         0.176
##   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.070
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   MH =~                                                                 
##     MENHLTH           3.345    0.253   13.224    0.000    3.345    0.756
##     SOCFUNC           1.730    0.121   14.252    0.000    1.730    0.798
##     ROLEEMO           2.438    0.201   12.153    0.000    2.438    0.711
##     VITAL             2.436    0.193   12.596    0.000    2.436    0.730
##   PH =~                                                                 
##     PHSFUNC           3.160    0.307   10.308    0.000    3.160    0.643
##     ROLPHYS           2.818    0.234   12.046    0.000    2.818    0.728
##     BODPAIN           1.470    0.150    9.822    0.000    1.470    0.618
##     GENHLTH           2.442    0.254    9.617    0.000    2.442    0.607
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   MH ~~                                                                 
##     PH                0.846    0.038   22.119    0.000    0.846    0.846
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .MENHLTH           8.387    0.957    8.762    0.000    8.387    0.428
##    .SOCFUNC           1.713    0.214    8.013    0.000    1.713    0.364
##    .ROLEEMO           5.826    0.624    9.339    0.000    5.826    0.495
##    .VITAL             5.210    0.571    9.121    0.000    5.210    0.468
##    .PHSFUNC          14.202    1.515    9.375    0.000   14.202    0.587
##    .ROLPHYS           7.064    0.855    8.258    0.000    7.064    0.471
##    .BODPAIN           3.500    0.365    9.599    0.000    3.500    0.618
##    .GENHLTH          10.210    1.054    9.686    0.000   10.210    0.631
##     MH                1.000                               1.000    1.000
##     PH                1.000                               1.000    1.000

Standarized Parameter Estimates

standardizedSolution(fit)

##        lhs op     rhs est.std    se      z pvalue ci.lower ci.upper
## 1       MH =~ MENHLTH   0.756 0.034 22.416      0    0.690    0.822
## 2       MH =~ SOCFUNC   0.798 0.030 26.228      0    0.738    0.857
## 3       MH =~ ROLEEMO   0.711 0.037 18.985      0    0.637    0.784
## 4       MH =~   VITAL   0.730 0.036 20.338      0    0.659    0.800
## 5       PH =~ PHSFUNC   0.643 0.046 14.080      0    0.553    0.732
## 6       PH =~ ROLPHYS   0.728 0.040 18.279      0    0.650    0.806
## 7       PH =~ BODPAIN   0.618 0.047 13.046      0    0.525    0.711
## 8       PH =~ GENHLTH   0.607 0.048 12.629      0    0.513    0.701
## 9  MENHLTH ~~ MENHLTH   0.428 0.051  8.399      0    0.328    0.528
## 10 SOCFUNC ~~ SOCFUNC   0.364 0.049  7.504      0    0.269    0.459
## 11 ROLEEMO ~~ ROLEEMO   0.495 0.053  9.306      0    0.391    0.599
## 12   VITAL ~~   VITAL   0.468 0.052  8.930      0    0.365    0.570
## 13 PHSFUNC ~~ PHSFUNC   0.587 0.059 10.011      0    0.472    0.702
## 14 ROLPHYS ~~ ROLPHYS   0.471 0.058  8.128      0    0.357    0.584
## 15 BODPAIN ~~ BODPAIN   0.618 0.059 10.570      0    0.504    0.733
## 16 GENHLTH ~~ GENHLTH   0.631 0.058 10.815      0    0.517    0.746
## 17      MH ~~      MH   1.000 0.000     NA     NA    1.000    1.000
## 18      PH ~~      PH   1.000 0.000     NA     NA    1.000    1.000
## 19      MH ~~      PH   0.846 0.038 22.119      0    0.771    0.921

Residuals Analysis

Displays residual correlations and z-scores indicating the magnitude of the residuals. Typically, z values > |2.54| are interpreted as significant and indicate where the model is doing a poor job of reproducing the observed correlations. Note that the largest residual correlation and |z| is between MENHLTH & PHSFUNC.

lavResiduals(fit)

## $type
## [1] "cor.bentler"
## 
## $cov
##         MENHLT SOCFUN ROLEEM  VITAL PHSFUN ROLPHY BODPAI GENHLT
## MENHLTH  0.000                                                 
## SOCFUNC -0.007  0.000                                          
## ROLEEMO  0.096 -0.018  0.000                                   
## VITAL    0.026 -0.004 -0.097  0.000                            
## PHSFUNC -0.212 -0.034 -0.074  0.034  0.000                     
## ROLPHYS -0.119  0.020  0.090  0.026  0.088  0.000              
## BODPAIN -0.060  0.034 -0.038 -0.017  0.065  0.002  0.000       
## GENHLTH  0.122  0.074  0.027  0.169 -0.031 -0.103 -0.038  0.000
## 
## $cov.z
##         MENHLT SOCFUN ROLEEM  VITAL PHSFUN ROLPHY BODPAI GENHLT
## MENHLTH  0.000                                                 
## SOCFUNC -0.492  0.000                                          
## ROLEEMO  3.736 -1.037  0.000                                   
## VITAL    1.211 -0.225 -4.729  0.000                            
## PHSFUNC -6.474 -1.156 -2.090  0.948  0.000                     
## ROLPHYS -4.326  0.791  2.692  0.861  3.310  0.000              
## BODPAIN -1.739  1.120 -1.049 -0.509  2.006  0.086  0.000       
## GENHLTH  3.183  2.327  0.754  4.582 -0.998 -4.531 -1.187  0.000
## 
## $summary
##                           cov
## srmr                    0.070
## srmr.se                 0.005
## srmr.exactfit.z         8.834
## srmr.exactfit.pvalue    0.000
## usrmr                   0.065
## usrmr.se                0.009
## usrmr.ci.lower          0.051
## usrmr.ci.upper          0.080
## usrmr.closefit.h0.value 0.050
## usrmr.closefit.z        1.769
## usrmr.closefit.pvalue   0.038

Modification Indices

The values below indicate the drop in chi-square (i.e., model improvement) anticipated if parameters currently fixed at 0.0 are freed to be estimated. Note that =~ indicates a crossloading of a subscale on its unintended factor, and ~~ indicates a correlation between the residual variances of two subscales. In this example they help us pinpoint why the fit of the hypothesized model was not great. GENHLTH “wants” to load on MH, and the correlation between MENHLTH & PHYSFUNC, which load on different factors, “wants” to be stronger than the model specifies.

modificationIndices(fit, sort. = TRUE, minimum.value = 3)

##        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
## 23      MH =~ GENHLTH 32.527  4.270   4.270    1.062    1.062
## 31 MENHLTH ~~ PHSFUNC 26.178 -4.306  -4.306   -0.395   -0.395
## 24      PH =~ MENHLTH 25.898 -3.608  -3.608   -0.816   -0.816
## 20      MH =~ PHSFUNC 21.195 -4.291  -4.291   -0.873   -0.873
## 29 MENHLTH ~~ ROLEEMO 20.032  2.690   2.690    0.385    0.385
## 41 ROLEEMO ~~   VITAL 17.170 -1.890  -1.890   -0.343   -0.343
## 54 ROLPHYS ~~ GENHLTH 16.933 -3.053  -3.053   -0.360   -0.360
## 50 PHSFUNC ~~ ROLPHYS 14.352  3.468   3.468    0.346    0.346
## 32 MENHLTH ~~ ROLPHYS 13.047 -2.280  -2.280   -0.296   -0.296
## 34 MENHLTH ~~ GENHLTH 12.376  2.468   2.468    0.267    0.267
## 43 ROLEEMO ~~ ROLPHYS 11.784  1.741   1.741    0.271    0.271
## 49   VITAL ~~ GENHLTH  9.706  1.689   1.689    0.232    0.232
## 27      PH =~   VITAL  8.440  1.558   1.558    0.467    0.467
## 51 PHSFUNC ~~ BODPAIN  4.492  1.180   1.180    0.167    0.167
## 46   VITAL ~~ PHSFUNC  3.711  1.252   1.252    0.146    0.146

Reliability Analysis

The semTools package has a function that reads the lavaan output and returns the reliability estimates for each latent variable. Setting the options ‘tau.eq=T’ and ‘obs.var=T’ produces equivalent to Cronbach’s reliability coefficient alpha. Setting ‘tau.eq=F’ produces the omega coefficient of reliability.

relests<-semTools::compRelSEM(fit, tau.eq=T, obs.var=T, return.total=T)

The reliability of MH is 0.8138296.
The reliability of PH is 0.7192983.
The overall reliability of the measurement model is 0.846431.

---