Stata-Python交互-5:边际效应三维立体图示

发布时间:2021-02-22 阅读 5265

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh

课程详情 https://gitee.com/arlionn/Course   |   lianxh.cn

课程主页 https://gitee.com/arlionn/Course

编译: 张子登 (中山大学,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-

中文编译稿列表如下:

1. 引言

在前面四篇关于 Stata 和 Python 的推文中,我们展示了如何设置 Stata 与 Python 环境、三种在 Stata 中使用 Python 的方法、如何安装 Python 软件包和如何使用 Python 软件包。在浏览这篇推文前,如果您之前不太熟悉 Python,那么阅读之前的推文也许对您有所帮助。

现在,我们将重点转移到在 Stata 中 Python 的一些实际应用上。这篇推文将演示如何使用 Stata 来估计逻辑回归模型中的边际预测,以及如何使用 Python 绘制这些预测的三维表面图。

我们将会用到 NumpypandasMatplotlib 软件包,所以在开始之前,请您检查是否已安装它们。

2. 连续型变量相互交乘的概率预测

数年前,我在 Stata News 上撰写了一篇题为 连续型变量交乘项的可视化 的文章。在文章中,我用逻辑回归模型拟合了来自美国健康与营养调查 (NHANES) 的数据。

因变量 highbp 是衡量高血压的指标。同时,我纳入了两个连续型变量 ageweight 的交乘项作为主要影响因素。

. 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.

估计出的 ageweight 交乘项的 odds ratio 为 0.9994,对应的 p 值为 0.000。解释此结果是比较困难的,因为 odds ratio 基本等于无效值 1,但 p 值却基本等于 0。这个交乘项有意义吗? 通常情况下,与仅观察 odds ratio 相比,通过查看协变量不同水平上的预测概率结果更容易得到答案。

下面的代码块使用 margins 边际方法 来估计 ageweight 所有组合中的高血压预测概率。其中,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 中的 ageweight 变量。

. 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 绘制三维表面图

现在,我们准备开始制图。首先,将 NumpypandasMatplotlib 包导入 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. 总结

生成的三维表面图展示了 ageweight 的值预测出的高血压可能性。概率由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

相关课程

连享会-直播课 上线了!
http://lianxh.duanshu.com

免费公开课:


课程一览

支持回看

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]

Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。


关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
  • 连享会-主页知乎专栏,400+ 推文,实证分析不再抓狂。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh