Python 中是否有等效的 R-svydesign 来应用复杂的调查设计权重？

Question

有没有办法将复杂调查设计中的调查权重结合起来进行描述性统计（例如频率和交叉表）以及线性/逻辑回归和建模？

分析具有复杂调查设计的调查数据集通常涉及分层和聚类以代表更大的人群（例如，CDC 国家健康访谈调查。第 41 页）。

几篇文章Python 中的设计校正方差估计和如何在 Python 中使用分层进行 SAS Proc SurveyFreq提及 samplics、Quantipy 和 PandaSurvey 包。

不幸的是，这些软件包都无法真正实现调查数据分析，并且文档非常有限。

在 R 中，我知道这可以通过调查包实现，通过以下方式将调查设计应用于数据框：

# Readin data
df = read_dta('data.dta')

# Define survey design
design <- svydesign(id=~ID, weights=~analwt ,strata=~Final_strata, data=df)

# Apply survey weights
df_weighted <- svydesign(id=~ID, weights=~analwt, strata=~Final_strata, nest=TRUE, survey.lonely.psu="adjust",data=df)

# Crosstab
xtab = svytable(~var1+var2, df_weighted) %>%
  prop.table(1)*100

在 SAS 中，一个例子是：

proc surveyfreq data=PS_Data;
tables INT4a; strata
final_strata; cluster granteeid;
weight analwt;
run;

在 STATA 中，这是通过以下方式完成的：

Svyset GranteeID [pweight=ANALWT], strata(Final_strata) vce(linearized)
Then use the svy prefix for analysis commands.

我已经尝试过 pandas 中的交叉表功能，但它似乎无法解释 strat/clustering。

xtab_weighted = pd.crosstab(
    df['var1'], df['var2'], 
    df.analwt,
    aggfunc=sum, dropna=True, normalize=True)

Answer 1

这是 R 中的代码块，通过将 pd.crosstab() 中的标准化设置为“列”而不是“索引”，这似乎是一个虚拟错误。

# Load libraries 
library(pacman)
p_load(haven, survey, gmodels)

# Read data
df = read_dta('hcps_puf_2022.dta')

# Apply survey design
design <- svydesign(id=~GranteeID, weights=~analwt, strata=~Final_strata, data=df)

# Create the weighted dataset
df_weighted <- svydesign(id=~GranteeID, weights=~analwt,strata=~Final_strata,nest=TRUE,survey.lonely.psu="adjust",data=df)

# Print an unweighted crosstab
table = table(df$sex_at_birth, df$MED1)
print(table)
      -2   -1    1    2
  1    1    2  968  504
  2    0    5 1928 1000

prop = prop.table(table)
print(prop*100)
             -2          -1           1           2
  1  0.02268603  0.04537205 21.96007260 11.43375681
  2  0.00000000  0.11343013 43.73865699 22.68602541

# Print a weighted crosstab
xtab = svytable(~sex_at_birth+MED1, df_weighted) %>%
  prop.table(1)*100
print(xtab)
            MED1
sex_at_birth          -2          -1           1           2
           1  0.57676525  0.01002916 58.60754564 40.80565996
           2  0.00000000  0.05332145 62.64414195 37.30253660

Python 代码似乎可以解决这个问题，即将 Normalize='columns' 更改为 Normalize='index'。

import pandas as pd
import pyreadstat

# Load the data
df, meta = pyreadstat.read_dta ('hcps_puf_2022.dta')

# Print unweighted crosstab
table = pd.crosstab(df['sex_at_birth'], df['MED1'])
print(table)
MED1          -2  -1     1     2
sex_at_birth                    
1              1   2   968   504
2              0   5  1928  1000

prop = pd.crosstab(df['sex_at_birth'], df['MED1'],
                   normalize='all')
print(prop*100)
MED1                -2        -1          1          2
sex_at_birth                                          
1             0.022686  0.045372  21.960073  11.433757
2             0.000000  0.113430  43.738657  22.686025

# Print weighted crosstab
xtab = pd.crosstab(df['sex_at_birth'], df['MED1'], df['analwt'], 
                   dropna=True, aggfunc=sum, normalize='index')
print(xtab*100)
MED1                -2        -1          1          2
sex_at_birth                                          
1             0.576765  0.010029  58.607546  40.805660
2             0.000000  0.053321  62.644142  37.302537