加重一般化線形モデル¶
[1]:
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
重み付きGLM: ポアソン応答データ¶
データの読み込み¶
この例では、いくつかの外生変数を用いて不倫率を予測するために、アフェアデータセットを使用します。
重みは、freq_weightsがデータのレコードを繰り返すことに相当することを示すために生成されます。一方で、var_weightsはデータを集約することに相当します。
[2]:
print(sm.datasets.fair.NOTE)
::
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 relgious, 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-colloar; 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 |
以下では、主にポアソン分布を使用します。小数点の計算を使用することも可能ですが、カウント分布を持つために整数に変換します。
[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]:
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], dtype=int64)
観測の圧縮と集約¶
元のデータセットには 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.reset_index(inplace=True)
dc.rename(columns={"const": "freq"}, inplace=True)
print(dc.shape)
dc.head()
(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 |
一意の説明変数（exog）を持つデータセット¶
次のデータセットでは、説明変数の値が同じである観測値を組み合わせます。ただし、応答変数は組み合わせた観測値の間で異なる可能性があるため、すべての組み合わせた観測値について、応答変数の平均値と合計を計算します。
再度、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()
(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 |
観測値を組み合わせた後、467 のユニークな観測値を持つデータフレーム dc と、説明変数のユニークな値を持つ 130 の観測値を持つデータフレーム df_a が得られました。
[11]:
print("number of rows: \noriginal, with unique observations, with unique exog")
data.shape[0], dc.shape[0], df_a.shape[0]
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())
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: Wed, 28 Aug 2024 Deviance: 15375.
Time: 10:36:12 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]:
5.078702313363182
圧縮データ（頻度を持つユニークな観測値）¶
同一の観測値を統合し、観測値の重複を考慮するために頻度重みを使用することで、全く同じ結果が得られます。ただし、観測値の情報を求めている場合、または全ての同一観測値の集約について情報を求めている場合では、いくつかの結果属性が異なることがあります。例えば、残差は「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())
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: Wed, 28 Aug 2024 Deviance: 15375.
Time: 10:36:12 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]:
5.078702313363233
freq_weightsの代わりにvar_weightsを使用して圧縮¶
次に、var_weightsとfreq_weightsを比較します。内生変数が平均を反映し、同一の観測値ではない場合に、var_weightsを組み込むことは一般的な方法です。 なぜそれが同じ結果を生むのか（一般的に）は理論的な理由が見当たりません。
これにより同じ結果が得られますが、df_residはfreq_weightsの例と異なります。これは、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())
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: Wed, 28 Aug 2024 Deviance: 15375.
Time: 10:36:12 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のため不正確です。元の
df_residを使用すれば正しい結果が得られます。[17]:
res_fv.pearson_chi2 / res_fv.df_resid, res_f.pearson_chi2 / res_f.df_resid
[17]:
(68.45488160512053, 5.078702313363233)
集計または平均化されたデータ（説明変数のユニークな値）¶
これらのケースでは、説明変数の値が同じ観測値を組み合わせます。対応する従属変数は、合計または平均となります。
exposureの使用¶
もし従属変数が結合されたすべての観測値の反応の合計である場合、ポアソン分布の仮定の下では分布は変わりませんが、集計された観測値が表す個体数によって異なるexposure（曝露）が与えられます。
パラメータの推定値やパラメータの共分散は元のデータと同じですが、対数尤度、偏差、ピアソンのカイ二乗は異なります。
[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())
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: Wed, 28 Aug 2024 Deviance: 967.46
Time: 10:36:12 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]:
7.350789109179566
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())
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: Wed, 28 Aug 2024 Deviance: 967.46
Time: 10:36:12 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 に関連するウォルド推定は、バージョン間で一致しています。これを、バージョンごとの個別の結果属性を比較することで要約します。
パラメータ推定値 params、パラメータの標準誤差 bse、およびパラメータがゼロであるかどうかを検定するための pvalues はすべて一致します。ただし、尤度および適合度統計量である llf、deviance、pearson_chi2 は部分的にしか一致しません。具体的には、集約されたバージョンは元のデータを使用した結果と一致しません。
警告: llf、deviance、および pearson_chi2 の挙動は、今後のバージョンで変更される可能性があります。
説明変数の固有の値に対する応答変数の合計および平均は、適切な尤度解釈を持っています。しかし、この解釈はこれらの三つの統計量には反映されていません。計算的には、集約データを使用した場合に調整が欠落していることが原因かもしれません。しかし、理論的には、特に尤度分析が不適切で結果を準尤度推定として解釈すべき場合に、誤指定されたケースの var_weights に関しては、これを考慮することができます。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 |
尤度比型検定¶
上記で、尤度と関連統計量が集約データと元の個別データで一致しないことを見ました。次に、尤度比検定と逸脱度の差がバージョン間で一致することを示しますが、ピアソンのカイ二乗検定は一致しません。
以前と同様に：まだ十分に明確ではなく、変更される可能性があります。
テストケースとして、age変数を削除し、縮小または制約付きモデルと完全または非制約モデルの間の差として尤度比型統計量を計算します。
元の観測値と頻度重み¶
[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]:
(52.91343161885379, 45.726693322507344, -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]:
(52.91343161840632, 45.726693322507344, -22.863346661251853)
集約データ: exposure と var_weights¶
注: LR検定は元の観測結果と一致しますが、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]:
(-31.618527525110608, 45.72669332250632, -22.863346661252876)
[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]:
(-31.618527525115155, 45.72669332250621, -22.863346661251853)
ピアソンのカイ二乗統計量の調査¶
まず、pearson_chi2 と resid_pearson の計算に基本的なバグがないか、いくつかの簡単なチェックを行います。
[30]:
res_e2.pearson_chi2, res_e.pearson_chi2, (res_e2.resid_pearson ** 2).sum(), (
res_e.resid_pearson ** 2
).sum()
[30]:
(894.5809002315146, 926.1994277566253, 894.5809002315146, 926.1994277566251)
[31]:
res_e._results.resid_response.mean(), res_e.model.family.variance(res_e.mu)[
:5
], res_e.mu[:5]
[31]:
(6.864935970112968e-14,
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]:
926.1994277566253
[33]:
res_e2._results.resid_response.mean(), res_e2.model.family.variance(res_e2.mu)[
:5
], res_e2.mu[:5]
[33]:
(-2.7109938225923824e-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]:
894.5809002315146
[35]:
(res_e2._results.resid_response ** 2).sum(), (res_e._results.resid_response ** 2).sum()
[35]:
(51204.85737832321, 47104.64779595963)
誤った符号の可能な理由の一つは、異なる分母で割られた2次項を引いていることです。関連するいくつかのケースでは、文献で共通の分母を使用することが推奨されています。フルモデルと縮小モデルで同じ分散の仮定を使用して、ピアソンのカイ二乗統計量を比較することができます。
この場合、縮小モデルとフルモデルの間で、すべてのバージョンで同じピアソンカイ二乗スケールの差を得ることができます。（問題#3616は、これをさらに追跡することを目的としています。）
[36]:
(
(res_e2._results.resid_response ** 2 - res_e._results.resid_response ** 2)
/ res_e2.model.family.variance(res_e2.mu)
).sum()
[36]:
44.433141751219104
[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]:
44.43314175121908
[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]:
44.43314175121938
[39]:
(
(res_o2._results.resid_response ** 2 - res_o._results.resid_response ** 2)
/ res_o2.model.family.variance(res_o2.mu)
).sum()
[39]:
44.43314175121891
残り¶
ノートブックの残りの部分は、いくつかの追加的な確認を含んでおり、無視しても構いません。
[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], dtype=int64))
[41]:
res_e2.resid_pearson.sum() - res_e.resid_pearson.sum()
[41]:
-9.664817945864009
[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]:
(894.580900231515, 926.1994277566301, -42.347207135188135, -32.68238918932136)
[44]:
(
(res_a2._results.resid_response ** 2)
/ res_a2.model.family.variance(res_a2.mu)
* res_a2.model.var_weights
).sum()
[44]:
894.580900231515
[45]:
(
(res_a._results.resid_response ** 2)
/ res_a.model.family.variance(res_a.mu)
* res_a.model.var_weights
).sum()
[45]:
926.1994277566301
[46]:
(
(res_a._results.resid_response ** 2)
/ res_a.model.family.variance(res_a2.mu)
* res_a.model.var_weights
).sum()
[46]:
850.1477584802958
[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]:
(1.8390448369994542e-07, 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())
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: Wed, 28 Aug 2024 Deviance: 1013.2
Time: 10:36:13 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: Wed, 28 Aug 2024 Deviance: 967.46
Time: 10:36:13 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())
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: Wed, 28 Aug 2024 Deviance: 15420.
Time: 10:36:13 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: Wed, 28 Aug 2024 Deviance: 15375.
Time: 10:36:13 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
=================================================================================
最終更新日:
2025年01月28日