在 Python 中重现 Stata 的标准错误

在 Python 中重现 Stata 的标准错误

最近，我一直在使用 Python 重新创建 Stata 估计值。下面，我重现了 Stata（和 R 的调查包）在其 prop 和 mean 命令中使用的稳健标准错误。我正在使用 Python 的 vega_datasets 中的人口数据。

_# 核心包_  
 **进口** 麻木的 **作为** np  
 **进口** 熊猫 **作为** PD  
 **从** vega_datasets **进口** 数据  
  
 _# 加载数据_  
 df = data.population()  
 df['女性'] = df['性别']-1  
 df = df[['年份', '女性', '年龄', '人']]  
  
 df.head()

加权标准误差

Statsmodels 是一个用于生成加权统计数据的常用 Python 库。它使用最常规的方法来计算加权标准偏差（注意：标准误差只是标准偏差除以样本量的平方根）。不幸的是，这不是 Stata 或 R 在计算稳健修正时计算加权标准偏差的方式。下面，Python 计算出 0.021 的标准误差，而 Stata 计算出 0.028。

**从** statsmodels.stats.weightstats **进口** 描述统计  
  
 stat = DescrStatsW(df['female'], weights=df['people'])  
 print('% 女性')  
 打印（'意思：'，stat.mean）  
 print('se:', stat.std/np.sqrt(df.shape[0])) **％ 女性  
 平均值：0.5046361028351961  
 与：0.0209417951337017**

状态码：

稳健的加权标准误差

幸运的是，statsmodels 可以使用加权线性回归计算相同的稳健加权标准误差。

**进口** 状态模型.api **作为** sm  
  
 y = df['女']  
 x = np.ones(df.shape[0])  
 w = df['人']  
 fpc = 1  
  
 _＃ 计算_  
 stat = sm.WLS(y, x, weights=w).fit(cov_type='HC1')  
  
 _# 统计_  
 平均值 = stat.params[0]  
 se = stat.bse[0]*fpc  
  
 print('% 女性')  
 打印（'意思：'，意思）  
 打印（'se：'，se） **％ 女性  
 平均值：0.5046361028351962  
 与：0.028250181556609938**

手动计算

这是如何运作的？好吧，我研究了来自 statsmodels 加权线性回归的原始代码，这样您就不必这样做了。下面是我们如何手动计算它。

 n = df.shape[0]  
 y = df['女']  
 x = df['人']  
 fpc = 1  
  
 _# 道具/平均值_  
 平均值 = (y*x).sum() / (x).sum()  
  
 _# 东南_  
 het_scale = n/(n-1) * ( ((y-mean) * np.sqrt(x))**2 )  
 wx = np.array([ x*(1/np.sqrt(x)) ]).T  
 pinv_wx = np.linalg.pinv(wx)  
 H = np.dot(pinv_wx, het_scale[:, None] * pinv_wx.T).item()  
 se = np.sqrt(H)*fpc  
  
 print('% 女性')  
 打印（'意思：'，意思）  
 打印（'se：'，se） **％ 女性  
 平均值：0.5046361028351961  
 与：0.028250181556609938**

数学有点超出本文，但这里是步骤的快速摘要：
1.异方差鲁棒协方差矩阵（HC1）
2. 通过 Cholesky 分解 (cholsigmainv) 美白
3. 矩阵的 Moore-Penrose 伪逆
4.异方差校正（或“White-Huber”）协方差矩阵

自定义函数

我们可以使用这段代码作为重写 Stata 的加权 prop 和 mean 函数的基础。我编写了 stata_func() 以允许 y_var 是二进制“0, 1”比例或连续平均值——但不是多项式比例。过变量（即子总体）可以是二元的、多项式的或不包括在内的。可以添加权重和FPC，但不能添加strata（这会大大扩大功能）。

_# 重新创建Stata统计信息_  
 **定义** stata_func(y_var, df, over= **没有任何** , 重量= **没有任何** , fpc=1):  
    _# 库_  
    **进口** 麻木的 **作为** np  
    **进口** 熊猫 **作为** PD  
    **从** category_encoders.one_hot **进口** OneHotEncoder  
      
    _# 如果 weight==None 则删除权重_  
    **如果** 重量是 **没有任何** ：  
 df['weight'] = np.ones(df.shape[0])  
 重量='重量'  
    _# 如果 over==None 则移除 over_  
    **如果** 结束了 **没有任何** ：  
 df['over'] = np.ones(df.shape[0])  
 结束 = '结束'  
    _# 准备要循环的变量_  
 one_hot = OneHotEncoder(handle_missing='return_nan', use_cat_names= **真的** ).fit_transform(df[over].astype('O'))  
 x_vars = 排序（one_hot.columns）  
 df = pd.concat([df, one_hot], axis=1)  
      
    _# 计算统计_  
 统计 = {}  
 stats['stats'] = ('mean', 'se', 'n')  
    **为了** x_vars 中的 x_var：  
          
        _# 文本输出_  
        **如果** 结束 == '结束'：  
 list_name = y_var  
        **别的** ：  
 list_name = y_var + " 当 " + x_var  
          
        _# 参数_  
 cw = df[[y_var, x_var, weight]].dropna()  
 n = cw.shape[0]  
 subpop = cw[cw[x_var]==1]  
 y = subpop[y_var]  
 x = subpop [权重]  
          
        _# 道具/平均值_  
 平均值 = (y*x).sum() / (x).sum()  
          
        _# 东南_  
        _##没有权重_  
        **如果** len(np.unique(x)) == 1:  
 het_scale = (y-mean)**2 _# 异方差鲁棒协方差矩阵（HC0）_  
        _##带权重_  
        **别的** ：  
 het_scale = n/(n-1) * ( ((y-mean) * np.sqrt(x))**2 ) _# 异方差鲁棒协方差矩阵（HC1）_  
        _## 鲁棒协方差矩阵上的一些线性代数_  
 wx = np.array([ x*(1/np.sqrt(x)) ]).T _# 通过 Cholesky 分解 (cholsigmainv) 变白_  
 pinv_wx = np.linalg.pinv(wx) _# Moore-Penrose 矩阵的伪逆_  
 H = np.dot(pinv_wx, het_scale[:, **没有任何** ] * pinv_wx.T).item() _# 异方差校正（或“White-Huber”）协方差矩阵_  
 se = np.sqrt(H)*fpc  
          
 stats[list_name] = (mean, se, str(n))  
          
 stats = pd.DataFrame(stats)  
      
    **返回** （统计）

以下是它的三个实际应用示例：

_# 添加子种群_  
 猫 = [-np.inf, 1900, 1950, 1990, np.inf]  
 cat_labs = [“1850-1900”、“1901-1950”、“1951-1990”、“1991-2000”]  
 df['period'] = pd.cut(df['year'], bins=cats, labels=cat_labs)  
  
 stata_func('female', df, over='period', weight='people')

 stata_func('年龄', df, over='周期', weight='人')

 stata_func('age', df, over=['period', 'female'], weight='people')

结论

为什么标准误不同？

让我们从讨论标准偏差开始。传统的标准偏差是通过相对于样本量（除以样本量然后取平方根）的总体方差（值的总和减去均方）来计算的。添加权重的最简单方法是将值减去均值乘以权重，并将样本大小替换为权重的总和。

然而，传统的标准偏差只考虑相对于 n 值（样本大小或权重总和）的总方差，而不考虑方差的分布。稳健标准差考虑方差的分布。在数学上，稳健的公式在将方差与 n 值相关之前不会对方差求和。相反，它考虑每个值与 n 值的方差，对较大的方差分布增加惩罚，并使用一些花哨的线性代数转换为单个标准差值。

这是做什么的（用简单的英语）？

稳健标准偏差应始终增加与传统公式的标准偏差。我们称其为“稳健修正”。如果每个值与平均值之间的差异非常不一致（有时很大，有时很小），则鲁棒校正应该更大。如果每个值与平均值之间的差异通常相同，则稳健校正不应使标准偏差增加太多。这在概念上与线性回归中的异方差性相同。可以在此处找到稳健标准偏差的公式：

MacKinnon, JG 和 White, H. (1985)。一些具有改进的有限样本属性的异方差一致协方差矩阵估计器。计量经济学杂志，29（3），305–325。

但是我们不关心错误而不是偏差吗？

标准误差只是按比例缩放的标准偏差（按比例缩放到原始值）。我们可以通过将标准差除以 n 值的平方根来得到标准误差。这适用于传统标准误和稳健标准误。

在我的博客上找到更多 Python 和 R 技巧： siebelm.github.io

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明

本文链接：https://www.qanswer.top/15484/42570509