因子分析适合处理这样一类社会科学研究问题:问卷中有许多彼此相关的题项,但研究者希望将它们归并为少数几个具有理论含义的潜在维度。例如,多个关于学习投入、学业压力或组织认同的题项,可能分别反映更少数量的潜在因子。使用 R 完成因子分析时,真正困难的地方通常不在于运行一条函数,而在于确定哪些变量可以进入模型、提取几个因子、如何解释载荷,以及如何把分析结果转化为后续回归或结构模型中的变量。

下面以问卷题项为例,介绍一套可以写入论文方法部分的完整流程。示例中的变量名使用 q1q2 等占位符,实际分析时应替换为自己的题项名称。

一、先明确因子分析的研究目的

因子分析不是简单的“把变量压缩成几个变量”。在开始运行 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)))

如果题项采用五级或七级量表,应确认所有变量的编码方向一致。例如,有些题目分值越高表示认同程度越高,另一些反向题可能分值越高表示认同程度越低。反向题如果未处理,会降低相关性,甚至导致因子结构被错误分裂。

以五点量表为例,假设 q3q7 是反向题:

# 五点量表的反向计分
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 中重复运行分析,并让论文读者理解每一个因子结果是如何形成的。

声明:本站所有文章,如无特殊说明或标注,均为本站原创发布。任何个人或组织,在未征得本站同意时,禁止复制、盗用、采集、发布本站内容到任何网站、书籍等各类媒体平台。如若本站内容侵犯了原著者的合法权益,可联系我们进行处理。