1 Introduction
1.1 Background
This entry examines the mathematical models of exploratory factor analysis (EFA) and its applications to the initial response in milliseconds across a series of choice tasks. The data on initial response times were extracted from the 2016 predictive modeling competition in health preference research (Jakubczyk, Craig, et al. 2017). We used the same data set to provide examples of calculating statistics in the previous entries. Please see the entry, Correlations between initial response times, as well as the entry, Illustrating time to initial response in choice tasks (Okubo and Craig, 2022).
In the previous entries, we have been analyzing statistics of observed variables, the initial response times, such as covariance, correlations, mean, and IQR. On the other hand, in this entry, we are attempting to figure out latent relations among dependent variables, the initial response times, through EFA. EFA models are possible to indicate interdependencies among dependent variables as well as latent factors that can explain the variance of dependent variables. As a result, we can simplify the analysis of data sets by reducing the number of variables by utilizing factors that explain the variance of dependent variables. FA models are different from regression models in the following sense: FA models are utilized in the process of data reduction by explaining dependent variables by common latent factors; on the other hand, regression models examine how a dependent variable could be explained by independent variables. FA models are also different from principal component analyses (PCA) in the following sense: FA is based on a model that assumes the existence of unobservable factors and error terms while PCA is not based on a model with them.
1.2 Load libraries and source files
Notes: Change the working directory to the location of the source files on your computer. To do so, you need to replace the inside of setwd() with the location of data file on your computer.
setwd("C:\\Users\\aaa\\OneDrive\\USF\\Dr. Craig\\entry 3")
# Please change the inside of quotes of 'setwd("")' to set your working directory.
#library(knitr) # This is for 'kable,' which makes well-organized tables.
library(tidyverse) # This is for 'read_csv.'
library(magrittr) # This is for '%>%.
library(reshape2) # This is for 'dcast.'
library(ggcorrplot) # This is for 'ggcorrplot.'
library(psych) # For factor analysis 'fa' and 'factanal'
library(GPArotation) # For rotation
library(gridExtra) # This is for 'grid.arrange'
data1 <- read_csv("resp1wave1_220723.csv")
# By 'read_csv("")' you can load the data on initial response times.We reorganize the loaded data, data1, which is a long format, into a wide format data, data2, by using the command, “dcast.”
data2 <- dcast(data1, survey_id ~ task)
data2 <- data2[,-1]
colnames(data2) <- paste0("Task", 1:max(data1$task))
rownames(data2) <- paste0("ID ", 1:max(data1$survey_id))2 Math model of exploratory factor analysis
\[ \begin{aligned} & x_1-\mu_1=\lambda_{11} f_1+\lambda_{12} f_2+\ldots+\lambda_{1m} f_{m}+\varepsilon_1 \\ & x_2-\mu_2=\lambda_{21} f_1+\lambda_{22} f_2+\cdots+\lambda_{2 n} f_m+\varepsilon_2 \\ \vdots\\ & x_p-\mu_p=\lambda_{p1} f_1+\lambda_{p2} f_2+\cdots+\lambda_{p m} f_m+\varepsilon_p \\ \\ \end{aligned} \tag{} \]
\[ \begin{equation} \vec{x}-\vec{\mu}=\Lambda \vec{f}+\vec{\varepsilon} \end{equation} \tag{} \]
where, $$ \[\begin{aligned} & \vec{x}=\left[\begin{array}{c} x_1 \\ \vdots \\ x_p \end{array}\right], \ \vec{\mu} =\left[\begin{array}{c} \mu_1 \\ \vdots \\ \mu_p \end{array}\right], \ \Lambda=\left[\begin{array}{cccc} \lambda_{11} & \lambda_{1 2} & \cdots & \lambda_{1 m} \\ \lambda_{21} & \lambda_{22} & \ldots & \vdots \\ \vdots & & & \vdots\\ \lambda_{p1} & \cdots & \cdots & \lambda_{pm} \end{array}\right], \ \vec{\varepsilon}=\left[\begin{array}{c} \varepsilon_1 \\ \vdots \\ \varepsilon_p \end{array}\right], \ \vec{f}=\left[\begin{array}{c} f_1 \\ \vdots \\ f_m \end{array}\right] \\ & \end{aligned} \tag{}\]$$ Actually, \(x_1, x_2, \cdots, x_p\) are vectors that store observations of each individual (sample size \(N\)).
Assumptions for the factor analysis models: \[ \begin{equation} \begin{aligned} & \cdot \mathbb{E}[\vec{f}]=\overrightarrow{0} \\ & \cdot \mathbb{E}[\vec{\varepsilon}]=\overrightarrow{0} \\ & \cdot \operatorname{Cov}(\vec{\varepsilon})=\left[\begin{array}{cccc} \psi_1 & 0 & \cdots & 0 \\ 0 & \psi_2 & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & \psi_p \end{array}\right]=\operatorname{Diag}\left(\psi_1\ \psi_2\ \ldots\ \psi_p\right)=\Psi \\ & \cdot \operatorname{Cov}(\vec{\varepsilon}, \vec{f})=0 \\ & \cdot \operatorname{Cov}(\vec{f})=I \end{aligned} \end{equation} \tag{} \] Specifically, the condition, \(\operatorname{Cov}(\vec{f})=I\), is needed for orthogonal factor models.
$$ \[\begin{equation} \begin{aligned} \operatorname{Cov}(\vec{x}) & =\mathbb{E}\left[(\vec{x}-\vec{\mu})(\vec{x}-\vec{\mu})^{\prime}\right]=\mathbb{E}\left[(\Lambda \vec{f}+\vec{\varepsilon})(\Lambda \vec{f}+\vec{\varepsilon})^{\prime}\right] \\ & =\mathbb{E}\left[(\Lambda \vec{f}+\vec{\varepsilon})\left(\vec{f}^{\prime} \Lambda^{\prime}+\vec{\varepsilon}^{\prime}\right)\right] \\ & =\mathbb{E}\left[\Lambda \vec{f} \vec{f}^{\prime} \Lambda^{\prime}+\Lambda \vec{f} \vec{\varepsilon}^{\prime}+\vec{\varepsilon} \vec{f}^{\prime} \Lambda^{\prime}+\vec{\varepsilon} \vec{\varepsilon}^{\prime}\right] \\ & =\Lambda \mathbb{E}\left[\vec{f} \vec{f}^{\prime}\right] \Lambda^{\prime}+\Lambda \mathbb{E}\left[\vec{f}^{\prime} \vec{\varepsilon}^{\prime}\right]+\mathbb{E}\left[\vec{\varepsilon} \vec{f}^{\prime}\right] \Lambda^{\prime}+\mathbb{E}\left[\vec{\varepsilon} \vec{\varepsilon}^{\prime}\right] \\ & =\Lambda I \Lambda^{\prime}+\Lambda \operatorname{Cov}(\vec{\varepsilon}, \vec{f})+\operatorname{Cov}(\vec{\varepsilon}, \vec{f}) \Lambda^{\prime}+\operatorname{Cov}(\vec{\varepsilon}) \\ & =\Lambda \Lambda^{\prime}+0+0+\Psi\\ & =\Lambda \Lambda^{\prime}+\Psi \end{aligned} \end{equation}\] \[ \] \[\begin{aligned} \operatorname{Cov}(\vec{x}, \vec{f}) & =\operatorname{Cov}(\Lambda \vec{f}+\vec{\varepsilon}, \vec{f}) \\ & =\mathbb{E}\left[(\Lambda \vec{f}+\vec{\varepsilon}) \vec{f}^{\prime}\right]=\mathbb{E}\left[\Lambda \vec{f} \vec{f}^{\prime}+\vec{\varepsilon} \vec{f}^{\prime}\right] \\ & =\Lambda \mathbb{E}\left[\vec{f} \vec{f}^{\prime}\right]+\mathbb{E}\left[\vec{\varepsilon} \vec{f}{ }^{\prime}\right] \\ & =\Lambda \operatorname{Cov}(\vec{f})+\operatorname{Cov}(\vec{\varepsilon}, \vec{f}) \\ & =\Lambda I +0 \\ & =\Lambda \end{aligned} \tag{}\]$$
\[ \begin{aligned} & \operatorname{Cov}(\vec{x})=\left[\begin{array}{cccc} \sigma_{11} & \sigma_{12} & & \sigma_{1 p} \\ \sigma_{21} & \sigma_{22} & \cdots & \sigma_{2 p} \\ \vdots & & & \vdots \\ \sigma_{p 1} & \sigma_{p 2} & \cdots & \sigma_{p p} \end{array}\right]=\Sigma \\ \end{aligned} \tag{} \] $$ \[\begin{aligned} & \Lambda \Lambda^{\prime}=\left[\begin{array}{ccc} \lambda_{11} & \cdots & \lambda_{1 m} \\ \vdots & & \vdots \\ \lambda_{p 1} & \cdots & \lambda_{p m} \end{array}\right]\left[\begin{array}{ccc} \lambda_{11} & \cdots & \lambda_{p 1} \\ \vdots & & \vdots \\ \lambda_1 m & \cdots & \lambda_{p m} \end{array}\right] \\ \\ & =\left[\begin{array}{cccc} \sum_{k=1}^m \lambda_{1 k}^2 & \sum_{k=1}^m \lambda_{1 k} \lambda_{2 k} & \cdots & \sum_{k=1}^m \lambda_{1 k} \lambda_{p k} \\ \sum_{k=1}^m \lambda_{2 k} \lambda_{1 k} & \sum_{k=1}^m \lambda_{2 k}^2 & \cdots &\vdots \\ \vdots & & &\vdots\\ \sum_{k=1}^m \lambda_{p k} \lambda_{1 k} & \cdots & \cdots & \sum_{k=1}^m \lambda_{p k}^2 \\ \end{array}\right] \\ & \end{aligned} \tag{}\]$$
\[ \begin{aligned} \sum & =\Lambda\Lambda^{\prime}+\Psi \\ \end{aligned} \tag{} \]
$$ \[\begin{aligned} \sigma_{ij} &= \sum_{k=1}^{m} \lambda_{i k} \lambda_{j k}\\ \sigma_{ii}&= \sum_{k=1}^m \lambda_{i k}^2+\psi_i \\ & =\underbrace{\left(\lambda{i 1}^2+\lambda_{i 2}^2+\cdots+\lambda_{i m}^2\right)}_{\text {Communality }}\ +\underbrace{\psi_i}_{Specific\ Variance}\\ &=h_i^2 + \psi_i \end{aligned} \tag{}\]$$
3 Scree plot
3.1 KMO
KMO(data2) # check if there exists at least one latent variable (factor). ## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data2)
## Overall MSA = 0.91
## MSA for each item =
## Task1 Task2 Task3 Task4 Task5 Task6 Task7 Task8 Task9 Task10 Task11
## 0.90 0.91 0.89 0.93 0.91 0.93 0.92 0.90 0.93 0.93 0.93
## Task12 Task13 Task14 Task15 Task16 Task17 Task18 Task19 Task20
## 0.90 0.89 0.89 0.87 0.92 0.89 0.88 0.87 0.89# Ideally, overall MSA > 0.6.3.2 Scree plot
VSS.scree(data2, main = "Graph: Scree Plot" ) #Create a scree plot. 
fa.parallel(data2, main = "Graph:Parallel Analysis Scree Plot") # Create a scree plot by parallel analysis. 
## Parallel analysis suggests that the number of factors = 5 and the number of components = 34 Loading matrix
##Use 3 factors ----
### use "fa" ---
f1.3.1 <- fa(data2, nfactors = 3, fm = "pa", rotate="none") # factor analysis by 'fa'
f1.3.2 <- fa(data2, nfactors = 3, fm = "pa", rotate ="varimax")
f1.3.3 <- fa(data2, nfactors = 3, fm = "pa", rotate = "promax")
f1.3.4 <- fa(data2, nfactors = 3, fm = "pa", rotate = "oblimin")
f2.3.1 <- fa(data2, nfactors = 3, fm = "ml", rotate="none") # factor analysis by 'fa'
f2.3.2 <- fa(data2, nfactors = 3, fm = "ml", rotate ="varimax")
f2.3.3 <- fa(data2, nfactors = 3, fm = "ml", rotate = "promax")
f2.3.4 <- fa(data2, nfactors = 3, fm = "ml", rotate = "oblimin")
loading1.3.1 <- cbind(round(f1.3.1$loadings, digits = 4)) %>% data.frame()
loading1.3.2 <- cbind(round(f1.3.2$loadings, digits = 4)) %>% data.frame()
loading1.3.3 <- cbind(round(f1.3.3$loadings, digits = 4)) %>% data.frame()
loading1.3.4 <- cbind(round(f1.3.4$loadings, digits = 4)) %>% data.frame()
loading2.3.1 <- cbind(round(f2.3.1$loadings, digits = 4)) %>% data.frame()
loading2.3.2 <- cbind(round(f2.3.2$loadings, digits = 4)) %>% data.frame()
loading2.3.3 <- cbind(round(f2.3.3$loadings, digits = 4)) %>% data.frame()
loading2.3.4 <- cbind(round(f2.3.4$loadings, digits = 4)) %>% data.frame()r1 <- as.matrix(loading1.3.1)
fig1 <-ggcorrplot(r1, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Principal Factor Solution and No Rotation")
r2 <- as.matrix(loading1.3.2)
fig2 <- ggcorrplot(r2, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Principal Factor Solution and Varimax Rotation")
r3 <- as.matrix(loading1.3.3)
fig3 <-ggcorrplot(r3, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Principal Factor Solution and Promax Rotation")
r4 <- as.matrix(loading1.3.4)
fig4 <-ggcorrplot(r4, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Principal Factor Solution and Oblimin Rotation")gridExtra::grid.arrange(fig1, fig2, fig3, fig4, nrow = 4)
ml1 <- as.matrix(loading2.3.1)
figml1 <-ggcorrplot(ml1, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Maximum Likelihood and No Rotation")
ml2 <- as.matrix(loading2.3.2)
figml2 <- ggcorrplot(ml2, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Maximum Likelihood and Varimax Rotation")
ml3 <- as.matrix(loading2.3.3)
figml3 <-ggcorrplot(ml3, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Maximum Likelihood and Promax Rotation")
ml4 <- as.matrix(loading2.3.4)
figml4 <-ggcorrplot(ml4, hc.order =F, lab =TRUE, tl.cex = 7, lab_size = 2.5, ) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1),
name = "Factor loadings", space = "Lab") +
labs(title = "Figure: Loading Matrix for Maximum Likelihood and Oblimin Rotation")gridExtra::grid.arrange(figml1, figml2, figml3, figml4, nrow = 4)
5 Illustrate EFA
par(mfrow = c(1,2))
chart1.3.1 <- fa.diagram(f1.3.1, main = "Chart:EFA with Principal Factor Solution and No Rotation")
fa.diagram(f2.3.1, main = "Chart:EFA with Maximum Likelihood and No Rotation")
chart1.3.2 <- fa.diagram(f1.3.2, main = "Chart:EFA with Principal Factor Solution and Varimax Rotation")
fa.diagram(f2.3.2, main = "Chart:EFA with Maximum Likelihood and Varimax Rotation")
chart1.3.3 <- fa.diagram(f1.3.3, main = "Chart:EFA with Principal Factor Solution and Promax Rotation")
fa.diagram(f2.3.3, main = "Chart:EFA with Maximum Likelihood and Promax Rotation")
chart1.3.4 <- fa.diagram(f1.3.4, main = "Chart:EFA with Principal Factor Solution and Oblimin Rotation")
fa.diagram(f2.3.4, main = "Chart:EFA with Maximum Likelihood and Oblimin Rotation")
6 Reference
Jakubczyk, M., Craig, B. M., Barra, M., Groothuis-Oudshoorn, C. G. M., Hartman, J. D., Huynh, E., Ramos-Goñi, J. M., Stolk, E. A., & Rand, K. (2017). Choice defines value: A predictive modeling competition in Health Preference Research. Value in Health, 21(2), 229–238. https://doi.org/10.1016/j.jval.2017.09.016
Mukhopadhyay, P. (2009). Multivariate Statistical Analysis. World Scientific.
Raykov, T., & Marcoulides, G. A. (2008). An Introduction to Applied Multivariate Analysis (1st ed.). Routledge.
7 How to cite this entry
Okubo, S. (yyyy, month dd). Title. R4HPR.