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

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

New！ `lianxh` 命令发布了：

`. ssc install lianxh`

`. help lianxh`

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-

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

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

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

## 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
-----------------------------------------------------
``````

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

``````. 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 创造一个数字列表

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

## 5. 使用 Matplotlib 绘制三维表面图

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

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

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

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

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

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

## 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 +`

`ssc install lianxh, replace`

## 相关课程

http://lianxh.duanshu.com

### 课程一览

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

#### 关于我们

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

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

New！ `lianxh` 命令发布了：

`. ssc install lianxh`

`. help lianxh`