Logit模型拟合实战案例（Biogeme）
前⾔：本⽂详细介绍如何利⽤ Biogeme 拟合⼆项 Logit 模型，包括数据准备、分类变量的处理、Biogeme 中效⽤函数 的设定等内容。
在【DCM 笔记】系列⽂章中，我们已经尝试过⽤ SAS 和 Python 的 statsmodels 包去拟合⼆项 Logit 模型；相关案例：
·
·
相⽐较之下，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
ps切图后变模糊
当然，还可以做⼀下交叉频数分析，粗略地观察（离散的）⾃变量和因变量之间关系。例如，根据下图我们就可以看出：在样本中，当申请者所在的学校排名越⾼时（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
python入门教程(非常详细)书
如果⽤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的表达式中均不包含个⼈属性变量——这显然是不对的。更详细的讨论可以参见《》⼀⽂。简⽽⾔之：对于个⼈属性，由于它们对所有的⽅案都是⼀样的，因此⼀般会将某⼀个基准⽅案设置为 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 中定义的参
thinkphp6进销存数，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）。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。