二项回归模型： 你需要知道的一切

使用Python和statsmodels的分布实现 

本文分为两个部分。

1. 第一部分：二项式回归模型简介，我们介绍二项式回归模型，理解它是如何融入广义线性模型家族的，以及为什么它可以用来预测随机事件的概率。

2. 第二部分：使用二项式回归模型：我们使用 Python 和 statsmodels 库在大家所熟知的泰坦尼克号数据集上训练二项式回归模型。我们将解释为什么二项式回归模型是预测泰坦尼克号生存概率的正确模型。我们还将学习如何解释拟合模型的回归系数，这是一项必须学习的技能，在泰坦尼克号数据集的情况下，它产生了惊艳的结果。

第一部分 

简介：二项式回归模型

二项式回归模型可用于在给定回归变量的矢量的情况下，预测看到一个事件的几率（odds）。例如，考虑到当前的温度、湿度、气压、年份、地理位置、海拔高度等。数据科学家可以使用二项回归模型来预测未来2小时内下雨的几率。

在二项式回归模型中，因变量   是一个离散的随机变量，其数值为0、1、5、6、7等。每个值代表了在   次试验中观察到的 "成功"的次数。因此，  服从二项式分布。

下面的方程式给出了在   独立的伯努利试验中观察到   次成功的概率。每个伯努利试验的成功概率=p，失败概率=(1- p)。抛硬币是最简单的伯努利试验的例子，其中

二元分布的随机变量   的 ：

垂直括号内的 是组合的符号。它给你提供了从一组   可能的结果中选择   结果的不同方式的数量。

在回归模型中，我们将假设因变量   取决于一个   大小的回归变量矩阵  。  中的第   行可以表示为   ，这是一个大小为   的向量，它对应于第   个结果   。

人们通常把   取某个值的概率   表示为 以回归变量   取值   为条件。在符号形式上，它可以写成  ，可以读作   是   的概率，当以   等于   为条件。

现在我们可以在   对   进行回归的情况下，说明二项分布的   的概率分布如下：

• 在上述   方程的上，我们将用条件概率   取代无条件概率 。
• 我们将用条件概率    取代观察到成功的无条件概率  ，其中   是数据集第   行观察到成功的概率，即当回归向量   时。

通过这两个替换，二元分布   的PMF变成如下：

对   进行回归的二项分布   的条件概率分布

在上式中，对于某些 ，观察到成功的概率   ，通常表示为   的函数g(.)。用符号表示。

成功的概率表示为回归矩阵   的第   行的函数。

上述一组概念可以用一个插图来可视化，如下所示：

每组的 y 值的分布符合均值为   且   ,m 与 p 就决定了二项式的分布，回顾一下二项式的分布公式。

在上图中，  是10个二项分布的随机变量。每组的 y 值的分布符合均值为   且   ,m 与 p 就决定了二项式的分布，回顾一下二项式的分布公式。 它们也恰好是由因变量   的组成，它是一个（10 X 1）矩阵，如下所示：

y=  
[[y_1],  
 [y_2],  
 [y_3],  
  …,  
 [y_10]] 

在这种情况下，因为只涉及一个回归变量，所以，相应的回归变量矩阵   也正好是一个(10 X 1)的矩阵。

X =  
[[1],  
 [2],  
 [3],  
 [4],  
 …,  
 [10]]

因为   是一个随机变量，其分布范围为   ，该图显示了对于   的每个值，   可以取在其期望值   周围的任何二项分布值，其中   和   正如我们前面看到的是   的某个函数g(.)。因此，  的期望值是  ，且可以表示为   的某个函数。

二项式回归模型是广义线性模型系列的一部分。广义线性模型被用来模拟响应变量   的期望值与解释变量向量   的线性组合之间的关系。

 和   之间的关系是通过一个合适的 链接函数 来表达的，如下所示：

在上述方程中，g(.)是连接   在   上的[[条件概率]]与回归变量   的线性组合的链接函数。   是大小为   的回归变量矩阵，其中   为行数，   为每行的回归变量，   =   是这个大小为 的矩阵中的第   行，  是回归系数的   向量。

当   是二项分布时，我们感兴趣的是固定的一次伯努利试验的概率   与   的特定值之间的关系，即  ，或者简而言之， 。所以二项回归模型的GLM方程式可以写成如下：

在二项回归模型中，链接函数 g(.) 有以下四种形式（为了简单起见，我们不再提及每种形式中的条件符号  ，但它存在）。

Logistic(logit)链接函数，也称为对数-赔率函数。

Logistic被称为对数-赔率函数，它表示为成功概率与失败概率之比，即成功概率的对数。我们将在本文的后面使用这个链接函数。

probit链接函数:

probit（概率单位的简称）链接函数被用来模拟一个具有二元是/否结果的事件的发生。这个链接函数表示标准正态分布N(0,1)的累积分布函数Φ(.)的逆。

Colin Cameron和Pravin K. Trivedi的计数数据的回归分析一书在3.6节中对Probit链接函数做了很好的介绍。有序和其他离散选择模型。在书里，你还会发现一个非常清晰的推导，即为什么Probit模型的链接函数恰好是正态分布的CDFΦ(.)的逆。

对数-对数函数：

对数函数对于建模"类似泊松的计数过程 "很有用，其中概率分布的参数（通常包含平均值）位于概率分布公式的指数中，而且参数也被表达为回归变量的线性组合的指数。因此，它具有双指数格式： ，因此需要连续进行两次对数运算，以使   项下降到'地面上'。

互补的对数-对数联系函数:

互补对数-对数之所以被称为对数，是因为它的作用的对象是  ，即失败的概率，而不是 。

在这篇文章中，我们将使用logistic又称logit又称log-odds链接函数来建立我们的二项式回归模型。在这里，它再次出现，这次是以一种略微不同的方式表达。在R.H.S上，我用加粗的向量符号代替了求和。

对数-赔率链接函数

第二部分 

使用二项式回归模型来预测泰坦尼克号的生存几率

我们将以泰坦尼克号的数据集为例，了解适合二项式回归模型的各种使用情况。

这里有一个原始数据集的链接在文末。

泰坦尼克号数据集包含了命运多舛的泰坦尼克号上2229名乘客中的887人的信息。每个乘客的记录包含以下属性。

• 姓名 Name
• 乘客等级 Passenger class (1/2/3)
• 性别 Sex
• 年龄 Age
• 乘客是否有兄弟姐妹、父母或子女陪同
• 支付的费用
• 是否获救 survived (1=Survived, 0=Died)

使用Python和 Pandas 数据分析库，让我们把数据集加载到一个Pandas数据框中，并打印出前几行的数据。

import pandas as pd
df = pd.read_csv('titanic_dataset.csv', header=0)
df.head(10)

打印输出:

我们将把注意力集中在四个关键属性上。

• Pclass (乘客等级)
• Age
• Sex
• Survived (Whether they survived)

让我们从数据框中删除其余的列:

#删除不必要列
df = df.drop(['Name','Siblings/Spouses Aboard', 'Parents/Children Aboard', 'Fare'], axis=1)
#打印10个结果
df.head(10)

这就打印出了以下内容:

为二项式回归模型建立模型

我们将推测，当泰坦尼克号沉没时， [Pclass, Age, Sex]的组合极大地影响了乘客的生存几率。

请注意，Survived 列包含一个[0, 1] 伯努利随机变量。

那么我们是否可以这样说呢？

回归变量    = [Pclass, Age, Sex],  而且。

因变量是布尔变量    = [Survived]。

不，我们不能。让我们来看看为什么...

为什么不使用Probit回归模型？

由于   是一个布尔变量，所以使用Probit回归模型似乎是一个直接的demo。但是请注意，如果一个人不幸在泰坦尼克号这样的船上，他想知道的不是二元问题的答案：我是100%确定地活下来还是100%确定地死掉？相反，更有必要知道的是生存的概率。

例如，如果你是一个22岁的女人，在船的二等舱，你想知道你的生存几率是10分之1，4分之1，50分之1等等。

但是按照泰坦尼克号数据集的组织方式，响应变量 survived 有一个yes/no即1/0的格式。换句话说，survived 有一个伯努利分布，即 。

  ,

其中   是介于0和1之间的概率。

我们想要的是   表示几率，即在   个独立的、相同的试验中，成功（存活）与失败（死亡）的比率。

换句话说，我们希望   有一个 Log-Odds分布。

因此，我们要做的不是使用真或假、1或0类型的Probit回归模型，而是建立一个二项式回归模型，其中响应变量是二项式分布，链接函数是Logit即对数函数。

链接函数将允许我们把生存几率与回归变量   =[Pclass, Age & Sex] 的线性组合加上截距联系起来，如下所示:

如何将   从伯努利变量转化为二项式变量。

为了将响应变量   从伯努利变量转化为二项式变量，我们必须将共享相同的   值组合[Pclass, Age , Sex] 的数据集行分组。在我们做这件事之前，有一件小事我们需要处理，那就是对年龄属性进行分组。你看年龄在数据集中的表达方式，是一个连续变量，范围从0.42到80。

分别为0.42岁和0.67岁的婴儿会有不同的存活几率，这似乎很难说得通。对于年龄为26岁、27岁、28岁、29岁等的青年来说，同样的逻辑也成立，其他情况也是如此。

我们需要使年龄数据更加细化，以限制群体的数量。让我们通过把整个年龄范围分成大小为5岁的桶，并给每个桶贴上这样的标签来做到这一点。

(0, 5] → 5
(5, 10] → 10
(10, 15] → 15

pandas.cut()方法可以非常整齐地进行分桶。

#定义桶的bins
age_range_bins=[0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80]
#定义 标签 
Num labels = Num bins - 1age_range_labels=[5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80]
#切割 
columndf['Age_Range']=pd.cut(df['Age'],age_range_bins,labels=age_range_labels)
#Print
df.head(10)

我们得到以下输出。注意我们添加的新Age_Range列:

让我们删除 age 列，因为我们将使用Age_Range来代替。

#删除 age
df = df.drop(['Age'],axis=1)

现在我们已经按照我们想要的方式设置了解释变量，让我们开始工作，按照 [Pclass, Sex, Age_Range] 的组合将样本分组。

#Group by ['Pclass', 'Sex', 'Age_Range'] 
groups = df.groupby(['Pclass', 'Sex', 'Age_Range'])
#Get the counts for each group. This is the number of passengers in each group who have survived
df_grouped_counts_1 = groups.count().reset_index()
#Get the size (number of passengers) in each group
df_grouped_survived_1 = groups.sum().reset_index()

将每组的幸存者人数和乘客人数合并到每个分组的数据框中。(一会儿我们会看到这对我们有什么帮助）。

最后，让我们构建一个新的数据框架，其中包含。

• 分组的列Pclass, Sex, Age_Range。
• 每组中相应的幸存者人数。
• 每组中的乘客总数，即组的大小，以及。
• 每组中死亡的乘客人数。

#Create a new Data Frame
df_grouped = pd.DataFrame()
#Copy over the Pclass, Sex and Age Range columns
df_grouped['Pclass'] = df_grouped_counts_1['Pclass']
df_grouped['Sex'] = df_grouped_counts_1['Sex']
df_grouped['Age_Range'] = df_grouped_counts_1['Age_Range']
#Copy over the num passengers from the counts grouped Data Frame
df_grouped['Total'] = df_grouped_counts_1['Survived']
#Copy over the num survivors from the summation grouped Data Frame
df_grouped['Survived'] = df_grouped_survived_1['Survived']
#Add a column containing the number who died
df_grouped['Died'] = df_grouped['Total'] - df_grouped['Survived']

让我们打印出分组数据集的前20行。

df_grouped.head(20)

泰坦尼克号数据集的前20行按[pClass, Sex, Age_Range] 让我们看看分组后的数据集告诉我们什么。我突出强调了第9、14和19行，以示说明:

• 在第9行中，我们发现在年龄范围（45，50）内有10名妇女持有头等舱机票，其中9人幸存。因此，这个群体中的女性生存的几率是相当大的（9比1），特别是如果她占据了头等舱。

• 在第19行，我们看到有4名年龄在15-20岁的男性乘客，其中只有一人幸存下来。

• 在第14行，我们看到没有任何年龄在70-75岁的女性乘客占据头等舱。这就是为什么我们在该组的总列中看到了[1, female, 75]。

让我们从数据框中删除所有这样的NaN行。

df_grouped = df_grouped.query("Total !=0")

使用Python和statsmodels建立二项式回归模型

在我们建立二项式模型之前，让我们处理最后一项数据准备工作，即用整数1和2替换 "女性 "和 "男性 "的字符串。

df_grouped=df_grouped.replace(to_replace={'female': 1, 'male': 2})

使用Python和statsmodels建立二项式回归模型 

我们将使用statsmodels库提供的优秀支持来建立和训练二项回归模型。

让我们分出训练和测试数据集。

import numpy as np
#切割训练数据与测试数据
mask = np.random.rand(len(df_grouped)) < 0.85
df_train = df_grouped[mask]
df_test = df_grouped[~mask]

让我们使用patsy的语法来设置回归模型的公式。在下面提到的公式中，我们所说的因变量是一个由数据框架中的 Survived 和 Died 列组成的矩阵，而回归变量是 Pclass 、Age_Range 和 Sex 。

#构建二项式回归公式
formula = 'Survived + Died ~ Pclass + Age_Range + Sex'

使用这个公式，让我们从一分钟前创建的训练和测试数据框中划出   和   设计矩阵。

from patsy import dmatrices
#使用回归公式从训练数据中提取训练矩阵
y_train, X_train = dmatrices(formula, df_train,return_type='dataframe')
#使用回归公式从测试数据中提取训练矩阵
y_test, X_test = dmatrices(formula, df_test, return_type='dataframe')

接下来，我们将 X_train 和 y_train 送入二项回归模型类的一个实例，并训练模型。

import statsmodels.api as sm
binom_model = sm.GLM(y_train, X_train, family=sm.families.Binomial())
binom_model_results = binom_model.fit()

让我们打印出拟合模型的摘要:

print(binom_model_results.summary())

这就打印出了以下内容:

如何解释回归结果

在上面的输出中，statsmodels告诉我们，它已经训练了一个二项式的广义线性模型，因为，我们要求它这样做，它使用了对数链接函数，并且它使用了迭代重加权最小二乘法（IRLS）算法来训练我们的模型。

1. GLM Model=Binomial

1. Link function = log-odds

1. Training algo = IRLS

Statsmodels报告说，我们的模型有三个自由度。性别、Pclass和Age_Range，与预想相同的。 4. 自由度

对于二项式模型，statsmodels为你计算了三种拟合度量。最大对数可能性、偏差和皮尔逊齐次方。我们不会进一步检查它们，因为所有这三种措施只有在你比较两个或多个二项回归模型的拟合度时才有用，而在这种情况下，我们不是这样。

让我们看一下 拟合的系数

1. 拟合的系数

1. p values 所有的回归系数在0.1%的误差范围内都具有统计学意义，如P值都<0.001。

让我们看看每个系数都在告诉我们什么。

如何解释模型系数。

对于logit链接函数，拟合系数可以解释如下。

Age_Range: 它的系数是-0.0396。注意这个负值。解释这个系数值的方法是，在保持所有其他变量不变的情况下，乘客的年龄每增加一个单位，他们的生存几率就会 减少 一个系数=exp(-0.0396)=0.9564。也就是说，乘客的年龄每增加一个单位，就需要将他们的生存几率乘以0.9564，从而每次减少一定的生存几率。例如，如果一个12岁的二等舱男性乘客在灾难中的已知生存几率为8/9，那么一个22岁的二等舱男性乘客的生存几率为（8/9）* 0.9564¹⁰=大约6:10。

Pclass：Pclass的系数是-1.3095。再次注意到这个负号。乘客舱位等级的降低对乘客在泰坦尼克号上的生存几率有更大的影响。Pclass变量被编码为一等舱=1，二等舱=2，三等舱=3。因此，舱位等级每增加一个单位，即从头等舱到二等舱再到三等舱，在保持年龄和性别不变的情况下，生存几率就会减少一个系数exp(-1.3095)=0.30！也就是说，每降低一个单位，你的生存几率就会乘以0.30。例如，如果一个30岁的男性乘员在泰坦尼克号上有7：9的生存几率，只要把他降到2级，他的生存几率就减少到(7/9)/0.3=大约1：4。再往下跳到三等舱，其几率就降低到(7/9)0.30.3=7 : 100。

sex：最后，注意到性别变量的负系数非常大，为-2.5715。在沉没的泰坦尼克号上，与女性乘客相比，男性乘客的生存机会相当悲惨。在保持Pclass和年龄不变的情况下，男性乘客的生存几率仅为女性乘客的exp(-2.5715) = 7%。

预测

回顾一下，我们已经把测试数据集放在数据框_df_test_中。现在是时候在这个数据集上测试我们模型的性能了。要做到这一点，我们首先要在测试数据框中添加一个 Percentage Survived_列，我们将要求我们的模型预测这个值。

df_test['Pcnt_Survived'] = df_test['Survived']/df_test['Total']

我们将在结果对象上使用.predict()方法，并通过测试数据集获得预测的存活率。

predicted_survival_rate = binom_model_results.predict(X_test)

让我们绘制实际生存率与预测生存率的对比图。

import matplotlib.pyplot as plt
plt.xlabel('Actual Survival Rate')
plt.ylabel('Predicted Survival Rate')
plt.scatter(df_test['Pcnt_Survived'], 
   predicted_survival_rate, color = 'blue')plt.show()

实际与预测的存活率

正如你所看到的，当生存率接近范围的顶部，即1.0时，拟合就变得不可接受了。在很大程度上，预测的准确性是由样本大小决定的，即每组乘客的大小，由元组[Pclass, Sex, Age Range]分组。对于训练集中的某些组，组的规模太小，模型无法以有意义的方式进行训练。对于测试数据集中的这种组合，准确率会很低，这是可以理解的。

这里是本文中使用的泰坦尼克号数据链接集。

总结 

• 二项式回归模型可以用来预测一个事件的几率。
• 二项回归模型是广义线性模型家族中的一员，它使用一个合适的链接函数来建立响应变量   的条件期望与解释变量   的线性组合之间的关系。
• 在二项式回归模型中，我们通常使用对数函数作为链接函数。
• Logistic回归模型是二项式回归模型的一个特例，即数据集中每组解释变量的大小为1。