# =========================
# 1. 安装并加载需要的包
# =========================
packages <- c("readxl", "psych", "GPArotation", "dplyr", "writexl")

installed <- rownames(installed.packages())
for (p in packages) {
  if (!(p %in% installed)) install.packages(p)
}

library(readxl)
library(psych)
library(GPArotation)
library(dplyr)
library(writexl)

# =========================
# 2. 读取数据
# =========================
# 把这里改成你的文件路径
file_path <- "C:/Users/HP/Desktop/数据-英文.xlsx"

# 查看 sheet 名称
excel_sheets(file_path)

# 这里改成你真正用于分析的 sheet 名称
# 例如 "清洗后（分析用）" 或 "清洗后保留样本"
df <- read_excel(file_path, sheet = "data")

# 看一下列名
names(df)

# =========================
# 3. 如果你的列名还是中文题目，建议先手动确认题项位置
# =========================
# 下面这一步假设：
# 职业规划题项 = C1-C12
# 学习参与题项 = S1-S12
# 中介变量 = I1-I8
#
# 如果你的表里已经把列名改成 C1、C2……S1……I1……，可直接用下面代码
# 如果还没改列名，请先把列名改短，或者按列号提取

# ---------- 方案A：如果你已经改成短变量名 ----------
career_plan <- df %>% select(C1:C12)
learn_engage <- df %>% select(S1:S12)
identity_eff  <- df %>% select(I1:I8)

career_explore <- df %>% select(C1:C4)
career_clarity <- df %>% select(C5:C8)
career_action  <- df %>% select(C9:C12)

cognitive_eng  <- df %>% select(S1:S4)
behavior_eng   <- df %>% select(S5:S8)
practical_eng  <- df %>% select(S9:S12)

prof_identity  <- df %>% select(I1:I4)
acad_efficacy  <- df %>% select(I5:I8)

# ---------- 方案B：如果你还没改列名，可按列号提取 ----------
# 先用 names(df) 看清楚对应位置，再取消注释
# career_plan    <- df[, 12:23]
# learn_engage   <- df[, 24:35]
# identity_eff   <- df[, 36:43]
#
# career_explore <- df[, 12:15]
# career_clarity <- df[, 16:19]
# career_action  <- df[, 20:23]
#
# cognitive_eng  <- df[, 24:27]
# behavior_eng   <- df[, 28:31]
# practical_eng  <- df[, 32:35]
#
# prof_identity  <- df[, 36:39]
# acad_efficacy  <- df[, 40:43]

# =========================
# 4. 确保题项都是数值型
# =========================
to_numeric_df <- function(x) {
  x %>% mutate(across(everything(), ~ as.numeric(.)))
}

career_plan    <- to_numeric_df(career_plan)
learn_engage   <- to_numeric_df(learn_engage)
identity_eff   <- to_numeric_df(identity_eff)

career_explore <- to_numeric_df(career_explore)
career_clarity <- to_numeric_df(career_clarity)
career_action  <- to_numeric_df(career_action)

cognitive_eng  <- to_numeric_df(cognitive_eng)
behavior_eng   <- to_numeric_df(behavior_eng)
practical_eng  <- to_numeric_df(practical_eng)

prof_identity  <- to_numeric_df(prof_identity)
acad_efficacy  <- to_numeric_df(acad_efficacy)

# =========================
# 5. 定义一个函数：输出信度分析结果
# =========================
run_alpha <- function(data, scale_name) {
  cat("\n============================\n")
  cat("量表名称：", scale_name, "\n")
  cat("============================\n")
  
  a <- psych::alpha(data)
  
  cat("Cronbach's alpha =", round(a$total$raw_alpha, 3), "\n")
  cat("标准化 alpha =", round(a$total$std.alpha, 3), "\n\n")
  
  cat("题项分析（含 CITC）:\n")
  print(round(a$item.stats[, c("r.drop", "mean", "sd")], 3))
  
  return(a)
}

# =========================
# 6. 信度分析
# =========================
alpha_cp_total   <- run_alpha(career_plan, "职业规划总量表")
alpha_ce         <- run_alpha(career_explore, "职业探索")
alpha_cc         <- run_alpha(career_clarity, "职业目标清晰度")
alpha_ca         <- run_alpha(career_action, "职业规划执行力")

alpha_le_total   <- run_alpha(learn_engage, "学习参与总量表")
alpha_cog        <- run_alpha(cognitive_eng, "认知投入")
alpha_beh        <- run_alpha(behavior_eng, "行为投入")
alpha_pra        <- run_alpha(practical_eng, "实践投入")

alpha_pi         <- run_alpha(prof_identity, "专业认同")
alpha_ae         <- run_alpha(acad_efficacy, "学习自我效能感")

# =========================
# 7. KMO 和 Bartlett 检验
# =========================
run_validity_test <- function(data, scale_name) {
  cat("\n============================\n")
  cat("效度检验：", scale_name, "\n")
  cat("============================\n")
  
  cor_mat <- cor(data, use = "pairwise.complete.obs")
  
  kmo_res <- KMO(cor_mat)
  bart_res <- cortest.bartlett(cor_mat, n = nrow(data))
  
  cat("KMO =", round(kmo_res$MSA, 3), "\n")
  cat("Bartlett 球形检验:\n")
  print(bart_res)
  
  return(list(kmo = kmo_res, bartlett = bart_res))
}

valid_cp   <- run_validity_test(career_plan, "职业规划总量表")
valid_le   <- run_validity_test(learn_engage, "学习参与总量表")
valid_ie   <- run_validity_test(identity_eff, "中介变量量表（I1-I8）")
valid_all  <- run_validity_test(cbind(career_plan, learn_engage, identity_eff), "全量表（32题）")

# =========================
# 8. 探索性因子分析（EFA）
# =========================
# 先看建议提取因子数
fa.parallel(career_plan, fa = "fa", main = "职业规划平行分析")
fa.parallel(learn_engage, fa = "fa", main = "学习参与平行分析")
fa.parallel(identity_eff, fa = "fa", main = "I量表平行分析")

# 根据理论结构提取因子
efa_cp <- fa(career_plan, nfactors = 3, rotate = "oblimin", fm = "ml")
efa_le <- fa(learn_engage, nfactors = 3, rotate = "oblimin", fm = "ml")
efa_ie <- fa(identity_eff, nfactors = 2, rotate = "oblimin", fm = "ml")

cat("\n============================\n职业规划 EFA 结果\n============================\n")
print(efa_cp$loadings, cutoff = 0.30)

cat("\n============================\n学习参与 EFA 结果\n============================\n")
print(efa_le$loadings, cutoff = 0.30)

cat("\n============================\nI量表 EFA 结果\n============================\n")
print(efa_ie$loadings, cutoff = 0.30)

# =========================
# 9. 计算各维度均值（后续回归会用到）
# =========================
df_result <- df %>%
  mutate(
    career_explore_mean = rowMeans(career_explore, na.rm = TRUE),
    career_clarity_mean = rowMeans(career_clarity, na.rm = TRUE),
    career_action_mean  = rowMeans(career_action, na.rm = TRUE),
    career_plan_mean    = rowMeans(career_plan, na.rm = TRUE),
    
    cognitive_mean      = rowMeans(cognitive_eng, na.rm = TRUE),
    behavior_mean       = rowMeans(behavior_eng, na.rm = TRUE),
    practical_mean      = rowMeans(practical_eng, na.rm = TRUE),
    learn_engage_mean   = rowMeans(learn_engage, na.rm = TRUE),
    
    prof_identity_mean  = rowMeans(prof_identity, na.rm = TRUE),
    acad_efficacy_mean  = rowMeans(acad_efficacy, na.rm = TRUE)
  )

# =========================
# 10. 导出均值变量数据到桌面
# =========================
desktop_path <- "C:/Users/HP/Desktop"

write_xlsx(
  df_result,
  file.path(desktop_path, "数据_含量表均值变量.xlsx")
)

# =========================
# 11. 汇总主要信度结果，导出表格到桌面
# =========================
alpha_summary <- data.frame(
  Scale = c(
    "职业规划总量表", "职业探索", "职业目标清晰度", "职业规划执行力",
    "学习参与总量表", "认知投入", "行为投入", "实践投入",
    "专业认同", "学习自我效能感"
  ),
  Alpha = c(
    alpha_cp_total$total$raw_alpha,
    alpha_ce$total$raw_alpha,
    alpha_cc$total$raw_alpha,
    alpha_ca$total$raw_alpha,
    alpha_le_total$total$raw_alpha,
    alpha_cog$total$raw_alpha,
    alpha_beh$total$raw_alpha,
    alpha_pra$total$raw_alpha,
    alpha_pi$total$raw_alpha,
    alpha_ae$total$raw_alpha
  )
)

kmo_summary <- data.frame(
  Scale = c("职业规划总量表", "学习参与总量表", "中介变量量表(I1-I8)", "全量表(32题)"),
  KMO = c(
    valid_cp$kmo$MSA,
    valid_le$kmo$MSA,
    valid_ie$kmo$MSA,
    valid_all$kmo$MSA
  )
)

write_xlsx(
  list(
    alpha_summary = alpha_summary,
    kmo_summary = kmo_summary
  ),
  file.path(desktop_path, "信效度分析结果汇总.xlsx")
)

cat("\n全部分析完成，文件已导出到桌面。\n")