第11章 二值选择模型

11.1 二值选择模型的例子

• 解释变量是离散的，不影响回归。
  • 比如虚拟变量
• 被解释变量是离散的，不适合进行OLS回归。
  • 离散选择模型、定性反应模型
  • 最常见的：二值选择行为

定义 线性概率模型（Linear Probility Model）$$\left { \begin{array}

P(y=1|x)= F(x,\beta) \
P(y=0|x)= 1-F(x,\beta)
\end{array} \right.$$

• 连接函数：\(F(x,\beta)\)
• 通过选择合适的链接函数，可以保证 \(0\le \hat y \le1\) 。例如：
  • Probit模型：正态分布的累积函数$$\Phi(x'\beta)\equiv \int_{- \infty}^{x'\beta}\phi(t)dt$$
  • Logit模型：逻辑分布的累积函数$$\Lambda(x'\beta)\equiv\frac{exp(x'\beta)}{1+exp(x'\beta)}$$

Logit模型优势：逻辑分布的CDF有解释表达式，计算方便，回归系数更具经济意义

11.2 最大似然估计的原理

非线性模型，无法通过变量转换转为线性模型，常使用最大似然估计法（Maxinum Likeihood Estimate）

定义 最大似然估计法（Maxinum Likeihood Estimate）

似然函数（Likehood function）$$L(\theta\ ;\ y_1,\cdots,y_n)=\prod_{i=1}^{n}f(y_i\ ;\ \theta)$$
对数似然函数（log-Likehood function）$$\ln L(\theta\ ;y_1,\cdots,y_n) = \sum_{i=1}^{n} \ln f(y_i \ ;\ \theta)$$
给定样本取值后，该样本最有可能来自参数 \(\theta\) 为何值的总体。换言之，寻找 \(\hat \theta_{ML}\) 使得观测到样本数据的可能性最大，即最大化对数似然函数：$$\max_{\theta} \ln L(\theta\ ;\ y_1,\cdots,y_n)$$
无约束极值问题的一阶条件为：$$\frac{\partial L(\theta\ ;y_1,\cdots,y_n)}{\partial \theta} = 0$$
求解此一阶条件，即可得到最大似然估计量 \(\hat \theta_{ML}\)

MLE估计量具有良好的大样本性质，可照常进行大样本统计推断。

• 是一致估计，\(p\lim_{n \to \infty} \hat\theta_{ML} = \theta\)
• 服从渐近正太分布
• 大样本下，渐近方差最小

非线性方程，通常没有解析解，只能寻找数值解。常用迭代法求数值解，高斯-牛顿法（Gauss-Newton Method）

• 解不唯一
• 求得的可能是局部最大值

11.3 二值选择模型的MLE估计

11.4 三种边际效应

线性模型回归系数：边际效应
非线性模型回归系数：通常不是常数，随x变化

• 平均边际效应
• 样本均值处边际效应
• 某代表值处边际效应

11.5 回归系数的经济意义

对于Logit模型，回归系数意味着：

• 几率(odds)：$$\hat\beta = \frac{p}{1-p}$$
  在实际运用中还可能运用到：
• 对数几率(log-odds):$$\ln(\hat\beta)=\ln(\frac{p}{1-p})$$
• 几率比(odds-ratio):$$exp(\hat\beta_j) = \frac{\frac{p^*}{1-P*}}{\frac{p}{1-p}}$$

11.6 拟合优度

• 准\(R^2\)(qusai-R-square)
• 正确预测百分比（precent correctly predicted）

11.7 准最大似然估计

11.8 三类渐近等价的大样本检验

• （1）沃尔德检验（wald test）
• （2）似然比检验（LR）
• （3）拉格朗日算子检验（LM）

11.9 二值选择模型的Stata命令及实例

[[Chapter_11.ipynb]]

1. 导入数据，查看各变量的统计特征

freq字段表示，数据出现的频次

• 这种类型的数据需要还原原始数据，不如会严重影响后续的回归结果。

import pandas as pd
import numpy as np
from cq import describe_bcmodel

# 读取数据
df = pd.read_stata('../2_Data/Data-2e/titanic.dta')
des = describe_bcmodel(df, frequency='freq')

为简化后续的步骤，借助ai生成了一个函数describe_bcmodel():

def describe_bcmodel(df, frequency, target=None, condition_col=None, condition=None):
    '''describe_bcmodel 二值模型的描述性统计，返回原始数据
    Arguments:
        df:dataframe -- 含有频次的数据集
        frequency:str --  频次的字段名
    Keyword Arguments:
        target:str -- 观测对象的字段名 (default: {None})
        condition_col:str -- 条件变量的字段名 (default: {None})
        condition:any -- 条件值 (default: {None})
    Returns:
         -- 按频次还原后的数据集
    '''
    result = df.loc[np.repeat(df.index.values, df[frequency])].drop(frequency, axis=1).reset_index(drop=True)
    if (target is None) and (condition_col is None):
        print(result.describe().T)
    else:
        result = result[[target, condition_col]][result[condition_col] == condition].drop(condition_col, axis=1).reset_index(drop=True)
        print(f'when {condition_col} is {condition}:')
        print(result.describe().T)
        return result
    return result

函数说明
np.repeat(df.index.values, df['freq'])：

• df.index.values 返回数据框的索引值，这是一个代表行号的数组。
• df['freq'] 返回'freq'列的值，这是一个代表信息重复次数的数组。
• np.repeat函数将行索引根据'freq'列的值进行重复，以便在最终结果中重复出现对应次数。

df.loc[]：

• df.loc 是用于按标签选择行和列的方法。在这里，它使用重复后的索引来选择数据框中的行。

.drop('freq', axis=1)：

• drop 方法用于删除数据框中的列。在这里，它删除了名为'freq'的列。参数axis=1表示删除列

.reset_index(drop=True)：

• reset_index 方法用于重置索引。参数drop=True表示删除原始索引，使新索引从零开始。这样可以确保最终结果的索引是连续的整数序列。

综合起来，这行代码的作用是将数据框中的行根据'freq'列的值重复多次，然后丢弃'freq'列，并重置索引，以得到非'freq'列的信息按照出现次数重复的结果。

2.观察不同特征下的存活率

for col in des.drop('survive', axis=1).columns:
    describe_bcmodel(df,
                     'freq',
                     target='survive',
                     condition_col=col,
                     condition=1)

3.构建OLS参照系

import statsmodels.api as sm

X = des[['class1','class2','class3','child','female']]
y = des['survive']
X = sm.add_constant(X)
model_ols = sm.OLS(y,X)
results_ols = model_ols.fit()
print(results_ols.summary())

4.使用Logit模型进行估计

sm.logit(endog, exog).fit(disp=0)

• disp = 0 不现实迭代过程，只显示结果
• disp = 1 显示迭代过程

model_logit = sm.Logit(y, X)
result_logit = model_logit.fit(disp=1)
print(result_logit.summary())

5.使用稳健标准误进行Logit估计

# print(model_logit.fit(cov_type='HC0').summary())
# print(model_logit.fit(cov_type='HC1').summary())
# print(model_logit.fit(cov_type='HC2').summary())
print(model_logit.fit(cov_type='HC3').summary())

6.显示Logit回归的几率比

import numpy as np

odds_ratios = np.exp(result_logit.params)
result_logit_or = pd.DataFrame({'odds ratio': odds_ratios,
                                'std err': result_logit.bse,
                                'z':result_logit.tvalues,
                                'p>|z|':result_logit.pvalues,
                                },
                               index=result_logit.params.index)
pd.set_option('display.float_format', '{:.4f}'.format)
result_logit_or

7.计算Logit模型的平均边际效应

mfx = result_logit.get_margeff()
print(mfx.summary())

8.计算均值处的平均边际效应

mfx = result_logit.get_margeff(at='mean')
print(mfx.summary())

9.准确度测量

用模型预测值与实际值进行比较，计算预测值与实际值相符的比例

predicted_classes = result_logit.predict(X) > 0.5
# 计算准确率
accuracy = (predicted_classes == y).mean()
print(f"Accuracy of the model: {accuracy*100:.2f}%")

10.数据预测

msrose = pd.DataFrame([1,1,0,0,0,1],
                      index=result_logit.params.index,columns=['MS-ROSE'])
# 两种不同的赋值方式
mrjack = pd.DataFrame({'const':1,
                       'class1':0,
                       'class2':0,
                       'class3':1,
                       'child':0,
                       'female':0},
                      index=['MR-Jack'],columns=result_logit.params.index.T)

print(result_logit.predict(msrose.T))
print(result_logit.predict(mrjack))

11.使用Probit模型进行回归

model_probit = sm.Probit(y,X)
results_probit = model_probit.fit()
print(results_probit.summary())

# 计算边际效用
mfx_probit = results_probit.get_margeff()
print(mfx_probit.summary())

# 计算准确率
predicted_classes_probit = results_probit.predict(X) > 0.5
accuracy = (predicted_classes_probit == y).mean()
print(f"Accuracy of the model: {accuracy*100:.2f}%")

# 对比 logit 和 probit 模型
df = pd.DataFrame(np.corrcoef(predicted_classes,predicted_classes_probit),index=['logit','probit'],columns=['logit','probit'])
df

11.10 其他离散选择模型

略