R包lavaan的验证性因子分析（CFA）
图注：常见精神疾病的因子结构示例，似乎可以归纳为两个潜在的维度，即内在性和外在性疾病。
因子分析主要存在两种形式：
探索性因子分析（Exploratory Factor Analysis，EFA）无需事先获知潜在因子数量，即无需先验知识，需要提取的因子数量通过数据判断，并使用因子载荷了解因子和变量间的关系。EFA目标是确定变量和因子之间的最相关关系。
验证性因子分析（Confirmatory Factor Analysis，CFA）是结构方程模型（SEM）的一种方法，试图确定因子的数量以及观测变量在因子上的载荷是否符合基于既定理论的预期。CFA从先验知识出发，预先确定一组简单的因子结构，包括假定潜在的因子数量和组成、以及因子间的影响方式等，并通过拟合优度对假设结构进行数据相关性评估。
EFA和CFA的最主要区别在于，既定因子结构是否是由先验知识获得的。
lavaan的CFA方法
EFA已在前文作了简介，本篇简介R包lavaan的CFA方法。
数据集
lavaan包的内置数据集HolzingerSwineford1939，由来自两个不同学校（Pasteur和Grant-White）的七年级和八年级孩子的智力测验得分组成（变量X1-Xn），以及孩子的年龄、性别等信息。
library(lavaan)
#数据集，详情 ?HolzingerSwineford1939
hz <- HolzingerSwineford1939
head(hz)
已知测验（变量）X1、X2、X3的分数和视觉（visual）能力有关， X4、X5、X6的分数和写作（writing）能力有关，X7、X8、X9的分数和数学（maths）能力有关。
因此，这三类学习能力可视为代表该数据集智力测验得分的潜在因子。
CFA模型
接下来对该数据集执行CFA，看使用假定的三种学习能力是否能够有效代表原数据集中9种测验得分间的关系。
#因子和变量关系的书写格式，=~ 表示“由...表示”
hz.model <- '
visual =~ x1 + x2 + x3
writing =~ x4 + x5 + x6
maths =~ x7 + x8 + x9'
#执行 CFA 的一个简单方法，详情 ?cfa
hz.fit <- cfa(model = hz.model, data = hz)
summary(hz.fit, standardized = TRUE)
概要中主要包含如下部分，包含了拟合度指标和统计检验结果：
参数估计（Parameter Estimates），观测变量到潜在因子的标准化权重；
潜在变量（Latent Variables），估计的潜在因子的因子得分；
因子协方差（Covariances），潜在因子之间的协方差；
误差方差（Variances），每个观测变量的误差方差的估计值。
通过树状图评估因子和变量间的关系。
例如，边的粗细分别代表参数估计值或标准化的参数估计值。对于变量和因子，边越粗代表变量在因子上的载荷越高；对于因子，边越粗代表因子间的相关性越强。
#可视化因子和变量间的关系，详情 ?semPathslibrary(semPlot)
par(mfrow = c(1, 2))
semPaths(hz.fit, what = 'est', layout = 'tree')
semPaths(hz.fit, what = 'std', layout = 'tree')
拟合优度评估
评估CFA的整体拟合度，可通过fitmeasures()函数实现。
#模型拟合度，详情 ?fitmeasures
fitmeasures(hz.fit, c('cfi', 'rmsea', 'rmsea.ci.upper', 'bic'))
比较拟合指数（comparative fit index，CFI）介于0和1之间，根据一般经验，良好拟合模型的阈值通常要求CFI>0.9。
近似值的均方根误差（Root Mean Square Error of Approximation，RMSEA）值等于0.01、0.05和0.08分别表示优秀、良好和中等的拟合度， RMSEA<0.05通常用作合理拟合模型的临界值。
整体来看，本次CFA还不错，尽管拟合指数表明（特别是RMSEA值高于0.05），模型可以改进。
其它情况，拟合度较差时的可能原因：
变量的数量，尽管RMSEA趋势会随变量数量增加而变差，但其它拟合统计量一般会提升；
模型的复杂性，不同的统计量在简约模型中的改变程度不同；
样本数量，不同的统计量存在差异，有些会随样本数量增加而提升，但有些会变差；
数据的非正态性（趋于）会使拟合度变差，一般而言，CFA要求输入数据尽可能满足多元正态性假设。
模型优化
如果期望优化现有的CFA结果，首先确定需要修改哪些变量与哪些因子之间的潜在关联。
#检查修改索引
mf <- modificationindices(hz.fit)
mf <- mf[order(mf$mi, mf$epc, decreasing = TRUE), ] #按 MI 值排序
mf
lhs、op、rhs指示了推荐修改的因子与变量间关系。MI值越高，代表这种关系更有利于改善现有模型。
例如这里将最大MI的对应路径添加到模型中。
#追加 X9 与 visual 的关系hz.model.2 <- '
visual =~ x1 + x2 + x3 + x9
writing =~ x4 + x5 + x6
maths =~ x7 + x8 + x9'
hz.fit.2 <- cfa(model = hz.model.2, data = hz)
#summary(hz.fit.2, standardized = TRUE)
#semPaths(hz.fit.2, what = 'std', layout = 'tree')
fitmeasures(hz.fit.2, c('cfi', 'rmsea', 'rmsea.ci.upper', 'bic'))
与上述初始建立的CFA相比，CFI和RMSEA值都有所改善，如果期望想进一步优化模型，可尝试其它改进方法。
推荐一次优化步骤只尝试一种推荐的路径，因为每次修改因子和变量的关系对时，MI值都会改变，根据MI值排序后的顺序可能并非和先前的一致。
在模型优化过程中，切记注意不要强行建立不相干变量与因子间的潜在关联，因为MI可能会建议没有实质意义的路径。尽管它们在数值上、数学上统计出相关性，但也必须注意这种相关是否是真实合理的。
此外，也尽可能避免很多路径的出现，使模型产生过拟合问题。
参考资料