代码单元 1
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
泊松响应数据的加权 GLM
载入数据
这个例子使用 Fair 数据集,通过少量解释变量建立模型,考察婚外关系相关指标。这里的重点是权重与数据组织方式如何影响估计结果,并非因果推断,也不据此判断某个具体人的行为。
原 notebook 构造权重,说明 freq_weights 如何对应重复观测;再比较 var_weights 在聚合数据中的作用。两者不能仅凭名称相近而互换。
后面的输出均保留自官方 notebook 的已保存结果,属于原文示例输出。变量名和模型输出字段保留英文,便于与代码一一对应。
代码单元 2
print(sm.datasets.fair.NOTE)
原文输出(单元 2):
::
Number of observations: 6366
Number of variables: 9
Variable name definitions:
rate_marriage : How rate marriage, 1 = very poor, 2 = poor, 3 = fair,
4 = good, 5 = very good
age : Age
yrs_married : No. years married. Interval approximations. See
original paper for detailed explanation.
children : No. children
religious : How religious, 1 = not, 2 = mildly, 3 = fairly,
4 = strongly
educ : Level of education, 9 = grade school, 12 = high
school, 14 = some college, 16 = college graduate,
17 = some graduate school, 20 = advanced degree
occupation : 1 = student, 2 = farming, agriculture; semi-skilled,
or unskilled worker; 3 = white-collar; 4 = teacher
counselor social worker, nurse; artist, writers;
technician, skilled worker, 5 = managerial,
administrative, business, 6 = professional with
advanced degree
occupation_husb : Husband's occupation. Same as occupation.
affairs : measure of time spent in extramarital affairs
See the original paper for more details.
将数据载入 pandas DataFrame。
代码单元 3
data = sm.datasets.fair.load_pandas().data
因变量(内生变量)是 affairs。
代码单元 4
data.describe()
原文输出(单元 4):
| rate_marriage | age | yrs_married | children | religious | educ | occupation | occupation_husb | affairs | |
|---|---|---|---|---|---|---|---|---|---|
| count | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 | 6366.000000 |
| mean | 4.109645 | 29.082862 | 9.009425 | 1.396874 | 2.426170 | 14.209865 | 3.424128 | 3.850141 | 0.705374 |
| std | 0.961430 | 6.847882 | 7.280120 | 1.433471 | 0.878369 | 2.178003 | 0.942399 | 1.346435 | 2.203374 |
| min | 1.000000 | 17.500000 | 0.500000 | 0.000000 | 1.000000 | 9.000000 | 1.000000 | 1.000000 | 0.000000 |
| 25% | 4.000000 | 22.000000 | 2.500000 | 0.000000 | 2.000000 | 12.000000 | 3.000000 | 3.000000 | 0.000000 |
| 50% | 4.000000 | 27.000000 | 6.000000 | 1.000000 | 2.000000 | 14.000000 | 3.000000 | 4.000000 | 0.000000 |
| 75% | 5.000000 | 32.000000 | 16.500000 | 2.000000 | 3.000000 | 16.000000 | 4.000000 | 5.000000 | 0.484848 |
| max | 5.000000 | 42.000000 | 23.000000 | 5.500000 | 4.000000 | 20.000000 | 6.000000 | 6.000000 | 57.599991 |
代码单元 5
data[:3]
原文输出(单元 5):
| rate_marriage | age | yrs_married | children | religious | educ | occupation | occupation_husb | affairs | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 3.0 | 32.0 | 9.0 | 3.0 | 3.0 | 17.0 | 2.0 | 5.0 | 0.111111 |
| 1 | 3.0 | 27.0 | 13.0 | 3.0 | 1.0 | 14.0 | 3.0 | 4.0 | 3.230769 |
| 2 | 4.0 | 22.0 | 2.5 | 0.0 | 1.0 | 16.0 | 3.0 | 5.0 | 1.400000 |
接下来主要使用泊松 GLM。虽然程序能够对小数形式的响应进行拟合,原例为了演示计数响应,将 affairs 向上取整。这是教学中的数据变换,不意味着原始指标天然就是符合泊松分布的事件次数。
代码单元 6
data["affairs"] = np.ceil(data["affairs"])
data[:3]
原文输出(单元 6):
| rate_marriage | age | yrs_married | children | religious | educ | occupation | occupation_husb | affairs | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 3.0 | 32.0 | 9.0 | 3.0 | 3.0 | 17.0 | 2.0 | 5.0 | 1.0 |
| 1 | 3.0 | 27.0 | 13.0 | 3.0 | 1.0 | 14.0 | 3.0 | 4.0 | 4.0 |
| 2 | 4.0 | 22.0 | 2.5 | 0.0 | 1.0 | 16.0 | 3.0 | 5.0 | 2.0 |
代码单元 7
(data["affairs"] == 0).mean()
原文输出(单元 7):
np.float64(0.6775054979579014)
代码单元 8
np.bincount(data["affairs"].astype(int))
原文输出(单元 8):
array([4313, 934, 488, 180, 130, 172, 7, 21, 67, 2, 0,
0, 17, 0, 0, 0, 3, 12, 8, 0, 0, 0,
0, 0, 2, 2, 2, 3, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1])
压缩与聚合观测
原始数据有 6366 条观测。如果只考虑选定的变量,其中不少行会完全相同。下面采用两种方式合并:第一种要求响应与解释变量全部相同;第二种只要求解释变量相同。
完全相同观测的压缩数据
用 pandas 的 groupby 合并所选变量完全相同的行,再创建 freq 列,记录每一行代表了多少条原始观测。
代码单元 9
data2 = data.copy()
data2["const"] = 1
dc = (
data2["affairs rate_marriage age yrs_married const".split()]
.groupby("affairs rate_marriage age yrs_married".split())
.count()
)
dc = dc.reset_index()
dc = dc.rename(columns={"const": "freq"})
print(dc.shape)
dc.head()
原文输出(单元 9):
(476, 5)
原文输出(单元 9):
| affairs | rate_marriage | age | yrs_married | freq | |
|---|---|---|---|---|---|
| 0 | 0.0 | 1.0 | 17.5 | 0.5 | 1 |
| 1 | 0.0 | 1.0 | 22.0 | 2.5 | 3 |
| 2 | 0.0 | 1.0 | 27.0 | 2.5 | 1 |
| 3 | 0.0 | 1.0 | 27.0 | 6.0 | 5 |
| 4 | 0.0 | 1.0 | 27.0 | 9.0 | 1 |
解释变量相同的聚合数据
接下来只按解释变量分组。同组观测的响应可能不同,因此同时计算响应的均值、总和与数量。
继续使用 pandas 的 groupby 和聚合操作,再把产生的多层列索引 MultiIndex 展平成普通列名。
代码单元 10
gr = data["affairs rate_marriage age yrs_married".split()].groupby(
"rate_marriage age yrs_married".split()
)
df_a = gr.agg(["mean", "sum", "count"])
def merge_tuple(tpl):
if isinstance(tpl, tuple) and len(tpl) > 1:
return "_".join(map(str, tpl))
else:
return tpl
df_a.columns = df_a.columns.map(merge_tuple)
df_a.reset_index(inplace=True)
print(df_a.shape)
df_a.head()
原文输出(单元 10):
(130, 6)
原文输出(单元 10):
| rate_marriage | age | yrs_married | affairs_mean | affairs_sum | affairs_count | |
|---|---|---|---|---|---|---|
| 0 | 1.0 | 17.5 | 0.5 | 0.000000 | 0.0 | 1 |
| 1 | 1.0 | 22.0 | 2.5 | 3.900000 | 39.0 | 10 |
| 2 | 1.0 | 27.0 | 2.5 | 3.400000 | 17.0 | 5 |
| 3 | 1.0 | 27.0 | 6.0 | 0.900000 | 9.0 | 10 |
| 4 | 1.0 | 27.0 | 9.0 | 1.333333 | 4.0 | 3 |
合并后,dc 包含 476 条不同的观测组合,df_a 包含 130 条不同的解释变量组合。原文此处写成 467,与它紧邻的代码输出不一致;这里依据原文输出改为 476。
代码单元 11
print("number of rows: \noriginal, with unique observations, with unique exog")
data.shape[0], dc.shape[0], df_a.shape[0]
原文输出(单元 11):
number of rows: original, with unique observations, with unique exog
原文输出(单元 11):
(6366, 476, 130)
模型分析
下面比较原始数据的泊松 GLM,以及用权重或暴露量表示重复与聚合观测的模型。
原始数据
代码单元 12
glm = smf.glm(
"affairs ~ rate_marriage + age + yrs_married",
data=data,
family=sm.families.Poisson(),
)
res_o = glm.fit()
print(res_o.summary())
原文输出(单元 12):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs No. Observations: 6366
Model: GLM Df Residuals: 6362
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -10351.
Date: Thu, 27 Aug 2026 Deviance: 15375.
Time: 06:50:44 Pearson chi2: 3.23e+04
No. Iterations: 6 Pseudo R-squ. (CS): 0.2420
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
代码单元 13
res_o.pearson_chi2 / res_o.df_resid
原文输出(单元 13):
np.float64(5.078702313363289)
用频数权重拟合压缩数据
合并完全相同的观测后,使用 freq_weights 表示每条记录的重复次数,会得到相同的参数估计与相关推断。不过,面向单条存储记录的结果属性可能不同。例如,残差数组不会因为频数权重而自动展开为原始数据的行数。
代码单元 14
glm = smf.glm(
"affairs ~ rate_marriage + age + yrs_married",
data=dc,
family=sm.families.Poisson(),
freq_weights=np.asarray(dc["freq"]),
)
res_f = glm.fit()
print(res_f.summary())
原文输出(单元 14):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs No. Observations: 476
Model: GLM Df Residuals: 6362
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -10351.
Date: Thu, 27 Aug 2026 Deviance: 15375.
Time: 06:50:44 Pearson chi2: 3.23e+04
No. Iterations: 6 Pseudo R-squ. (CS): 0.9754
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
代码单元 15
res_f.pearson_chi2 / res_f.df_resid
原文输出(单元 15):
np.float64(5.078702313363198)
将频数改传给 var_weights
接下来比较 var_weights 和 freq_weights。方差权重通常用于响应为均值的情形,而不是表示完全相同观测的重复次数。
原作者指出,不能据本例直接断言两种权重在一般情形下有相同理论结果。本例确实得到相同的参数估计,但残差自由度 df_resid 不同,因为 var_weights 不会像频数权重那样增加有效观测数。
代码单元 16
glm = smf.glm(
"affairs ~ rate_marriage + age + yrs_married",
data=dc,
family=sm.families.Poisson(),
var_weights=np.asarray(dc["freq"]),
)
res_fv = glm.fit()
print(res_fv.summary())
原文输出(单元 16):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs No. Observations: 476
Model: GLM Df Residuals: 472
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -10351.
Date: Thu, 27 Aug 2026 Deviance: 15375.
Time: 06:50:44 Pearson chi2: 3.23e+04
No. Iterations: 6 Pseudo R-squ. (CS): 0.9754
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
如果本意是恢复原始重复观测的离散程度,直接用这一方差权重模型的 df_resid 作分母,会得到与原始分析不对应的结果。采用原始观测的残差自由度,才得到原例预期的离散程度。这里是在比较不同权重语义,不应把两者自由度不同本身视为程序错误。
代码单元 17
res_fv.pearson_chi2 / res_fv.df_resid, res_f.pearson_chi2 / res_f.df_resid
原文输出(单元 17):
(np.float64(68.45488160512006), np.float64(5.078702313363198))
解释变量相同的总和与均值
现在按解释变量相同的观测进行聚合,响应取组内总和或均值。
总和配合 exposure
若响应是组内所有响应的总和,在相应的独立泊松模型假设下,总和仍为泊松分布。每组包含的人数不同,因此用 exposure 表示每条聚合记录代表的观测数量。
这里的参数估计与参数协方差与原始数据相同,但对数似然、偏差和 Pearson 卡方统计量有所不同。
代码单元 18
glm = smf.glm(
"affairs_sum ~ rate_marriage + age + yrs_married",
data=df_a,
family=sm.families.Poisson(),
exposure=np.asarray(df_a["affairs_count"]),
)
res_e = glm.fit()
print(res_e.summary())
原文输出(单元 18):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs_sum No. Observations: 130
Model: GLM Df Residuals: 126
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -740.75
Date: Thu, 27 Aug 2026 Deviance: 967.46
Time: 06:50:44 Pearson chi2: 926.
No. Iterations: 6 Pseudo R-squ. (CS): 1.000
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
代码单元 19
res_e.pearson_chi2 / res_e.df_resid
原文输出(单元 19):
np.float64(7.350789109179558)
均值配合 var_weights
也可以把组内响应均值作为因变量。这时,均值的方差与该组总暴露量成反比,因此用代表组内数量的方差权重。
代码单元 20
glm = smf.glm(
"affairs_mean ~ rate_marriage + age + yrs_married",
data=df_a,
family=sm.families.Poisson(),
var_weights=np.asarray(df_a["affairs_count"]),
)
res_a = glm.fit()
print(res_a.summary())
原文输出(单元 20):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs_mean No. Observations: 130
Model: GLM Df Residuals: 126
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -5954.2
Date: Thu, 27 Aug 2026 Deviance: 967.46
Time: 06:50:44 Pearson chi2: 926.
No. Iterations: 5 Pseudo R-squ. (CS): 1.000
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
对比四种结果
前面的汇总显示,params、cov_params 及相关 Wald 推断一致。下面逐项比较不同数据表示下的结果属性。
参数估计 params、参数标准误 bse,以及检验参数为零的 pvalues 都一致。然而,对数似然 llf、偏差 deviance 和 pearson_chi2 只部分一致。聚合后的统计量与原始逐条观测的结果不同。
原文提醒:这些统计量的处理方式在以后版本中仍可能改变。
当模型设定适当时,对解释变量相同的响应取总和或均值,都可以建立合理的似然解释。不过,当前这些结果统计量没有完整反映这种对应关系;计算上可能涉及聚合数据所需的调整。
理论上还需要区分另一类情形:特别是把 var_weights 当作方差调整、而分布模型并不完全正确时,普通似然分析可能不适用,应从准似然角度解释估计。方差权重既可以表示一个正确指定模型中的均值精度,也可以只是准似然模型中的方差调整,因此存在语义差异。原作者在这里并不尝试使所有似然定义完全匹配。
下一节将展示:在假定底层模型设定正确的前提下,各种聚合表示仍会得到相同的似然比型比较结果。不能因某一组绝对统计量不同,就直接跨表示比较模型优劣。
代码单元 21
results_all = [res_o, res_f, res_e, res_a]
names = "res_o res_f res_e res_a".split()
代码单元 22
pd.concat([r.params for r in results_all], axis=1, keys=names)
原文输出(单元 22):
| res_o | res_f | res_e | res_a | |
|---|---|---|---|---|
| Intercept | 2.715533 | 2.715533 | 2.715533 | 2.715533 |
| rate_marriage | -0.495180 | -0.495180 | -0.495180 | -0.495180 |
| age | -0.029914 | -0.029914 | -0.029914 | -0.029914 |
| yrs_married | -0.010763 | -0.010763 | -0.010763 | -0.010763 |
代码单元 23
pd.concat([r.bse for r in results_all], axis=1, keys=names)
原文输出(单元 23):
| res_o | res_f | res_e | res_a | |
|---|---|---|---|---|
| Intercept | 0.107360 | 0.107360 | 0.107360 | 0.107360 |
| rate_marriage | 0.011874 | 0.011874 | 0.011874 | 0.011874 |
| age | 0.004471 | 0.004471 | 0.004471 | 0.004471 |
| yrs_married | 0.004294 | 0.004294 | 0.004294 | 0.004294 |
代码单元 24
pd.concat([r.pvalues for r in results_all], axis=1, keys=names)
原文输出(单元 24):
| res_o | res_f | res_e | res_a | |
|---|---|---|---|---|
| Intercept | 3.756282e-141 | 3.756280e-141 | 3.756282e-141 | 3.756282e-141 |
| rate_marriage | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| age | 2.221918e-11 | 2.221918e-11 | 2.221918e-11 | 2.221918e-11 |
| yrs_married | 1.219200e-02 | 1.219200e-02 | 1.219200e-02 | 1.219200e-02 |
代码单元 25
pd.DataFrame(
np.column_stack([[r.llf, r.deviance, r.pearson_chi2] for r in results_all]),
columns=names,
index=["llf", "deviance", "pearson chi2"],
)
原文输出(单元 25):
| res_o | res_f | res_e | res_a | |
|---|---|---|---|---|
| llf | -10350.913296 | -10350.913296 | -740.748534 | -5954.219866 |
| deviance | 15374.679054 | 15374.679054 | 967.455734 | 967.455734 |
| pearson chi2 | 32310.704118 | 32310.704118 | 926.199428 | 926.199428 |
似然比型检验
前面已经看到,聚合数据与原始数据的似然和部分拟合优度统计量并不相同。下面展示:本例的似然差与偏差差在各个数据表示下相同,但 Pearson 卡方的差并不相同。原作者强调,这一部分仍有需要澄清之处,未来实现可能改变。
具体做法是去掉 age,形成约束较多的简化模型,再计算“简化模型减去完整模型”的统计量差。代码最后一个数是对数似然差;常规似然比统计量是它的 -2 倍,不能把负的对数似然差直接当作卡方统计量。
原始观测与频数权重
代码单元 26
glm = smf.glm(
"affairs ~ rate_marriage + yrs_married", data=data, family=sm.families.Poisson()
)
res_o2 = glm.fit()
# print(res_f2.summary())
res_o2.pearson_chi2 - res_o.pearson_chi2, res_o2.deviance - res_o.deviance, res_o2.llf - res_o.llf
原文输出(单元 26):
(np.float64(52.91343161856639), np.float64(45.726693322505525), np.float64(-22.863346661253672))
代码单元 27
glm = smf.glm(
"affairs ~ rate_marriage + yrs_married",
data=dc,
family=sm.families.Poisson(),
freq_weights=np.asarray(dc["freq"]),
)
res_f2 = glm.fit()
# print(res_f2.summary())
res_f2.pearson_chi2 - res_f.pearson_chi2, res_f2.deviance - res_f.deviance, res_f2.llf - res_f.llf
原文输出(单元 27):
(np.float64(52.913431618668255), np.float64(45.726693322505525), np.float64(-22.863346661251853))
聚合数据:exposure 与 var_weights
本例的似然比型比较仍与原始观测一致,但 pearson_chi2 的差值不仅不同,还出现了原作者所关注的负号。
代码单元 28
glm = smf.glm(
"affairs_sum ~ rate_marriage + yrs_married",
data=df_a,
family=sm.families.Poisson(),
exposure=np.asarray(df_a["affairs_count"]),
)
res_e2 = glm.fit()
res_e2.pearson_chi2 - res_e.pearson_chi2, res_e2.deviance - res_e.deviance, res_e2.llf - res_e.llf
原文输出(单元 28):
(np.float64(-31.61852752510879), np.float64(45.72669332250621), np.float64(-22.86334666125299))
代码单元 29
glm = smf.glm(
"affairs_mean ~ rate_marriage + yrs_married",
data=df_a,
family=sm.families.Poisson(),
var_weights=np.asarray(df_a["affairs_count"]),
)
res_a2 = glm.fit()
res_a2.pearson_chi2 - res_a.pearson_chi2, res_a2.deviance - res_a.deviance, res_a2.llf - res_a.llf
原文输出(单元 29):
(np.float64(-31.618527525113905), np.float64(45.72669332250621), np.float64(-22.863346661252763))
检查 Pearson 卡方统计量
先进行一些基本一致性检查,确认 pearson_chi2 与 Pearson 残差 resid_pearson 的计算没有明显矛盾。以下单元访问了 _results 等内部属性,只保留为原 notebook 的诊断过程;它们不应被当作稳定的业务接口依赖。
代码单元 30
res_e2.pearson_chi2, res_e.pearson_chi2, (res_e2.resid_pearson**2).sum(), (
res_e.resid_pearson**2
).sum()
原文输出(单元 30):
(np.float64(894.5809002315154), np.float64(926.1994277566242), np.float64(894.5809002315157), np.float64(926.1994277566242))
代码单元 31
res_e._results.resid_response.mean(), res_e.model.family.variance(res_e.mu)[
:5
], res_e.mu[:5]
原文输出(单元 31):
(np.float64(-2.4049138748803392e-15), array([ 5.42753476, 46.42940306, 19.98971769, 38.50138978, 11.18341883]), array([ 5.42753476, 46.42940306, 19.98971769, 38.50138978, 11.18341883]))
代码单元 32
(res_e._results.resid_response**2 / res_e.model.family.variance(res_e.mu)).sum()
原文输出(单元 32):
np.float64(926.1994277566242)
代码单元 33
res_e2._results.resid_response.mean(), res_e2.model.family.variance(res_e2.mu)[
:5
], res_e2.mu[:5]
原文输出(单元 33):
(np.float64(3.1045251839364374e-14), array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538]), array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538]))
代码单元 34
(res_e2._results.resid_response**2 / res_e2.model.family.variance(res_e2.mu)).sum()
原文输出(单元 34):
np.float64(894.5809002315154)
代码单元 35
(res_e2._results.resid_response**2).sum(), (res_e._results.resid_response**2).sum()
原文输出(单元 35):
(np.float64(51204.85737832324), np.float64(47104.64779595964))
差值出现负号的一个可能原因是:相减的两个平方项使用了不同的方差分母。相关分析中,可以比较在完整模型与简化模型里使用同一个方差假设的结果。
下面统一使用简化模型的方差作为分母。本例中,四种数据表示下得到相同的缩放平方差。后续讨论见官方 issue #3616。这里是一项诊断性比较,并未给出一个可普遍替代似然比检验的新检验。
代码单元 36
(
(res_e2._results.resid_response**2 - res_e._results.resid_response**2)
/ res_e2.model.family.variance(res_e2.mu)
).sum()
原文输出(单元 36):
np.float64(44.43314175121902)
代码单元 37
(
(res_a2._results.resid_response**2 - res_a._results.resid_response**2)
/ res_a2.model.family.variance(res_a2.mu)
* res_a2.model.var_weights
).sum()
原文输出(单元 37):
np.float64(44.43314175121905)
代码单元 38
(
(res_f2._results.resid_response**2 - res_f._results.resid_response**2)
/ res_f2.model.family.variance(res_f2.mu)
* res_f2.model.freq_weights
).sum()
原文输出(单元 38):
np.float64(44.43314175122029)
代码单元 39
(
(res_o2._results.resid_response**2 - res_o._results.resid_response**2)
/ res_o2.model.family.variance(res_o2.mu)
).sum()
原文输出(单元 39):
np.float64(44.43314175121979)
补充检查
notebook 的剩余部分是额外的一致性检查,首次阅读可以略过。这里完整保留,便于按顺序核对。
其中出现的固定卡方输入值是原笔记中的附加计算,不应自动当作前面模型比较得到的统计量;请分别查看变量来源。
代码单元 40
np.exp(res_e2.model.exposure)[:5], np.asarray(df_a["affairs_count"])[:5]
原文输出(单元 40):
(array([ 1., 10., 5., 10., 3.]), array([ 1, 10, 5, 10, 3]))
代码单元 41
res_e2.resid_pearson.sum() - res_e.resid_pearson.sum()
原文输出(单元 41):
np.float64(-9.66481794586187)
代码单元 42
res_e2.mu[:5]
原文输出(单元 42):
array([ 4.77165474, 44.4026604 , 22.2013302 , 39.14749309, 10.54229538])
代码单元 43
res_a2.pearson_chi2, res_a.pearson_chi2, res_a2.resid_pearson.sum(), res_a.resid_pearson.sum()
原文输出(单元 43):
(np.float64(894.5809002315154), np.float64(926.1994277566293), np.float64(-42.34720713518796), np.float64(-32.68238918932164))
代码单元 44
(
(res_a2._results.resid_response**2)
/ res_a2.model.family.variance(res_a2.mu)
* res_a2.model.var_weights
).sum()
原文输出(单元 44):
np.float64(894.5809002315154)
代码单元 45
(
(res_a._results.resid_response**2)
/ res_a.model.family.variance(res_a.mu)
* res_a.model.var_weights
).sum()
原文输出(单元 45):
np.float64(926.1994277566293)
代码单元 46
(
(res_a._results.resid_response**2)
/ res_a.model.family.variance(res_a2.mu)
* res_a.model.var_weights
).sum()
原文输出(单元 46):
np.float64(850.1477584802964)
代码单元 47
res_e.model.endog[:5], res_e2.model.endog[:5]
原文输出(单元 47):
(array([ 0., 39., 17., 9., 4.]), array([ 0., 39., 17., 9., 4.]))
代码单元 48
res_a.model.endog[:5], res_a2.model.endog[:5]
原文输出(单元 48):
(array([0. , 3.9 , 3.4 , 0.9 , 1.33333333]), array([0. , 3.9 , 3.4 , 0.9 , 1.33333333]))
代码单元 49
res_a2.model.endog[:5] * np.exp(res_e2.model.exposure)[:5]
原文输出(单元 49):
array([ 0., 39., 17., 9., 4.])
代码单元 50
res_a2.model.endog[:5] * res_a2.model.var_weights[:5]
原文输出(单元 50):
array([ 0., 39., 17., 9., 4.])
代码单元 51
from scipy import stats
stats.chi2.sf(27.19530754604785, 1), stats.chi2.sf(29.083798806764687, 1)
原文输出(单元 51):
(np.float64(1.8390448369994542e-07), np.float64(6.931421143170174e-08))
代码单元 52
res_o.pvalues
原文输出(单元 52):
Intercept 3.756282e-141 rate_marriage 0.000000e+00 age 2.221918e-11 yrs_married 1.219200e-02 dtype: float64
代码单元 53
print(res_e2.summary())
print(res_e.summary())
原文输出(单元 53):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs_sum No. Observations: 130
Model: GLM Df Residuals: 127
Model Family: Poisson Df Model: 2
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -763.61
Date: Thu, 27 Aug 2026 Deviance: 1013.2
Time: 06:50:44 Pearson chi2: 895.
No. Iterations: 6 Pseudo R-squ. (CS): 1.000
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.0754 0.050 41.512 0.000 1.977 2.173
rate_marriage -0.4947 0.012 -41.743 0.000 -0.518 -0.471
yrs_married -0.0360 0.002 -17.542 0.000 -0.040 -0.032
=================================================================================
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs_sum No. Observations: 130
Model: GLM Df Residuals: 126
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -740.75
Date: Thu, 27 Aug 2026 Deviance: 967.46
Time: 06:50:44 Pearson chi2: 926.
No. Iterations: 6 Pseudo R-squ. (CS): 1.000
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
代码单元 54
print(res_f2.summary())
print(res_f.summary())
原文输出(单元 54):
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs No. Observations: 476
Model: GLM Df Residuals: 6363
Model Family: Poisson Df Model: 2
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -10374.
Date: Thu, 27 Aug 2026 Deviance: 15420.
Time: 06:50:44 Pearson chi2: 3.24e+04
No. Iterations: 6 Pseudo R-squ. (CS): 0.9729
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.0754 0.050 41.512 0.000 1.977 2.173
rate_marriage -0.4947 0.012 -41.743 0.000 -0.518 -0.471
yrs_married -0.0360 0.002 -17.542 0.000 -0.040 -0.032
=================================================================================
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: affairs No. Observations: 476
Model: GLM Df Residuals: 6362
Model Family: Poisson Df Model: 3
Link Function: Log Scale: 1.0000
Method: IRLS Log-Likelihood: -10351.
Date: Thu, 27 Aug 2026 Deviance: 15375.
Time: 06:50:44 Pearson chi2: 3.23e+04
No. Iterations: 6 Pseudo R-squ. (CS): 0.9754
Covariance Type: nonrobust
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.7155 0.107 25.294 0.000 2.505 2.926
rate_marriage -0.4952 0.012 -41.702 0.000 -0.518 -0.472
age -0.0299 0.004 -6.691 0.000 -0.039 -0.021
yrs_married -0.0108 0.004 -2.507 0.012 -0.019 -0.002
=================================================================================
来源、版本与许可
原文:Weighted Generalized Linear Models,statsmodels 0.15.0 文档。对应的 Show Source notebook 保留了 2026-08-27 的执行结果;项目中的 notebook 源文件 位于 statsmodels 仓库,采用仓库的 BSD 三条款许可。
本文汉化全部说明,保留 54 个代码单元和原始输出;校正了原文“467 条”的笔误,并补充了数据取整、权重语义、似然比倍数和内部属性的边界说明。输出中的日期和统计值均属于原 notebook。
文档版权:© 2009–2025 Josef Perktold、Skipper Seabold、Jonathan Taylor、statsmodels-developers。以下保留项目许可证的版权、条件与免责声明:
Copyright (C) 2006, Jonathan E. Taylor All rights reserved. Copyright (c) 2006-2008 Scipy Developers. All rights reserved. Copyright (c) 2009-2018 statsmodels Developers. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: a. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. b. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. c. Neither the name of statsmodels nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.











暂无评论内容