Statsmodels与R的survey/srvyr包加权逻辑回归结果差异的原因排查

我有一份虚构的加权调查数据集，包含受访者的汽车颜色和他们对“我喜欢开快车”这个问题的回答。我想做回归分析，看看轻微同意这个问题的可能性是否和受访者开黑色汽车有关。（这不是严肃的分析，只是用来对比R和Python中加权回归的输出差异。）

我先用R的survey和srvyr包做了加权逻辑回归，黑色汽车系数的检验统计量是-1.18，p值是0.238。但用statsmodels做加权逻辑回归时，这个系数的检验统计量是-1.35，p值是0.177。我想搞清楚为什么这两个检验统计量不一样，是不是我在其中某一个的设置里犯了错导致的差异。

值得一提的是，当我去掉两个测试中的权重部分后，检验统计量和p值几乎完全一致。所以看起来这两个实现对调查权重的处理方式不同。

下面是我的代码（注意：我用rpy2在同一个Notebook里运行R测试，这样就能和Python代码放在一起执行。只要你配置好了rpy2，在Jupyter Notebook里应该能复现我的输出。）

import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
%load_ext rpy2.ipython
%R library(dplyr)
%R library(srvyr)
%R library(survey)
%R library(broom)
import pandas as pd

df_car_survey = pd.read_csv(
    'https://raw.githubusercontent.com/ifstudies/\
carsurveydata/refs/heads/main/car_survey.csv')

# Adding dummy columns for independent and dependent variables:

for column in ['Car_Color', 'Enjoy_Driving_Fast']:
    df_car_survey = pd.concat([df_car_survey, pd.get_dummies(
        df_car_survey[column], dtype = 'int', prefix = column)],
             axis = 1)

df_car_survey.columns = [column.replace(' ', '_') for column in 
                         df_car_survey.columns]


# Loading DataFrame into R and creating a survey design object:
# See https://tidy-survey-r.github.io/tidy-survey-book/c10-sample-designs-replicate-weights.html
# for more details.
# This book was also inval
%Rpush df_car_survey
%R df_sdo <- df_car_survey %>% as_survey_design(\
    weights = 'Weight')

print("Survey design object:")
%R print(df_sdo)
    
# Logistic regression in R:
# (This code was based on that found in 
# https://tidy-survey-r.github.io/tidy-survey-book/c07-modeling.html )
%R logit_result <- svyglm(\
    Enjoy_Driving_Fast_Slightly_Agree ~ Car_Color_Black, \
design = df_sdo, na.action = na.omit,\
family = quasibinomial())
print("\n\n Logistic regression results from survey package within R:")
%R print(tidy(logit_result))

# Logistic regression within Python:
# (Based on StupidWolf's response at
# https://stackoverflow.com/a/62798889/13097194 )
glm = smf.glm("Enjoy_Driving_Fast_Slightly_Agree ~ Car_Color_Black",
        data = df_car_survey, 
        family = sm.families.Binomial(),
        freq_weights = df_car_survey['Weight'])
model_result = glm.fit()
print("\n\nStatsmodels logistic regression results:")
print(model_result.summary())

输出结果： #

Survey design object:
Independent Sampling design (with replacement)
Called via srvyr
Sampling variables:
  - ids: `1` 
  - weights: Weight 
Data variables: 
  - Car_Color (chr), Weight (dbl), Enjoy_Driving_Fast (chr), Count (int),
    Response_Sort_Map (int), Car_Color_Black (int), Car_Color_Red (int),
    Car_Color_White (int), Enjoy_Driving_Fast_Agree (int),
    Enjoy_Driving_Fast_Disagree (int), Enjoy_Driving_Fast_Slightly_Agree (int),
    Enjoy_Driving_Fast_Slightly_Disagree (int),
    Enjoy_Driving_Fast_Strongly_Agree (int),
    Enjoy_Driving_Fast_Strongly_Disagree (int)


 Logistic regression results from survey package within R:
# A tibble: 2 × 5
  term            estimate std.error statistic  p.value
  <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       -2.08      0.145    -14.3  8.85e-43
2 Car_Color_Black   -0.293     0.248     -1.18 2.38e- 1


Statsmodels logistic regression results:
                         Generalized Linear Model Regression Results                         
=============================================================================================
Dep. Variable:     Enjoy_Driving_Fast_Slightly_Agree   No. Observations:                 1059
Model:                                           GLM   Df Residuals:                     1057
Model Family:                               Binomial   Df Model:                            1
Link Function:                                 Logit   Scale:                          1.0000
Method:                                         IRLS   Log-Likelihood:                -345.24
Date:                               Tue, 04 Feb 2025   Deviance:                       690.48
Time:                                       10:07:53   Pearson chi2:                 1.06e+03
No. Iterations:                                    5   Pseudo R-squ. (CS):           0.001763
Covariance Type:                           nonrobust                                         
===================================================================================
                      coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------
Intercept          -2.0834      0.125    -16.693      0.000      -2.328      -1.839
Car_Color_Black    -0.2931      0.217     -1.350      0.177      -0.719       0.133
===================================================================================

备注：内容来源于stack exchange，提问作者KBurchfiel