温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
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 工具包、如何进行机器学习分析等。
中文编译稿列表如下:
在前面四篇关于 Stata 和 Python 的推文中,我们展示了如何设置 Stata 与 Python 环境、三种在 Stata 中使用 Python 的方法、如何安装 Python 软件包和如何使用 Python 软件包。在浏览这篇推文前,如果您之前不太熟悉 Python,那么阅读之前的推文也许对您有所帮助。
现在,我们将重点转移到在 Stata 中 Python 的一些实际应用上。这篇推文将演示如何使用 Stata 来估计逻辑回归模型中的边际预测,以及如何使用 Python 绘制这些预测的三维表面图。
我们将会用到 Numpy、pandas 和 Matplotlib 软件包,所以在开始之前,请您检查是否已安装它们。
数年前,我在 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 绘制预测概率的三维表面图。
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
-----------------------------------------------------
我们还需要创建数字列表,以将刻度线放置在图形的轴上。我们使用别名 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)
。
现在,我们准备开始制图。首先,将 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
生成的三维表面图展示了 age 和 weight 的值预测出的高血压可能性。概率由z轴上的表面高度和表面的颜色表示。蓝色表示较低的高血压概率,红色表示较高的高血压概率。这可能是解释回归模型中两个连续协变量之间相互作用的有效方法。
我将 Stata 命令和 Python 代码块整合到了下面的单个do文件中。并且,我添加了一些注释,以提醒您每个命令集和语句集的目的。请注意,Stata 注释以 "//" 开头,Python 注释以 "#" 开头。
// 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
Note:产生如下推文列表的命令为:
lianxh Stata Python +
安装最新版lianxh
命令:
ssc install lianxh, replace
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟,课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能,时长:2 小时,课程主页 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh