R实现Malmquist指数及分解

R
科研
发布日期:   2024-07-02
文章字数:   2k
阅读时长:   10 分
R实现Malmquist指数计算与分析流程 问题定义 计算和分析Malmquist生产率指数 处理过程 输入: Excel 数据 (Data/数据模板.xlsx) 初始化 (加载R包) 核心计算: Malmquist指数 (使用 malm 函数, VRS, Input-Oriented) 结果提取与格式化 (effch, tech, pech, sech, tfpch) 聚合计算 (按年均/DMU均/总体均) 输出结果 Excel 报告 (年度/DMU均值) Word 报告 (格式化表格) 可视化图表 (ggplot2 点图/线图) 控制台输出 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
# ---- 初始化 ----

library(productivity) library(openxlsx) library(dplyr) library(readxl)
library(officer) library(flextable) library(ggrepel)

# ---- 数据导入

datafile \<- read_excel("Data/数据模板.xlsx") data \<-
as.data.frame(datafile)

# ---- 计算Malmquist指数 ----

Malmquist \<- malm(data, id.var = "DMU", time.var = "Year", x.vars =
c("Input1", "Input2", "Input3"), #投入指标 y.vars = c("Output1",
"解释变量2"), rts = "vrs", orientation = "in")#产出指标 BCC 模型

## ---- 提取并四舍五入效率值 ----

results \<- data.frame( DMU = Malmquist$Changes$DMU, Year.0 =
Malmquist$Changes$Year.0, Year.1 = Malmquist$Changes$Year.1, effch =
round(Malmquist$Changes$effch, 3), tech = round(Malmquist$Changes$tech,
3), pure_inp_effch = round(Malmquist$Changes$pure.inp.effch, 3),
inp_scalech = round(Malmquist$Changes$inp.scalech, 3), malmquist =
round(Malmquist$Changes$malmquist, 3) )

## ---- 计算相邻年份的平均值 ----

averages \<- results %\>% group_by(Year.1, Year.0) %\>% summarise(
effch_avg = round(mean(effch, na.rm = TRUE), 3), tech_avg =
round(mean(tech, na.rm = TRUE), 3), pure_inp_effch_avg =
round(mean(pure_inp_effch, na.rm = TRUE), 3), inp_scalech_avg =
round(mean(inp_scalech, na.rm = TRUE), 3), malmquist_avg =
round(mean(malmquist, na.rm = TRUE), 3) ) %\>% mutate(Years =
paste(pmin(Year.0, Year.1), pmax(Year.0, Year.1), sep = "-"))

## ---- 准备结果数据框 ----

results_avg \<- data.frame( Years = averages$Years,
  effch = averages$effch_avg, tech = averages$tech_avg,
  pure_inp_effch = averages$pure_inp_effch_avg, inp_scalech =
averages$inp_scalech_avg,
  malmquist = averages$malmquist_avg )

## ---- 计算总体平均值 ----

overall_avg \<- results %\>% summarise( Years = "平均值", effch =
round(mean(effch, na.rm = TRUE), 3), tech = round(mean(tech, na.rm =
TRUE), 3), pure_inp_effch = round(mean(pure_inp_effch, na.rm = TRUE),
3), inp_scalech = round(mean(inp_scalech, na.rm = TRUE), 3), malmquist =
round(mean(malmquist, na.rm = TRUE), 3) )

## ---- 将结果与总体平均值合并 ----

results_avg \<- bind_rows(results_avg, overall_avg)

# ---- 计算DMU平均值 ----

dmu_avg \<- results %\>% group_by(DMU) %\>% summarise( effch_avg =
round(mean(effch, na.rm = TRUE), 3), tech_avg = round(mean(tech, na.rm =
TRUE), 3), pure_inp_effch_avg = round(mean(pure_inp_effch, na.rm =
TRUE), 3), inp_scalech_avg = round(mean(inp_scalech, na.rm = TRUE), 3),
malmquist_avg = round(mean(malmquist, na.rm = TRUE), 3) ) %\>% ungroup()

# 添加“平均值”行

dmu_avg \<- bind_rows( dmu_avg, dmu_avg %\>% summarise( DMU = "平均值",
effch_avg = round(mean(effch_avg, na.rm = TRUE), 3), tech_avg =
round(mean(tech_avg, na.rm = TRUE), 3), pure_inp_effch_avg =
round(mean(pure_inp_effch_avg, na.rm = TRUE), 3), inp_scalech_avg =
round(mean(inp_scalech_avg, na.rm = TRUE), 3), malmquist_avg =
round(mean(malmquist_avg, na.rm = TRUE), 3) ) )

## ---- 准备 DMU 结果数据框 ----

results_DMU \<- dmu_avg

## ---- 打印 DMU 结果数据框 ----

print(results_DMU)

# ---- 将结果写入Excel ----

wb \<- createWorkbook() addWorksheet(wb, "平均效率指数") writeData(wb,
sheet = "平均效率指数", x = results_avg)

# 添加一个新的工作表用于DMU平均值

addWorksheet(wb, "DMU平均指数") writeData(wb, sheet = "DMU平均指数", x =
dmu_avg)

saveWorkbook(wb, file = "Table/Average_Indices_Across_Years.xlsx",
overwrite = TRUE)

# ---- 将结果写入Word ----

#Word创建Flextable ft1 \<- flextable(results_avg) %\>% theme_booktabs()
%\>% autofit()

ft2 \<- flextable(dmu_avg) %\>% theme_booktabs() %\>% autofit()

## ---- 准备中文Flextable ----

results_avg_cn \<- results_avg %\>% mutate( Years = ifelse(Years ==
"Overall Average", "平均值", Years) )

colnames(results_avg_cn) \<- c("年份", "技术效率变化指数",
"技术进步效率指数", "纯技术效率指数", "规模效率指数", "全要素生产率")

dmu_avg_cn \<- dmu_avg %\>% mutate( DMU = ifelse(DMU == "平均值",
"平均值", DMU) )

colnames(dmu_avg_cn) \<- c("地区", "技术效率变化指数",
"技术进步效率指数", "纯技术效率指数", "规模效率指数", "全要素生产率")

ft3 \<- flextable(results_avg_cn) %\>% theme_booktabs() %\>% autofit()

ft4 \<- flextable(dmu_avg_cn) %\>% theme_booktabs() %\>% autofit()

## ---- 创建并保存Word文档 ----

doc \<- read_docx() %\>% body_add_par("表 （面板数据期间）我国
xx效率分年度变动指数及其分解指数", style = "centered") %\>%
body_add_flextable(ft1) %\>% body_add_par(" ", style = "Normal") %\>%
body_add_par("表 （面板数据期间）我国 xx效率分年度变动指数及其分解指数",
style = "centered") %\>% body_add_flextable(ft3) %\>% body_add_par(" ",
style = "Normal") %\>% body_add_par("表
（面板数据期间）各地区效率整体变动指数及其分解指数", style = "centered")
%\>% body_add_flextable(ft2) %\>% body_add_par(" ", style = "Normal")
%\>% body_add_par("表
（面板数据期间）各地区效率整体变动指数及其分解指数", style = "centered")
%\>% body_add_flextable(ft4)

print(doc, target =
"Table/Malmquist指数及其分解指数结果（三线表）.docx")

# ---- 绘制图表 ----

# 过滤掉“总体平均”数据点

results_avg_filtered \<- results_avg[results_avg\$Years != "平均值", ]
\## 绘制图表1---- library(ggrepel)

# 绘制图表

p \<- ggplot(data = results_avg_filtered, aes(x = factor(Years), group =
1)) + geom_point(aes(y = effch, shape = "技术效率变化指数"), size = 3,
color = "black") + geom_point(aes(y = tech, shape = "技术进步效率指数"),
size = 3, color = "black", alpha = 0.2) + geom_point(aes(y =
pure_inp_effch, shape = "纯技术效率指数"), size = 3, color = "black",
alpha = 0.7) + geom_point(aes(y = inp_scalech, shape = "规模效率指数"),
size = 3, color = "black", alpha = 0.5) + geom_point(aes(y = malmquist,
shape = "全要素生产率指数"), size = 3, color = "black") +
geom_text_repel(aes(y = effch, label = round(effch, 3)), color =
"black") + geom_text_repel(aes(y = tech, label = round(tech, 3)), color
= "black") + geom_text_repel(aes(y = pure_inp_effch, label =
round(pure_inp_effch, 3)), color = "black") + geom_text_repel(aes(y =
inp_scalech, label = round(inp_scalech, 3)), color = "black") +
geom_text_repel(aes(y = malmquist, label = round(malmquist, 3)), color =
"black") + scale_shape_manual(name = " ", values = c(16, 17, 0, 3, 4)) +
labs( title = "2013-2022 年民族八省区Malmquist指数均数情况", \# 添加标题
x = "年份", y = "Malmquist指数" ) + theme_minimal() + theme( axis.text.x
= element_text(angle = 45, hjust = 1), legend.position = "bottom",
plot.margin = margin(30, 30, 50, 30), \# 增加顶部边距 legend.box.margin
= margin(10, 10, 10, 10), \# 增加图例区域的边距 plot.title =
element_text(hjust = 0.5, vjust = 1) \# 标题居中 ) + guides( shape =
guide_legend(nrow = 2, byrow = TRUE, override.aes = list(size = 4, color
= "black")) \# 增加图例形状的大小和颜色 )

# 绘制图表

print(p)

## 绘制图表 p2----

# 过滤掉“总体平均”数据点

results_avg_filtered \<- results_avg[results_avg\$Years != "平均值", ]
library(ggrepel) p2 \<- ggplot(data = results_avg_filtered, aes(x =
factor(Years), group = 1)) + geom_line(aes(y = effch, linetype =
"技术效率变化指数"), size = 1, color = "black") + geom_line(aes(y =
tech, linetype = "技术进步效率指数"), size = 1, color = "black", alpha =
0.2) + geom_line(aes(y = pure_inp_effch, linetype = "纯技术效率指数"),
size = 1, color = "black", alpha = 0.7) + geom_line(aes(y = inp_scalech,
linetype = "规模效率指数"), size = 1, color = "black", alpha = 0.5) +
geom_line(aes(y = malmquist, linetype = "全要素生产率指数"), size = 1,
color = "black") + geom_text_repel(aes(y = effch, label = round(effch,
3)), color = "black") + geom_text_repel(aes(y = tech, label =
round(tech, 3)), color = "black") + geom_text_repel(aes(y =
pure_inp_effch, label = round(pure_inp_effch, 3)), color = "black") +
geom_text_repel(aes(y = inp_scalech, label = round(inp_scalech, 3)),
color = "black") + geom_text_repel(aes(y = malmquist, label =
round(malmquist, 3)), color = "black") + scale_linetype_manual(name = "
", values = c("solid", "twodash", "dashed", "dotted", "dotdash"), labels
= c("技术效率变化指数", "技术进步效率指数", "纯技术效率指数",
"规模效率指数", "全要素生产率指数")) + labs( title =
"效率变动及其分解指数变化趋势图",\
x = "年份", y = "效率值" ) + theme_bw() + theme( legend.position =
"bottom", legend.title = element_blank(), \# 去掉图例标题 legend.text =
element_text(size = 9), \# 设置图例文本大小 legend.key.width = unit(2,
"cm"), \# 调整图例键的宽度 legend.key.height = unit(0.5, "cm"), \#
调整图例键的高度 axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5, vjust = 1) \# 标题居中 ) +
guides(linetype = guide_legend(nrow = 2, byrow = TRUE))

# 绘制图表

print(p2)

# ---- 打印分年度变动指数 ----

print(results_avg) \# ---- 打印整体变动指数 ---- print(results_DMU)
文章作者: 徐老师
文章链接: https://xulaoshi.com/r-shi-xian-malmquist-zhi-shu-ji-fen-jie/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来源 徐老师 !
R
赏
  目录