11_二值選擇模型

王大桃zzZ發表於2024-05-03
模型

第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 其他離散選擇模型

略

相關文章