因子分析适合处理这样一类社会科学研究问题:问卷中有许多彼此相关的题项,但研究者希望将它们归并为少数几个具有理论含义的潜在维度。例如,多个关于学习投入、学业压力或组织认同的题项,可能分别反映更少数量的潜在因子。使用 R 完成因子分析时,真正困难的地方通常不在于运行一条函数,而在于确定哪些变量可以进入模型、提取几个因子、如何解释载荷,以及如何把分析结果转化为后续回归或结构模型中的变量。
下面以问卷题项为例,介绍一套可以写入论文方法部分的完整流程。示例中的变量名使用 q1、q2 等占位符,实际分析时应替换为自己的题项名称。
一、先明确因子分析的研究目的
因子分析不是简单的“把变量压缩成几个变量”。在开始运行 R 之前,应先回答两个问题。
第一,题项是否具有共同的潜在结构?如果每个题项都在测量完全不同的内容,强行进行因子分析并不能产生有意义的因子。
第二,因子结果将如何用于后续研究?如果因子分数要进入回归、相关分析或组间比较,就需要在提取因子时关注稳定性、解释性和可复现性,而不能只选择能够最大限度解释方差的方案。
社会科学问卷中的题项通常来自理论构念,因此因子数量不能完全交给统计软件决定。统计指标可以帮助判断合理范围,但最终方案应同时考虑理论预期、题项内容和载荷结果。
二、导入数据并检查变量编码
假设数据保存在一个数据文件中,且每一行代表一名受访者,每一列代表一个题项。首先导入数据,并检查数据的基本结构。
# 读取数据
dat <- read.csv("survey_data.csv",
stringsAsFactors = FALSE,
na.strings = c("", "NA", "N/A", "缺失"))
# 查看数据结构
str(dat)
dim(dat)
head(dat)
# 查看每个变量的缺失数量
colSums(is.na(dat))
因子分析所使用的变量应当是数值型。如果问卷答案被读成字符或因子,需要先转换为数值。需要特别小心的是,直接使用 as.numeric() 转换因子变量,可能得到内部编码,而不是题项原本的分值。
# 查看变量类型
sapply(dat, class)
# 如果原始变量是字符型,确认其内容确实为数值后再转换
items <- c("q1", "q2", "q3", "q4", "q5", "q6",
"q7", "q8", "q9", "q10")
dat[items] <- lapply(dat[items], function(x) {
as.numeric(trimws(x))
})
# 检查转换后的取值
lapply(dat[items], function(x) sort(unique(x)))
如果题项采用五级或七级量表,应确认所有变量的编码方向一致。例如,有些题目分值越高表示认同程度越高,另一些反向题可能分值越高表示认同程度越低。反向题如果未处理,会降低相关性,甚至导致因子结构被错误分裂。
以五点量表为例,假设 q3 和 q7 是反向题:
# 五点量表的反向计分
dat$q3_r <- 6 - dat$q3
dat$q7_r <- 6 - dat$q7
# 后续使用反向处理后的变量
items <- c("q1", "q2", "q3_r", "q4", "q5",
"q6", "q7_r", "q8", "q9", "q10")
反向计分公式必须根据量表范围调整。对于最低分为 min_score、最高分为 max_score 的量表,反向后的分数通常可按下式计算:
dat$reverse_item <- max_score + min_score - dat$original_item
不要在没有核对问卷编码说明的情况下直接反向处理题项。题目措辞中的否定表达,并不一定意味着需要反向计分,最终应以量表设计和编码规则为准。
三、筛选进入模型的题项
检查取值范围和异常编码
先检查每个题项的最小值、最大值、均值和缺失情况。
item_summary <- data.frame(
variable = items,
n = sapply(dat[items], function(x) sum(!is.na(x))),
missing = sapply(dat[items], function(x) sum(is.na(x))),
min = sapply(dat[items], function(x) min(x, na.rm = TRUE)),
max = sapply(dat[items], function(x) max(x, na.rm = TRUE)),
mean = sapply(dat[items], function(x) mean(x, na.rm = TRUE)),
sd = sapply(dat[items], function(x) sd(x, na.rm = TRUE))
)
item_summary
如果量表范围为一到五,却出现了九十九、负值或其他异常编码,应先将这些编码转换为缺失值,再进行分析。
invalid_codes <- c(99, 999, -9)
dat[items] <- lapply(dat[items], function(x) {
x[x %in% invalid_codes] <- NA
x
})
不能把异常编码当成普通分值参与相关矩阵计算,否则相关系数和因子载荷都可能失真。
检查题项之间的相关性
因子分析要求题项之间存在一定程度的相关。如果所有题项之间都几乎没有相关性,提取出的因子通常缺少解释意义。
cor_mat <- cor(dat[items],
use = "pairwise.complete.obs",
method = "pearson")
round(cor_mat, 2)
如果使用成对删除处理缺失值,不同相关系数可能基于不同的样本量。缺失较多时,应在报告中说明缺失处理方式,并检查相关矩阵是否仍然稳定。对于题项缺失较少的情况,也可以使用完整案例进行分析:
dat_complete <- dat[complete.cases(dat[items]), ]
cor_mat_complete <- cor(dat_complete[items],
method = "pearson")
round(cor_mat_complete, 2)
不能仅凭某一个相关系数决定删除题项。题项之间的相关性还要结合题目内容、理论归属和后续载荷共同判断。
四、评估数据是否适合因子分析
计算 KMO 指标
KMO 指标用于比较题项之间的普通相关与偏相关。整体 KMO 越高,说明题项之间更可能存在共同因子结构。下面使用 R 代码计算整体 KMO 和各题项的 KMO。
kmo_result <- function(R) {
R_inv <- solve(R)
# 由逆相关矩阵计算偏相关矩阵
P <- -R_inv / sqrt(outer(diag(R_inv), diag(R_inv)))
diag(P) <- 1
r2 <- R^2
p2 <- P^2
# 去除对角线
diag(r2) <- 0
diag(p2) <- 0
overall_kmo <- sum(r2) / (sum(r2) + sum(p2))
item_kmo <- rowSums(r2) / (rowSums(r2) + rowSums(p2))
list(
overall = overall_kmo,
item = item_kmo
)
}
kmo <- kmo_result(cor_mat_complete)
kmo$overall
kmo$item
KMO 较低时,先不要急于删除变量。应检查是否存在以下问题:题项内容本身不属于同一构念、反向题没有正确处理、样本量过少、变量取值几乎没有变化,或者相关矩阵存在严重异常。
计算 Bartlett 球形检验
Bartlett 球形检验用于检验相关矩阵是否显著不同于单位矩阵。其原假设是各题项之间不存在总体相关。
bartlett_test <- function(R, n) {
p <- ncol(R)
determinant_value <- det(R)
chi_square <- -(n - 1 - (2 * p + 5) / 6) *
log(determinant_value)
df <- p * (p - 1) / 2
p_value <- pchisq(chi_square, df = df, lower.tail = FALSE)
data.frame(
chi_square = chi_square,
df = df,
p_value = p_value
)
}
bartlett_test(
R = cor_mat_complete,
n = nrow(dat_complete)
)
如果相关矩阵接近单位矩阵,因子分析通常缺乏基础。相反,检验显著也不能单独证明因子分析一定合理,因为样本量较大时,很小的相关也可能达到显著水平。因此,KMO、相关矩阵、题项内容和后续因子载荷应结合判断。
五、确定因子数量
因子数量是因子分析中最需要研究者判断的环节之一。可以先比较理论预期和不同因子数量下的模型表现,再观察结果是否具有清晰解释。
使用特征值进行初步判断
eigen_result <- eigen(cor_mat_complete)
eigen_values <- eigen_result$values
data.frame(
factor = seq_along(eigen_values),
eigenvalue = eigen_values,
explained_percent = eigen_values / sum(eigen_values) * 100,
cumulative_percent = cumsum(eigen_values / sum(eigen_values) * 100)
)
特征值可以帮助观察相关矩阵中主要维度的数量,但不应机械地把特征值大于一的维度全部保留。该标准容易在题项较多时提取过多因子,也可能产生难以解释的结构。
比较不同因子数量的可解释性
假设理论上可能存在二到四个因子,可以分别运行模型:
fa_2 <- factanal(
x = dat_complete[items],
factors = 2,
rotation = "promax",
scores = "regression"
)
fa_3 <- factanal(
x = dat_complete[items],
factors = 3,
rotation = "promax",
scores = "regression"
)
fa_4 <- factanal(
x = dat_complete[items],
factors = 4,
rotation = "promax",
scores = "regression"
)
fa_2
fa_3
fa_4
比较时重点关注以下问题:
- 每个因子是否至少有若干个具有清晰载荷的题项;
- 同一个因子中的题项是否具有共同的理论含义;
- 是否存在一个因子只有一个题项或极少题项;
- 是否有大量题项在多个因子上同时具有较高载荷;
- 增加因子后,解释是否变得更清晰,还是只是把相近题项机械拆开。
如果多个因子之间理论上可能相关,应优先考虑斜交旋转,例如 promax。如果研究者有充分理由认为因子彼此独立,可以考虑正交旋转,例如 varimax。社会科学构念往往存在相关,因此不应默认所有因子完全独立。
六、运行因子分析并选择旋转方式
使用最大似然法提取因子
factanal() 可以使用最大似然方法提取因子。下面以三个因子为例:
fa_model <- factanal(
x = dat_complete[items],
factors = 3,
rotation = "promax",
scores = "regression"
)
print(fa_model,
cutoff = 0.30,
sort = TRUE)
cutoff = 0.30 只影响打印结果,不会改变模型本身。为了判断题项在不同因子上的真实表现,也可以直接查看载荷矩阵:
loadings_matrix <- as.matrix(fa_model$loadings)
round(loadings_matrix, 3)
载荷可以理解为题项与因子的关系强度。载荷绝对值较高,通常说明题项与该因子联系较紧密,但不能只根据某个统一数值机械删除题项。载荷的解释还需要结合题项内容、样本量、因子数量和理论预期。
比较正交旋转与斜交旋转
fa_varimax <- factanal(
x = dat_complete[items],
factors = 3,
rotation = "varimax",
scores = "regression"
)
fa_promax <- factanal(
x = dat_complete[items],
factors = 3,
rotation = "promax",
scores = "regression"
)
round(as.matrix(fa_varimax$loadings), 3)
round(as.matrix(fa_promax$loadings), 3)
如果使用斜交旋转,应同时查看因子之间的相关矩阵:
fa_promax$Phi
因子之间存在相关并不代表模型失败。相反,如果理论上相关的构念在斜交旋转后显示出一定相关,往往比强行使用正交旋转更符合研究对象。报告中应说明提取方法、旋转方法和选择该方法的理由。
七、读取并解释因子载荷
为了更清楚地识别每个题项主要归属的因子,可以整理载荷矩阵。
loading_df <- as.data.frame(unclass(fa_model$loadings))
loading_df$item <- rownames(loading_df)
loading_df <- loading_df[, c("item",
setdiff(names(loading_df), "item"))]
loading_df
还可以为每个题项标记绝对载荷最大的因子:
factor_columns <- setdiff(names(loading_df), "item")
loading_df$primary_factor <- factor_columns[
apply(abs(loading_df[factor_columns]), 1, which.max)
]
loading_df$primary_loading <- apply(
abs(loading_df[factor_columns]),
1,
max
)
loading_df
解释载荷时,至少要关注三种情况。
第一,题项在某个因子上的主要载荷较高,且在其他因子上的载荷较低。这类题项通常有较清晰的归属。
第二,题项在两个或多个因子上的载荷都较高。这说明题项可能同时反映多个构念,也可能存在题目表述过于宽泛、反向题处理错误或因子数量设定不合适等问题。此时不能只看统计结果,还要回到题目内容判断是否保留。
第三,题项在所有因子上的载荷都较低。这类题项可能没有充分反映共同构念,但也应先检查数据编码、缺失处理和量表方向,再考虑删除。
因子命名应根据同一因子下题项的共同含义完成,而不是根据软件自动生成的编号决定。因子编号本身没有实质意义,旋转或重新运行模型后,因子编号可能改变。
例如,如果某个因子上的题项都围绕“主动参与课堂、投入学习时间和完成学习任务”,可以将其命名为“学习投入”。论文中应说明命名依据,并列出构成该因子的主要题项。
八、处理题项删除与模型重估
删除题项不能只因为它的载荷低。一个较稳妥的处理过程是:先发现问题,再核对题目内容,删除后重新估计模型,最后比较新旧结构是否更合理。
# 假设经过理论判断后,删除 q5 和 q8
items_reduced <- setdiff(items, c("q5", "q8"))
dat_reduced <- dat_complete[items_reduced]
fa_reduced <- factanal(
x = dat_reduced,
factors = 3,
rotation = "promax",
scores = "regression"
)
print(fa_reduced,
cutoff = 0.30,
sort = TRUE)
删除题项后,应重新检查相关矩阵、KMO、因子数量和载荷结构。不能先反复删除变量,直到输出结果看起来“漂亮”为止。过度根据当前样本调整题项,会增加结果对特定样本的依赖,也可能削弱量表的内容效度。
如果删除某个题项只是为了提高某个统计指标,却使理论构念的重要内容缺失,通常不应删除。论文中还应记录删除依据,而不是只报告最终保留下来的题项。
九、提取因子得分并构建后续变量
完成因子结构确认后,可以提取每位受访者的因子得分。使用 scores = "regression" 后,因子得分通常保存在模型对象的 scores 中。
factor_scores <- as.data.frame(fa_model$scores)
head(factor_scores)
summary(factor_scores)
将因子得分合并回原始数据:
dat_analysis <- cbind(
dat_complete,
factor_scores
)
head(dat_analysis)
为了让变量名更容易用于后续分析,可以重命名:
names(factor_scores) <- c(
"learning_engagement",
"academic_pressure",
"institutional_support"
)
dat_analysis <- cbind(
dat_complete,
factor_scores
)
因子得分通常是标准化形式,均值接近零。它们适合用于后续回归或相关分析,但解释时应明确说明:系数反映的是潜在因子得分变化与结果变量之间的关系,而不是原始问卷分值的直接变化。
除了模型估计的因子得分,也可以根据题项均值或总分构建量表得分。这种方法更容易解释,但前提是题项已经确认属于同一个因子,并且缺失处理规则已经确定。
engagement_items <- c("q1", "q2", "q4")
dat_analysis$engagement_mean <- rowMeans(
dat_analysis[engagement_items],
na.rm = TRUE
)
使用 rowMeans(..., na.rm = TRUE) 时要注意:如果某一行所有题项都是缺失,结果可能不符合预期。可以先设置最低有效题项数量。
row_mean_min <- function(x, min_valid = 2) {
valid_n <- rowSums(!is.na(x))
result <- rowMeans(x, na.rm = TRUE)
result[valid_n < min_valid] <- NA
result
}
dat_analysis$engagement_mean <- row_mean_min(
dat_analysis[engagement_items],
min_valid = 2
)
因子得分和题项均值不是同一个变量。前者由模型估计,考虑了题项之间的关系;后者直接反映题项答案的平均水平。研究者应根据研究目的和论文中的变量定义选择一种,并在全文保持一致。
十、将因子结果用于后续分析
如果后续需要把因子结果用于回归分析,可以直接使用已经生成的因子变量。下面的代码仅展示变量进入模型的方式,具体结果应根据自己的研究设计确定。
model <- lm(
outcome ~ learning_engagement +
academic_pressure +
institutional_support +
age +
gender,
data = dat_analysis
)
summary(model)
如果因子之间存在较高相关,应在解释回归结果时保持谨慎。高度相关的因子可能使回归系数不稳定,也可能导致单个因子的解释与零阶相关不一致。此时应回到因子相关矩阵和研究理论,确认这些因子是否确实代表可以区分的构念。
在构建后续变量之前,还应检查因子得分是否存在大量缺失、极端分布或与原始题项方向不一致。
sapply(dat_analysis[c(
"learning_engagement",
"academic_pressure",
"institutional_support"
)], function(x) {
c(
missing = sum(is.na(x)),
mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE)
)
})
十一、常见错误及排查方法
报错“相关矩阵不是正定矩阵”
这通常意味着题项之间存在完全重复、近乎重复,或者某些变量之间存在严重线性关系。可以先检查相关矩阵:
round(cor_mat_complete, 3)
重点查看是否有两个题项的相关系数接近一。如果存在,应核对这两个题项是否重复、是否被错误复制,或是否实际上测量了完全相同的内容。也要检查变量方差:
sapply(dat_complete[items], var, na.rm = TRUE)
方差为零或极低的题项没有足够的信息参与因子分析。
factanal() 无法收敛
无法收敛可能与因子数量过多、样本信息不足、相关矩阵异常或题项之间高度共线有关。可以先减少因子数量,检查数据和相关矩阵,再比较不同旋转方式。如果某些题项存在极端分布或几乎没有变化,也应优先处理这些变量。
不要把增加迭代次数当成唯一解决办法。即使模型最终收敛,如果因子数量和题项结构本身不合理,结果仍然可能无法解释。
结果中出现 Heywood 情况
如果某个题项的独特性接近零,或者模型估计出的共同度异常偏高,可能存在 Heywood 情况。这通常提示模型设定不合适、题项高度重复、因子数量不合理或数据规模不足。
fa_model$uniquenesses
发现异常时,应检查相关矩阵、因子数量和题项内容,必要时重新设定模型。不能简单把异常结果当作“题项质量很好”。
反向题载荷方向相反
因子载荷为负不一定意味着模型错误。它可能反映题项方向与其他题项相反,也可能是反向计分未完成。先核对原始问卷编码和反向处理过程,再决定是否需要调整。
如果只是整个因子的方向相反,而同一因子内各题项关系一致,可以在解释时统一方向;如果同一因子内部载荷正负混杂,则应重点检查题项编码和因子结构。
缺失值导致样本量变化
complete.cases() 会删除任一分析题项缺失的受访者,题项较多时可能损失较多样本。pairwise.complete.obs 虽然可以保留更多相关系数,但不同相关系数的样本量可能不同,相关矩阵有时会变得不稳定。
因此,论文中应明确报告缺失值处理方法、最终分析样本量,以及该处理方式可能带来的限制。不能在不同步骤中随意切换缺失处理方式而不记录。
十二、如何撰写可复现的分析报告
一份可复现的因子分析报告,应让读者知道数据如何进入模型、模型如何设定,以及最终因子如何被解释和使用。建议至少交代以下内容:
- 分析对象、题项数量和题项的计分范围;
- 反向题的处理方式;
- 缺失值和异常编码的处理方式;
- 题项进入或退出模型的依据;
- KMO 和 Bartlett 球形检验结果;
- 因子数量的判断依据;
- 因子提取方法和旋转方法;
- 主要题项在各因子上的载荷;
- 因子命名的理论依据;
- 因子得分或量表均值如何构建;
- 因子变量如何进入后续分析。
为了保留分析过程,可以将关键代码和模型对象保存下来:
analysis_record <- list(
items = items,
correlation_matrix = cor_mat_complete,
eigenvalues = eigen_values,
factor_model = fa_model,
loadings = as.matrix(fa_model$loadings),
factor_scores = fa_model$scores
)
saveRDS(analysis_record, file = "factor_analysis_record.rds")
也可以保存整理后的载荷表:
write.csv(
loading_df,
file = "factor_loadings.csv",
row.names = FALSE
)
论文中不必逐行展示全部代码,但应保留能够复现主要结果的脚本、变量筛选记录和模型设定。尤其要避免只报告最终因子名称,却不说明题项删除和因子数量选择过程。
因子分析的最终目标不是得到一张载荷表,而是建立一个有理论依据、统计结构清晰、能够用于后续研究的变量体系。只要数据编码、题项筛选、因子数量、旋转方式和变量构建都被完整记录,研究者就能在 R 中重复运行分析,并让论文读者理解每一个因子结果是如何形成的。


