R实现DEA-Tobit 两阶段分析流程

代码解读 (中文)

这段 R 代码执行了一个完整的 DEA-Tobit 两阶段分析流程，主要用于评估一组决策单元（DMU，例如公司、地区等）的效率，并探究影响这些效率值的因素。

主要步骤：

1. 环境准备与数据加载:

   • 加载必要的 R 包，涵盖 DEA 计算 (deaR)、数据处理 (dplyr, magrittr)、文件读写 (readxl, openxlsx)、可视化 (ggplot2, GGally)、字体支持 (showtext)、Word/Excel 输出 (officer, flextable) 以及 Tobit 回归 (AER)。
   • 设置中文字体 (showtext, font_add) 以确保图表中的中文能正确显示。
   • 从 Excel 文件 (“Data/数据模板.xlsx”) 读取数据。

2. 第一阶段：DEA 效率测算

   • 提取数据中的所有年份。
   • 初始化一个列表 (efficiency_results) 来存储不同 DEA 模型计算出的效率值。
   • 逐年循环计算效率：
     • 对每一年份的数据：
       • 定义投入 (inputs) 和产出 (outputs) 指标。
       • 使用 make_deadata 准备 DEA 数据格式。
       • 计算多种 DEA 模型的效率值：
         • CRS (固定规模报酬)
         • VRS (可变规模报酬)
         • SBM (非径向、非角度的松弛测量模型)
         • Super (超效率模型，用于区分效率值为1的DMU)
       • 计算规模效率 (Scale) = CRS 效率 / VRS 效率。
       • 将计算得到的各模型效率值（保留3位小数）与年份、DMU标识符一起存储到 efficiency_results 列表中。
   • 合并结果：将所有年份、所有模型的效率结果纵向合并 (rbind) 成一个大数据框 (all_efficiencies)，并添加一列 (Model) 标明效率类型。
   • 数据整合：将计算出的效率值 (all_efficiencies) 与原始输入数据 (data) 根据年份和 DMU 合并 (merge) 成 merged_data，为第二阶段做准备。

3. 第二阶段：Tobit 回归分析

   • (可选) 探索性分析：使用 ggpairs 绘制散点图矩阵，检查解释变量之间的相关性。
   • (可选/注释掉) R方计算：代码中有计算 R² 的注释，但未执行。
   • 执行 Tobit 回归：
     • 对每种 DEA 模型 (CRS, VRS, Scale, SBM, Super) 的效率值作为因变量：
       • 筛选出对应模型的数据 (subset)。
       • 使用 tobit 函数进行回归分析。因变量是 Efficiency，自变量是 解释变量1, 解释变量2, Input3。由于效率值通常在 0 到 1 之间（或 SBM/Super 可能大于1，但这里设定了 right=1，可能需要根据实际情况调整），Tobit 模型适合处理这种受限因变量。
       • 提取回归结果（系数、标准误、z 值、p 值）并存储在 tobit_results 列表中。
   • 打印回归结果：在控制台输出 tobit_results。

4. 结果导出与可视化

   • 导出到 Excel：
     • 创建一个新的 Excel 工作簿 (createWorkbook)。
     • 将每种模型下的 Tobit 回归结果写入到对应名称的工作表 (addWorksheet, writeData)。
     • 保存为 Excel 文件 (“Table/DEA_Tobit_Results.xlsx”)。
   • (注释) 三线表：代码中有创建三线表的注释，但未实现具体代码。
   • **打印回归结果 (重复)**：再次打印 tobit_results。
   • 模型诊断：绘制残差图：
     • 以 VRS 模型的 Tobit 回归为例。
     • 提取模型的残差 (residuals) 和拟合值 (fitted)。
     • 使用 ggplot2 绘制残差与拟合值的散点图，并添加平滑曲线，用于检查模型设定是否存在问题（如异方差）。
   • 计算并绘制边际效应图：
     • 定义一个函数 calculate_marginal_effects 来提取 Tobit 模型的系数（在此模型中近似于边际效应）。
     • 对每个模型再次运行 Tobit（代码中重复运行了），并计算指定自变量 (解释变量1, 解释变量2, Input3) 的边际效应。
     • 合并所有模型的边际效应结果。
     • 使用 ggplot2 绘制分组条形图，可视化不同解释变量在不同 DEA 模型效率下的边际影响。

总结: 该脚本自动化了 DEA 效率测算和后续使用 Tobit 模型分析效率影响因素的标准流程，并包含了结果输出和关键的可视化步骤（模型诊断和效应展示）。

# 加载相关的R包----
library(deaR)      # 用于数据包络分析（DEA）的R包
library(magrittr)  # 提供管道操作符 (%>%)，使代码更简洁
library(openxlsx)  # 用于读写Excel文件的R包
library(ggplot2)   # 用于数据可视化的强大R包
library(showtext)  # 用于在图形中添加和管理字体的R包
library(readxl)    # 用于读取Excel文件的R包
library(flextable) # 用于创建可灵活调整的表格的R包
library(officer)   # 用于操作Word和PowerPoint文档的R包
library(AER)       # 用于应用经济学研究的R包，包括Tobit回归等模型
library(GGally)    # 提供ggplot2扩展，用于绘制矩阵图和其他多变量图表
library(dplyr)


# 设置字体和中文支持
showtext_auto(enable = TRUE, record = TRUE)
font_add(family = "苹方",regular = "/System/Library/Fonts/PingFang.ttc")
#font_add(family = "黑体",regular = "/System/Library/Fonts/STHeiti Medium.ttc")

# 读取数据文件
data <- read_excel("Data/数据模板.xlsx")
# 提取年份列表
years <- unique(data$Year)

# 初始化存储所有年份的效率值的数据框
efficiency_results <- list(crs = data.frame(), vrs = data.frame(), scale = data.frame(), sbm = data.frame(), super = data.frame())

#一、循环计算每年的 DEA 效率-----
for (year in years) {
  # 提取该年的数据
  data_year <- subset(data, Year == year)
    # 定义投入和产出指标
  inputs <- c("Input1", "Input2", "Input3")
  outputs <- c("Output1", "Output2")
    # 创建 DEA 数据
  dea_data <- make_deadata(datadea = data_year, 
                           inputs = inputs,
                           outputs = outputs)
 # 构建 DEA 模型
  models <- list(
    crs = model_basic(dea_data, rts = "crs", orientation = "io"),
    vrs = model_basic(dea_data, rts = "vrs", orientation = "io"),
    sbm = model_sbmeff(dea_data, orientation = "io"),
    super = model_supereff(dea_data, rts = "crs", orientation = "io")
  )
  #计算规模效率
  scale_eff <- efficiencies(models$crs) / efficiencies(models$vrs)
  #计算效率，并存储结果
  efficiency_results$crs <- rbind(efficiency_results$crs, data.frame(Year = year, DMU = data_year$DMU, Efficiency = round(efficiencies(models$crs), 3)))
  efficiency_results$vrs <- rbind(efficiency_results$vrs, data.frame(Year = year, DMU = data_year$DMU, Efficiency = round(efficiencies(models$vrs), 3)))
  efficiency_results$scale <- rbind(efficiency_results$scale, data.frame(Year = year, DMU = data_year$DMU, Efficiency = round(scale_eff, 3)))
  efficiency_results$sbm <- rbind(efficiency_results$sbm, data.frame(Year = year, DMU = data_year$DMU, Efficiency = round(efficiencies(models$sbm), 3)))
  efficiency_results$super <- rbind(efficiency_results$super, data.frame(Year = year, DMU = data_year$DMU, Efficiency = round(efficiencies(models$super), 3)))
}

#合并所有模型的效率结果
all_efficiencies <- rbind(
  efficiency_results$crs %>% mutate(Model = "CRS"),
  efficiency_results$vrs %>% mutate(Model = "VRS"),
  efficiency_results$scale %>% mutate(Model = "Scale"),
  efficiency_results$sbm %>% mutate(Model = "SBM"),
  efficiency_results$super %>% mutate(Model = "Super")
)
#将效率结果与原数据合并
merged_data <- merge(data, all_efficiencies, by = c("Year", "DMU"))

#二、进行 Tobit 回归分析---------
##1.解释变量的相关度分析（散点图）---------
ggpairs(data, columns = c("解释变量1", "Input2", "Input3","Output1","解释变量2"))
#其中值越接近 1 相关性越显著

##2.方差被解释量计算----
# 计算相关系数
#r <- cor(fitted_values, model_data$apt)
#r_squared <- r^2#方差被解释量
#print(paste("VRS R-squared:", r_squared))
##3.Tobit 计算----
tobit_results <- list()
models <- c("CRS", "VRS", "Scale", "SBM", "Super")
for (model in models) {
  model_data <- subset(merged_data, Model == model)
  tobit_model <- tobit(Efficiency ~ 解释变量1 + 解释变量2 + Input3, data = model_data, left = 0, right = 1)#使用方式：在这里输入你的解释变量名称即可
  
  summary_tobit <- summary(tobit_model)
  tobit_results[[model]] <- data.frame(
    Model = model,
    Variable = rownames(summary_tobit$coefficients),
    Coefficient = summary_tobit$coefficients[, "Estimate"],
    StdError = summary_tobit$coefficients[, "Std. Error"],
    ZValue = summary_tobit$coefficients[, "z value"],
    PValue = summary_tobit$coefficients[, "Pr(>|z|)"]
  )
}
###打印回归结果
print(tobit_results)

#三、结果导出----
##1.创建表格----
wb <- createWorkbook()
for (model in models) {
  addWorksheet(wb, model)
  writeData(wb, sheet = model, tobit_results[[model]])
}
# 保存 Excel 文件
saveWorkbook(wb, "Table/DEA_Tobit_Results.xlsx", overwrite = TRUE)
##2.创建三线表，排除截距和 Log(scale)-----
#四、打印回归结果----
print(tobit_results)

#五、作图----
##1.创建残差图，检验模型与数据匹配程度------
# 提取 VRS 模型的 Tobit 回归模型
vrs_model_data <- subset(merged_data, Model == "VRS")
vrs_tobit_model <- tobit(Efficiency ~ 解释变量1 + 解释变量2, data = vrs_model_data, left = 0, right = 1)

# 提取残差和拟合值
tobit_residuals <- residuals(vrs_tobit_model)
fitted_values <- fitted(vrs_tobit_model)

residuals_data <- data.frame(
  Fitted = fitted_values,
  Residuals = tobit_residuals
)

# 生成残差图
ggplot(residuals_data, aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_smooth(method = "loess", col = "red") +
  labs(title = "DEA规模效率-Tobit 模型的残差图",
       x = "拟合值",
       y = "残差") +
  theme_minimal(base_family = "苹方")





#计算模型间的边际效应----
calculate_marginal_effects <- function(model, variable_name, data) {
  # 获取模型的系数
  coefs <- summary(model)$coefficients
  # 计算边际效应（对每个变量的系数）
  marginal_effects <- data.frame(
    Variable = variable_name,
    Effect = coefs[variable_name, "Estimate"]
  )
  return(marginal_effects)
}

#提取Tobit模型中的边际效应
marginal_effects_list <- list()
for (model in models) {
  model_data <- subset(merged_data, Model == model)
  tobit_model <- tobit(Efficiency ~ 解释变量1 + 解释变量2 + Input3, data = model_data, left = 0, right = 1)
  
  # 计算每个变量的边际效应
  effect1 <- calculate_marginal_effects(tobit_model, "解释变量1", model_data)
  effect2 <- calculate_marginal_effects(tobit_model, "解释变量2", model_data)
  effect3 <- calculate_marginal_effects(tobit_model, "Input3", model_data)
  
  marginal_effects_list[[model]] <- rbind(effect1, effect2, effect3)
}

# 合并所有模型的边际效应结果
all_marginal_effects <- do.call(rbind, marginal_effects_list)
all_marginal_effects$Model <- rep(models, each = nrow(all_marginal_effects) / length(models))

# 绘制边际效应图
ggplot(all_marginal_effects, aes(x = Variable, y = Effect, fill = Model)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(title = "边际效应图", x = "解释变量", y = "边际效应") +
  scale_fill_brewer(palette = "Set1") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))