在 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 版权协议，转载请附上原文出处链接和本声明
浙公网安备 33010602011771号