Logit模型拟合实战案例（Biogeme）

前言：本文详细介绍如何利用 Biogeme 拟合二项 Logit 模型，包括数据准备、分类变量的处理、Biogeme 中效用函数的设定等内容。

在【DCM 笔记】系列文章中，我们已经尝试过用 SAS 和 Python 的 statsmodels 包去拟合二项 Logit 模型；相关案例：

•Logit 模型拟合实战案例（SAS）——离散选择模型之五

•Logit 模型拟合实战案例（Python）——离散选择模型之六

相比较之下，Biogeme 在指定效用函数方面更加灵活——Biogeme 中可以非常方便地为每一个变量在不同的方案中指定不同的系数(即 Alternative Specific Coefficients)；此外，作为基于 Python 的开源软件包，Biogeme 的免费特性也吸引着很多的用户。

本篇将使用和上面两篇文章相同的数据，重点介绍如何使用 Biogeme 拟合二项 Logit 模型；并将结果和 SAS 对比。

Biogeme 相关资料：

•Biogeme 入门教程（中文版）：链接

•Biogeme 安装教程：链接

文末提供数据和代码的下载链接。

案例介绍

本例我们要研究的问题是：在申请的研究生的时候，什么样的学生更容易被录取。

原始数据保存在名为“Application.csv”的文件中（文件格式为.csv 格式），每一行代表一条申请者的记录；示例数据如下图所示：

原始数据示例

原始数据中包含 3 个自变量：

•申请者的 GRE 成绩，用变量gre表示；

•申请者的平均绩点，用变量gpa表示；

•申请者所在的本科院校的排名，用变量rank表示。

变量gre和gpa都是连续变量。rank为离散变量，只能取 1、2、3、4 中的某一个值；rank=1 对应的学校排名最高，而rank=4 对应的排名最低。

申请的结果只有两种情形：“录取”或者“拒绝”。我们用变量admit表示申请结果，显然，admit是一个二分类的变量——admit=1 表示“申请者被录取”，admit=0 表示“申请者被拒绝”。

软件准备

本例需要调用 Python 的 Biogeme 和 pandas 这两个包。运行 Python 代码之前，请确保已经正确安装相应的软件包。

数据探索

正式建模之前，可以先做一些描述性分析（Descriptive analysis）——看一看样本中各变量的均值、方差等等，以加强对数据的理解。具体实现步骤如下。

在 Python 中导入相应的包：

图 1

用 pandas 的 read_csv() 函数读取原始数据文件，并展示前 5 行：

图 2

由于 pandas 的 DataFrame 数据结构也有一个方法的名称为 rank，这容易与原始数据表中的列名rank产生混淆。将原始数据表中的列名 **rank **更改为 sch_rank。更改之后的数据表格如图 3 所示：

图 3

用 describe() 函数对样本中的各变量做描述性分析，结果如下面所示。我们可以得到每一个变量的出现的频数（count）、均值（mean）、标准差（std）、最大/小值（min/max）、百分位数（25%，50%，75%）等信息。这一步相当于 SAS 中的 Proc Means 和 Proc Freq。

图 4

当然，还可以做一下交叉频数分析，粗略地观察（离散的）自变量和因变量之间关系。例如，根据下图我们就可以看出：在样本中，当申请者所在的学校排名越高时（sch_rank=1），申请者被录取（admit=1）的比例也就越大。

图 5

BIOGEME 建模

本例中，我们的目标是建立一个二项 Logit模型，分析哪些因素会影响申请者是否被录取。

•为什么是二项 Logit 模型？因为因变量 admit 只能取 0、1 两个值；

•我们考虑 gre、gpa、sch_rank 这三个变量作为自变量——其背后的行为假设（Behavior Assumption） 是：某个申请者是否被录取，受他的成绩（由gre、gpa表示）和他所在的学校的排名（由sch_rank表示）影响。

首先，变量sch_rank是分类变量，而非连续变量；需要对其做哑变量处理。将 sch_rank=4 作为参照类（和上一篇的 SAS 案例保持一致），利用图 6 中的代码创建哑变量 sch_rank_1、sch_rank_2、sch_rank_3：

图 6

调用图 7 中的代码将一个 pandas DataFrame 对象转化成 Biogeme 的数据库对象；后面在参数拟合的时候需要用到这个数据对象：

图 7

查看 Biogeme 数据库对象的信息：样本量和变量名称：

图 8

调用图 9 中的代码，将 Biogeme 数据库中的变量转化成 Python 变量，方便后续直接调用：

图 9

Biogeme 需要指定每一个方案的效用函数。在二项 Logit 模型中，决策者 n 选择方案 i 的概率为：

Vi、Vj分别表示方案 i、方案 j 的效用函数的确定部分。

在本例中，我们将 admit=0 这一方案作为基准方案；用V1表示 admit=1 这一选项所对应的效用函数的确定部分，则：

•V1 = ASC1 + βgre gre + βgpa gpa + βsr1 sch_rank1 + βsr2 sch_rank2 + βsr3 sch_rank3

如果用V0表示 admit=0 这一选项所对用的效用函数的确定部分，那么问题来了：下面(1)、(2)两式中哪一个才是正确的的表达式呢？

•V0 = ASC0 + βgre gre + βgpa gpa + βsr1 sch_rank1 + βsr2 sch_rank2 + βsr3 sch_rank3 (1)

•V0 = 0 (2)

由于本例中的自变量都是个人属性，因此，基准方案 admit=0 中所有的自变量都应该为 0。因此，(2)式才是正确的答案。

也可换个角度理解——在所有方案的效用函数中同时增加一个相同的量并不会影响某个方案被选中的概率。

如果(1)式是正确的，此时，admit=0 和 admit=1 所对应的效用函数分别为：

•V0 = ASC0 + βgre gre + βgpa gpa + βsr1 sch_rank1 + βsr2 sch_rank2 + βsr3 sch_rank3

•V1 = ASC1 + βgre gre + βgpa gpa + βsr1 sch_rank1 + βsr2 sch_rank2 + βsr3 sch_rank3

约掉V0、V1中相同的项，上面的式子等价于：

•V0 = ASC0

•V1 = ASC1

此时，V0、V1的表达式中均不包含个人属性变量——这显然是不对的。更详细的讨论可以参见《如何将决策者的属性和方案属性同时放到MNL模型中？——离散选择模型之十八》一文。简而言之：对于个人属性，由于它们对所有的方案都是一样的，因此一般会将某一个基准方案设置为 0。

明确了、的表达式以后，我们就可以调用 Biogeme 中的 Beta class 来定义待估参数了，如图 10 所示。Beta()一共有 5 个参数：

•第 1 个参数是系数名称，一般和等号左边的变量名保持一致；

•第 2 个参数是系数的初始值，这里假设为 0；

•第 3 个参数是系数的下界，如果不知道可以用 None 替代；

•第 4 个参数是系数的下界，如果不知道可以用 None 替代；

•第 5 个参数指定是否需要估计该参数，值为 1 时相当于定义了一个常数；

图 10：定义待估参数

在 Biogeme 中定义效用函数，如图 11 所示。其中，ASC_1、B_GRE、B_GPA、B_SR1、B_SR2、B_SR2 是图 10 中定义的参数，gre、gpa、sch_rank_1、sch_rank_2、sch_rank_3 对应 Biogeme database 中相应的变量：

图 11：定义效用函数

如前所述，V0、V1分别对应admit=0、admit=1 这两个选项的效用函数；在 Biogeme 中，需要用下图中的代码将V0、V1和 admit变量的相应的值（0/1）对应起来：

图 12

调用 Biogeme 中的 models.loglogit()函数定义 Logit 模型：

图 13：定义 Logit 模型

调用 Biogeme 中的 bio.BIOGEME()函数，创建一个 biogeme 对象——这个对象既包含了之前创建的数据库对象，也包含了图 13 中定义的Logit 模型：

图 14

调用 biogeme 对象的 estimate 方法估计参数并展示结果，结果如下：

图 15

运行图 15 中的代码后，Biogeme 会生成一个 HTML 文件；里面包含了模型的总体信息以及相应的系数信息。

打开 HTML 文件。从图 16 的报告中可以看出，最终模型的对数似然函数值(LogL)为-229.2587，AIC 值为 470.5175；对比图 17，可以看出这些结果都是和 SAS 一致的（注意图 17 中给出的是-2LogL = -2* (-229.2587) = 458.5174）。

图 16：Biogeme 拟合结果摘要

图 17：SAS 拟合结果摘要

另外，对比图 18 和图 19，Biogeme 和 SAS 给出的系数都是一样的：

图 18：Biogeme 参数估计结果

图 19：SAS 参数估计结果

数据代码下载链接

百度网盘地址：

•链接：https://pan.baidu.com/s/1nQDbEROTi72rUVQVu1zjVQ

•提取码：ujj9（永久有效）

如果数据/代码不能下载，请您在后台留言告知，我会及时更新下载链接。

本篇完。

如果您觉得本篇干货满满，请您动动手指，点赞、留言、分享三连，让更多的朋友看到。谢谢！