Stata-Python交互-5：边际效应三维立体图示
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
New！
lianxh命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh
编译： 张子登 (中山大学，zhangzd5@mail2.sysu.edu.cn)
编译： 刘雅玄 (北京大学，lyxbamboo@163.com)
致谢： 本文摘自以下文章，特此感谢！
Source: Chuck Huber, The Stata Blog: Stata/Python integration part 5: Three-dimensional surface plots of marginal predictions -Link-
目录
Stata/Python 交互系列推文 源自 Stata 公司的统计项目总监 Chuck Huber 博士发表于 Stata 官网的系列博文，一共 9 篇。较为系统地介绍了 Stata 与 Python 的交互方式，包括：如何配置你的软件、如何实现 Stata 与 Python 数据集互通、如何调用 Python 工具包、如何进行机器学习分析等。
Part 1: Setting up Stata to use Python -Link- Part 2: Three ways to use Python in Stata -Link- Part 3: How to install Python packages -Link- Part 4: How to use Python packages -Link- Part 5: Three-dimensional surface plots of marginal predictions -Link- Part 6: Working with APIs and JSON data -Link- Part 7: Machine learning with support vector machines, -Link- Part 8: Using the Stata Function Interface to copy data from Stata to Python, -Link- Part 9: Using the Stata Function Interface to copy data from Python to Stata, -Link-
中文编译稿列表如下：
Stata-Python交互-9：将python数据导入Stata Stata-Python交互-8：将Stata数据导入Python Stata-Python交互-7：在Stata中实现机器学习-支持向量机 Stata-Python交互-6：调用APIs和JSON数据 Stata-Python交互-5：边际效应三维立体图示 Stata-Python交互-4：如何调用Python宏包 Stata-Python交互-3：如何安装Python宏包 Stata-Python交互-2：在Stata中调用Python的三种方式 Stata-Python交互-1：二者配合的基本设定
1. 引言
在前面四篇关于 Stata 和 Python 的推文中，我们展示了如何设置 Stata 与 Python 环境、三种在 Stata 中使用 Python 的方法、如何安装 Python 软件包和如何使用 Python 软件包。在浏览这篇推文前，如果您之前不太熟悉 Python，那么阅读之前的推文也许对您有所帮助。
现在，我们将重点转移到在 Stata 中 Python 的一些实际应用上。这篇推文将演示如何使用 Stata 来估计逻辑回归模型中的边际预测，以及如何使用 Python 绘制这些预测的三维表面图。
我们将会用到 Numpy、pandas 和 Matplotlib 软件包，所以在开始之前，请您检查是否已安装它们。
2. 连续型变量相互交乘的概率预测
数年前，我在 Stata News 上撰写了一篇题为 连续型变量交乘项的可视化 的文章。在文章中，我用逻辑回归模型拟合了来自美国健康与营养调查 (NHANES) 的数据。
因变量 highbp 是衡量高血压的指标。同时，我纳入了两个连续型变量 age 和 weight 的交乘项作为主要影响因素。
. webuse nhanes2, clear
. svy: logistic highbp c.age##c.weight
(running logistic on estimation sample)
Survey: Logistic regression
Number of strata = 31 Number of obs = 10,351
Number of PSUs = 62 Population size = 117,157,513
Design df = 31
F( 3, 29) = 418.97
Prob > F = 0.0000
--------------------------------------------------------------------------------
| Linearized
highbp | Odds Ratio Std. Err. t P>|t| [95% Conf. Interval]
---------------+----------------------------------------------------------------
age | 1.100678 .0088786 11.89 0.000 1.082718 1.118935
weight | 1.07534 .0063892 12.23 0.000 1.062388 1.08845
|
c.age#c.weight | .9993975 .0001138 -5.29 0.000 .9991655 .9996296
|
_cons | .0002925 .0001194 -19.94 0.000 .0001273 .0006724
--------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
估计出的 age 和 weight 交乘项的 odds ratio 为 0.9994，对应的 p 值为 0.000。解释此结果是比较困难的，因为 odds ratio 基本等于无效值 1，但 p 值却基本等于 0。这个交乘项有意义吗？ 通常情况下，与仅观察 odds ratio 相比，通过查看协变量不同水平上的预测概率结果更容易得到答案。
下面的代码块使用 margins 边际方法 来估计 age 和 weight 所有组合中的高血压预测概率。其中，age 从20到80岁(以5为增量)取值，weight 从40到180公斤(以5为增量)取值。
代码选项 save(predictions，replace) 将预测结果保存到名为 predictions.dta 的数据集。之后，打开predictions.dta 数据集，在数据集中重命名三个变量后再次保存。
. webuse nhanes2, clear
. svy: logistic highbp age weight c.age#c.weight
. quietly margins, at(age=(20(5)80) weight=(40(5)180)) ///
. vce(unconditional) saving(predictions, replace)
. use predictions, clear
. rename _at1 age
. rename _at2 weight
. rename _margin pr_highbp
. save predictions, replace
在我 之前的 Stata News 文章 中，我使用 双变量轮廓图 绘制了高血压概率预测的轮廓图。 在本文中，我将使用 Python 绘制预测概率的三维表面图。
3. 使用 pandas 将边际预测结果读入 Python
Python 必须能够读取储存在 predictions.dta 中的数据，才能进行后续三维表面图的绘制。我们首先使用别名 pd 将 pandas 包导入 Python；进而使用 pandas 包中的 read_stata() 函数将 predictions.dta 读取到名为 data 的 pandas 数据框中：
. python:
------------------ python (type end to exit) --------
>>> import pandas as pd
>>> data = pd.read_stata("predictions.dta")
>>> data
pr_highbp age weight
0 0.020091 20 40
1 0.027450 20 45
2 0.037401 20 50
3 0.050771 20 55
4 0.068580 20 60
.. ... ... ...
372 0.954326 80 160
373 0.958618 80 165
374 0.962523 80 170
375 0.966072 80 175
376 0.969296 80 180
[377 rows x 3 columns]
>>> end
-----------------------------------------------------
然后，我们可以通过键入 dataframe[‘varname’] 在数据框中引用变量。例如，我们可以通过键入 data ['age'] 在数据框 data 中引用变量 age。
. python:
------------------ python (type end to exit) --------
>>> import pandas as pd
>>> data = pd.read_stata("predictions.dta")
>>> data['age']
0 20
1 20
2 20
3 20
4 20
..
372 80
373 80
374 80
375 80
376 80
Name: age, Length: 377, dtype: int8
>>> end
-----------------------------------------------------
也可以用 dataframe [['varname1'，'varname2']] 来引用数据框中的多个变量。例如，我们可以通过键入 data [['age'，'weight']] 来引用数据框 data 中的 age 和 weight 变量。
. python:
------------------ python (type end to exit) --------
>>> import pandas as pd
>>> data = pd.read_stata("predictions.dta")
>>> data[['age', 'weight']]
age weight
0 20 40
1 20 45
2 20 50
3 20 55
4 20 60
.. ... ...
372 80 160
373 80 165
374 80 170
375 80 175
376 80 180
[377 rows x 2 columns]
>>> end
-----------------------------------------------------
4. 使用 NumPy 创造一个数字列表
我们还需要创建数字列表，以将刻度线放置在图形的轴上。我们使用别名 np 将 NumPy 包导入 Python。然后，我们可以使用 NumPy 包中的 arange() 函数创建列表。下面的示例创建了一个名为 mylist 的数字列表，该列表取值范围是 20 至 90 ，以 10 为增量。
. python:
------------------ python (type end to exit) --------
>>> import numpy as np
>>> mylist = np.arange(20,90, step=10)
>>> mylist
array([20, 30, 40, 50, 60, 70, 80])
>>> end
-----------------------------------------------------
您可能会感到惊讶：为什么列表中不包含数字 90 ？这并不是程序或人为的错误，而是 NumPy 中 arange() 函数的特性。如果想要在列表中包括90，则可以键入 np.arange(20,100，step = 10)。
5. 使用 Matplotlib 绘制三维表面图
现在，我们准备开始制图。首先，将 Numpy、pandas 和 Matplotlib 包导入 Python。此处，我们使用别名 plt 从 Matplotlib 包中导入 pyplot 模块。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.plot_trisurf(data['age'], data['weight'], data['pr_highbp'])
cmap=plt.cm.Spectral_r)
end
我在上面的代码块中添加了语句 matplotlib.use('TkAgg')。该语句能确保 Matplotlib 在 Windows 10 环境中的 Python 正常运行。Matplotlib 针对不同的目的和平台应使用不同的渲染引擎，故您可能需要在您的计算环境中使用其他渲染引擎。您可以在 Stata FAQ 的文章 如何用 Stata 使用 Matplotlib 中了解更多相关信息。
接下来，我们用 data = pd.read_stata("predictions.dta") 语句，调用 pandas 包中的 read_stata() 函数将 predictions.dta 读取到名为 data 的 pandas 数据框中。
然后，我们使用 ax = plt.axes(projection='3d') 语句，调取 pyplot 模块中的 axes() 函数来定义名为 ax 的三维轴集。
接下来，我们使用 ax.plot_trisurf(data['age'], data['weight'], data['pr_highbp']) 语句，调用 pyplot 模块中的 plot_trisurf() 函数来渲染三维表面图。
默认情况下，曲面图以纯蓝色渲染。我们使用 cmap = plt.cm.Spectral_r 选项将阴影添加到绘图中 (对应语句为： cmap=plt.cm.Spectral_r) )。配色方案 Spectral_r 将以蓝色显示较低的高血压概率，以红色显示较高的概率。
虽然默认的轴刻度看起来很合理，但我们可能希望对其进行自定义处理。y 轴看起来有些混乱，我们可以使用 axes 模块中的 set_yticks() 函数修改刻度之间的增量。NumPy 模块中的 arange() 函数定义了一个从40到200的、步长为40的数字列表。我们可以使用类似的语句向 x 和 z 轴添加自定义刻度线。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.plot_trisurf(data['age'], data['weight'], data['pr_highbp'],
cmap=plt.cm.Spectral_r)
# New added
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
end
接下来，我们可以使用 axes 模块中的 set_title() 函数为图添加标题。同样，我们可以分别使用 set_xlabel() 、 set_ylabel() 和 set_zlabel() 函数向 x，y 和 z 轴添加标签。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
# New added
ax.set_title("Probability of Hypertension by Age and Weight")
ax.set_xlabel("Age (years)")
ax.set_zlabel("Probability of Hypertension")
end
您可以使用 view_init() 函数调整视角。elev 选项可调整高程，而 azim 选项可调整方位角。两者均以度为单位。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
ax.set_title("Probability of Hypertension by Age and Weight")
ax.set_xlabel("Age (years)")
ax.set_ylabel("Weight (kg)")
ax.set_zlabel("Probability of Hypertension")
# New added
ax.view_init(elev=30, azim=240)
end
如果我希望从下至上查看 z 轴标题而不是从上至下，我可以使用 set_rotate_label（False） 函数和 set_xlabel（） 函数中 rotation = 90 选项的组合来更改此设置。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
ax.set_title("Probability of Hypertension by Age and Weight")
ax.set_xlabel("Age (years)")
ax.set_ylabel("Weight (kg)")
# New added
ax.zaxis.set_rotate_label(False)
ax.set_zlabel("Probability of Hypertension" `, rotation=90)
ax.view_init(elev=30, azim=240)
end
最后，我们可以使用 savefig() 函数将图形另存为每英寸1200点的便携式网络图形 (.png) 文件。
python:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
data = pd.read_stata("predictions.dta")
ax = plt.axes(projection='3d')
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
ax.set_title("Probability of Hypertension by Age and Weight")
ax.set_xlabel("Age (years)")
ax.set_ylabel("Weight (kg)")
ax.zaxis.set_rotate_label(False)
ax.set_zlabel("Probability of Hypertension", rotation=90)
ax.view_init(elev=30, azim=240)
# New added
plt.savefig("Margins3d.png", dpi=1200)
end
6. 总结
生成的三维表面图展示了 age 和 weight 的值预测出的高血压可能性。概率由z轴上的表面高度和表面的颜色表示。蓝色表示较低的高血压概率，红色表示较高的高血压概率。这可能是解释回归模型中两个连续协变量之间相互作用的有效方法。
我将 Stata 命令和 Python 代码块整合到了下面的单个do文件中。并且，我添加了一些注释，以提醒您每个命令集和语句集的目的。请注意，Stata 注释以 "//" 开头，Python 注释以 "#" 开头。
7. 代码示例(example.do)
// Fit the model and estimate the marginal predicted probabilities with Stata
webuse nhanes2, clear
logistic highbp c.age##c.weight
quietly margins, at(age=(20(5)80) weight=(40(5)180)) ///
saving(predictions, replace)
use predictions, clear
rename _at1 age
rename _at2 weight
rename _margin pr_highbp
save predictions, replace
// Create the three-dimensional surface plot with Python
python:
# Import the necessary Python packages
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('TkAgg')
# Read (import) the Stata dataset “predictions.dta”
# into a pandas data frame named “data”
data = pd.read_stata(“predictions.dta”)
# Define a 3-D graph named “ax”
ax = plt.axes(projection=’3d’)
# Render the graph
ax.plot_trisurf(data[‘age’], data[‘weight’], data[‘pr_highbp’],
cmap=plt.cm.Spectral_r)
# Specify the axis ticks
ax.set_xticks(np.arange(20, 90, step=10))
ax.set_yticks(np.arange(40, 200, step=40))
ax.set_zticks(np.arange( 0, 1.2, step=0.2))
# Specify the title and axis titles
ax.set_title(“Probability of Hypertension by Age and Weight”)
ax.set_xlabel(“Age (years)”)
ax.set_ylabel(“Weight (kg)”)
ax.zaxis.set_rotate_label(False)
ax.set_zlabel(“Probability of Hypertension”, rotation=90)
# Specify the view angle of the graph
ax.view_init(elev=30, azim=240)
# Save the graph
plt.savefig("Margins3d.png", dpi=1200)
end
8. 相关推文
Note：产生如下推文列表的命令为：
lianxh Stata Python +
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata入门 专题：Stata程序 专题：文本分析-爬虫 专题：Python-R-Matlab 专题：其它
相关课程
连享会-直播课 上线了！
http://lianxh.duanshu.com
免费公开课：
直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页 部分直播课 课程资料下载 (PPT，dofiles等)
课程一览
支持回看
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。 连享会-主页 和 知乎专栏，400+ 推文，实证分析不再抓狂。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD
New！
lianxh命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh
热门资讯
