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

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

library(readxl)
library(dplyr)
library(writexl)
library(lmtest)
library(sandwich)
library(psych)

# =========================
# 1. 读取数据
# =========================
file_path <- "C:/Users/HP/Desktop/山青院学报/数据-英文.xlsx"

# 看一下 sheet 名
excel_sheets(file_path)

# 这里改成真实 sheet 名，通常是 data
df <- read_excel(file_path, sheet = "data")

# 查看变量名
names(df)

# =========================
# 2. 构造量表均值变量
# =========================
df <- df %>%
  mutate(
    # 职业规划三维
    career_explore = rowMeans(select(., C1:C4), na.rm = TRUE),
    career_clarity = rowMeans(select(., C5:C8), na.rm = TRUE),
    career_action  = rowMeans(select(., C9:C12), na.rm = TRUE),
    career_plan    = rowMeans(select(., C1:C12), na.rm = TRUE),
    
    # 学习参与三维
    cognitive_eng  = rowMeans(select(., S1:S4), na.rm = TRUE),
    behavior_eng   = rowMeans(select(., S5:S8), na.rm = TRUE),
    practical_eng  = rowMeans(select(., S9:S12), na.rm = TRUE),
    learn_engage   = rowMeans(select(., S1:S12), na.rm = TRUE),
    
    # 中介变量
    prof_identity  = rowMeans(select(., I1:I4), na.rm = TRUE),
    acad_efficacy  = rowMeans(select(., I5:I8), na.rm = TRUE)
  )

# =========================
# 3. 处理控制变量
# =========================
# 这里默认你的控制变量名称如下：
# gender, grade, major, type, hometown,
# father_edu, mother_edu, rank, cadre, career_course, internship

df <- df %>%
  mutate(
    gender        = as.factor(gender),
    grade         = as.factor(grade),
    major         = as.factor(major),
    type          = as.factor(type),
    hometown      = as.factor(hometown),
    father_edu    = as.factor(father_edu),
    mother_edu    = as.factor(mother_edu),
    rank          = as.factor(rank),
    cadre         = as.factor(cadre),
    career_course = as.factor(career_course),
    internship    = as.factor(internship)
  )

# =========================
# 4. 描述统计与相关分析
# =========================
desc_vars <- df %>%
  select(career_plan, career_explore, career_clarity, career_action,
         learn_engage, cognitive_eng, behavior_eng, practical_eng,
         prof_identity, acad_efficacy)

desc_table <- psych::describe(desc_vars)

cor_table <- round(cor(desc_vars, use = "pairwise.complete.obs"), 3)

# =========================
# 5. 定义稳健标准误回归输出函数
# =========================
reg_robust <- function(model) {
  coeftest(model, vcov = vcovHC(model, type = "HC3"))
}

extract_reg <- function(model, model_name) {
  ct <- coeftest(model, vcov = vcovHC(model, type = "HC3"))
  out <- data.frame(
    model = model_name,
    term = rownames(ct),
    estimate = ct[,1],
    std_error = ct[,2],
    t_value = ct[,3],
    p_value = ct[,4],
    row.names = NULL
  )
  return(out)
}

# =========================
# 6. 基准回归：职业规划 -> 学习参与
# =========================
m1 <- lm(
  learn_engage ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

summary(m1)
reg_robust(m1)

# =========================
# 7. 分维度回归
# =========================
m2 <- lm(
  cognitive_eng ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

m3 <- lm(
  behavior_eng ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

m4 <- lm(
  practical_eng ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

summary(m2); reg_robust(m2)
summary(m3); reg_robust(m3)
summary(m4); reg_robust(m4)

# =========================
# 8. 三个职业规划维度同时进入模型
# =========================
m5 <- lm(
  learn_engage ~ career_explore + career_clarity + career_action +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

summary(m5)
reg_robust(m5)

# =========================
# 9. 中介效应回归（分步）
# =========================
# 9.1 职业规划 -> 专业认同
m6 <- lm(
  prof_identity ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

# 9.2 职业规划 -> 学习自我效能感
m7 <- lm(
  acad_efficacy ~ career_plan +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

# 9.3 加入两个中介后的模型
m8 <- lm(
  learn_engage ~ career_plan + prof_identity + acad_efficacy +
    gender + grade + major + type + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

summary(m6); reg_robust(m6)
summary(m7); reg_robust(m7)
summary(m8); reg_robust(m8)

# =========================
# 10. 异质性/调节效应：本科类型交互项
# =========================
# 如果 type 是分类变量，先确认基准组
levels(df$type)

m9 <- lm(
  learn_engage ~ career_plan * type +
    gender + grade + major + hometown +
    father_edu + mother_edu + rank + cadre +
    career_course + internship,
  data = df
)

summary(m9)
reg_robust(m9)

# =========================
# 11. 导出回归结果
# =========================
reg_results <- bind_rows(
  extract_reg(m1, "M1_基准模型_学习参与"),
  extract_reg(m2, "M2_认知投入"),
  extract_reg(m3, "M3_行为投入"),
  extract_reg(m4, "M4_实践投入"),
  extract_reg(m5, "M5_职业规划三维"),
  extract_reg(m6, "M6_职业规划_专业认同"),
  extract_reg(m7, "M7_职业规划_学习自我效能感"),
  extract_reg(m8, "M8_中介模型"),
  extract_reg(m9, "M9_交互项模型")
)

model_fit <- data.frame(
  model = c("M1_基准模型_学习参与", "M2_认知投入", "M3_行为投入", "M4_实践投入",
            "M5_职业规划三维", "M6_职业规划_专业认同", "M7_职业规划_学习自我效能感",
            "M8_中介模型", "M9_交互项模型"),
  R2 = c(summary(m1)$r.squared,
         summary(m2)$r.squared,
         summary(m3)$r.squared,
         summary(m4)$r.squared,
         summary(m5)$r.squared,
         summary(m6)$r.squared,
         summary(m7)$r.squared,
         summary(m8)$r.squared,
         summary(m9)$r.squared),
  Adj_R2 = c(summary(m1)$adj.r.squared,
             summary(m2)$adj.r.squared,
             summary(m3)$adj.r.squared,
             summary(m4)$adj.r.squared,
             summary(m5)$adj.r.squared,
             summary(m6)$adj.r.squared,
             summary(m7)$adj.r.squared,
             summary(m8)$adj.r.squared,
             summary(m9)$adj.r.squared),
  N = nrow(df)
)

# =========================
# 12. 导出到桌面
# =========================
desktop_path <- "C:/Users/HP/Desktop"

write_xlsx(
  list(
    描述统计 = as.data.frame(desc_table),
    相关矩阵 = as.data.frame(cor_table),
    回归结果 = reg_results,
    模型拟合 = model_fit,
    含均值变量数据 = df
  ),
  file.path(desktop_path, "回归分析结果汇总.xlsx")
)

cat("\n全部回归分析完成，结果已导出到桌面。\n")