R实现Dea 基础效率+超效率计算及图表输出

library(deaR)
library(magrittr)
library(openxlsx)
library(ggplot2)
library(showtext)
library(readxl)
library(flextable)
library(officer)
library(reshape2)

# 设置字体和中文支持
font_add_google("Noto Sans SC", "noto")
showtext_auto(enable = TRUE)

# 读取数据文件----
data <- read_excel("Data/数据模板.xlsx")
# 自定义函数来格式化数字保留 3 位小数但不包含尾随零
format_number <- function(x) {
  sprintf("%.3f", x) %>%
    sub("\\.?0+$", "", .)
}

# 提取年份列表
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 <- data_year[, c("Input1", "Input2", "Input3")]
  outputs <- data_year[, c("Output1", "Output2")]
  
  dea_data <- make_deadata(datadea = data_year, 
                           inputs = c("Input1", "Input2", "Input3"),
                           outputs = c("Output1", "Output2"))
  
  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)))
}


# 转换为宽格式并计算均值----
efficiency_wide <- lapply(efficiency_results, function(df) {
  df_wide <- reshape(df, idvar = "DMU", timevar = "Year", direction = "wide")
  names(df_wide)[-1] <- gsub("Efficiency.", "", names(df_wide)[-1])
  df_wide$均值 <- rowMeans(df_wide[,-1], na.rm = TRUE)
  df_wide
})

# 创建 Excel 文件并添加工作表----
wb <- createWorkbook()
sheet_names <- c("CRS基础教育资源配置", "VRS基础教育资源配置", "规模效率", "SBM效率", "超效率")

for (i in seq_along(sheet_names)) {
  addWorksheet(wb, sheet_names[i])
  writeData(wb, sheet = sheet_names[i], x = efficiency_wide[[i]])
}

saveWorkbook(wb, "Table/民族八省区DEA效率值表格.xlsx", overwrite = TRUE)

# 创建 Word 文档并添加表格----
doc <- read_docx()
table_titles <- c("表 1: CRS基础教育资源配置", "表 2: VRS基础教育资源配置", "表 3: 规模效率", "表 4: SBM效率", "表 5: 超效率")

for (i in seq_along(table_titles)) {
  doc <- doc %>%
    body_add_par(value = table_titles[i], style = "centered") %>%
    body_add_flextable(value = flextable(efficiency_wide[[i]])) %>%
    body_add_break()
}

print(doc, target = "Table/民族八省区DEA效率值报告.docx")

#绘图----
## ----------------- 箱线图和折线图 ----------------- 
ggplot(efficiency_results$super, aes(x = as.factor(Year), y = Efficiency)) +
  geom_boxplot() +
  labs(title = "超效率得分", x = "Year", y = "Efficiency") +
  theme_minimal()

ggplot(efficiency_results$super, aes(x = Year, y = Efficiency, group = DMU, color = DMU)) +
  geom_line() +
  geom_point() +
  labs(title = "超效率得分", x = "Year", y = "Efficiency") +
  theme_classic()
## ----------------- 绘制热图 ----------------- 
efficiency_crs_melt <- melt(efficiency_results$crs, id.vars = c("DMU", "Year"), variable.name = "Type", value.name = "Efficiency")
ggplot(efficiency_crs_melt, aes(x = as.factor(Year), y = DMU, fill = Efficiency)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(title = "CRS 技术效率得分", x = "Year", y = "DMU") +
  theme_minimal()
## ----------------- 有效效率值统计图 ----------------- 
# 计算各 DMU 三个效率值达到有效的数量
efficiency_effective <- data.frame(
  DMU = unique(efficiency_results$crs$DMU),
  综合技术效率值 = sapply(unique(efficiency_results$crs$DMU), function(dmu) sum(efficiency_results$crs$Efficiency[efficiency_results$crs$DMU == dmu] == 1)),
  纯技术效率值 = sapply(unique(efficiency_results$crs$DMU), function(dmu) sum(efficiency_results$vrs$Efficiency[efficiency_results$vrs$DMU == dmu] == 1)),
  规模效率值 = sapply(unique(efficiency_results$crs$DMU), function(dmu) sum(efficiency_results$scale$Efficiency[efficiency_results$scale$DMU == dmu] == 1))
)

# 转换为长格式
efficiency_effective_long <- melt(efficiency_effective, id = "DMU")
efficiency_effective_long$variable <- factor(efficiency_effective_long$variable, levels = c("综合技术效率值", "纯技术效率值", "规模效率值"))

# 绘制柱状图
ggplot(efficiency_effective_long, aes(x = value, y = DMU, fill = variable)) +
  geom_bar(stat = "identity", position = "dodge", color = "grey90", size = 0.2) +
  geom_text(aes(label = value), position = position_dodge(width = 0.9), hjust = -0.3, size = 3) +
  scale_fill_manual(values = c("grey30", "grey50", "grey70"), name = "效率值类型",
                    labels = c("综合技术效率值", "纯技术效率值", "规模效率值")) +
  labs(title = " 2013-2022年各DMU效率值统计图", x = "数量", y = "DMU") +
  theme_minimal() +
  theme(legend.position = "bottom",  # 将图例放在下方
        legend.justification = "left",  # 图例居中对齐
        text = element_text(family = "noto"), legend.title = element_text(size = 10), legend.text = element_text(size = 8))


## ------------------- 交互书图表 ------------------- 
library(plotly)
efficiency_vrs_long <- melt(efficiency_results$crs, id = c("Year", "DMU"), variable.name = "Type", value.name = "Efficiency")

p <- ggplot(efficiency_vrs_long, aes(x = as.factor(Year), y = Efficiency, color = DMU, group = DMU)) +
  geom_line() +
  geom_point() +
  labs(title = "CRS 技术效率得分", x = "Year", y = "Efficiency") +
  theme_minimal()

ggplotly(p)



#表