重み付き一般化線形モデル¶
[1]:
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
重み付きGLM：ポアソン応答データ¶
データの読み込み¶
この例では、少数の外生変数を使用して不倫率を予測するために、affairデータセットを使用します。
`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データフレームにロードします。
[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]:
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個の観測値があります。選択した変数のみを考慮すると、一意の観測値の数は少なくなります。以下では、2つの方法で観測値を組み合わせます。1つ目は、すべての変数の値が同一である観測値を組み合わせ、2つ目は、同じ説明変数を持つ観測値を組み合わせます。
一意の観測値を持つデータセット¶
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: Thu, 03 Oct 2024 Deviance: 15375.
Time: 15:45:07 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.078702313363214)
集約データ（頻度付きの一意の観測値）¶
同一の観測値を結合し、頻度ウェイトを使用して観測値の多重度を考慮に入れると、まったく同じ結果が得られます。すべての同一の観測値の集計ではなく、観測値に関する情報が必要な場合は、一部の結果属性が異なります。たとえば、残差は`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: Thu, 03 Oct 2024 Deviance: 15375.
Time: 15:45:07 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.078702313363196)
`freq_weights`の代わりに`var_weights`を使用して集約¶
次に、`var_weights`と`freq_weights`を比較します。内生変数が同一の観測値ではなく平均を反映している場合、`var_weights`を組み込むのが一般的です。それが（一般的に）同じ結果を生成する理論的な理由はありません。
これは同じ結果を生成しますが、`var_weights`は有効な観測値の数を変更しないため、`df_resid`は`freq_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: Thu, 03 Oct 2024 Deviance: 15375.
Time: 15:45:07 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]:
(np.float64(68.45488160512002), np.float64(5.078702313363196))
集計または平均化されたデータ（説明変数の一意の値）¶
これらの場合、説明変数の値が同じ観測値を組み合わせます。対応する応答変数は、合計または平均のいずれかです。
`exposure`の使用¶
従属変数がすべての結合された観測値の応答の合計である場合、ポアソン分布の仮定の下では、分布は同じままですが、1つの集約された観測値によって表される個人の数によって与えられる、変化する`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: Thu, 03 Oct 2024 Deviance: 967.46
Time: 15:45:07 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.35078910917951)
var_weightsの使用¶
従属変数のすべての結合された値の平均を使用することもできます。この場合、分散は、1つの結合された観測値によって反映される総エクスポージャーの逆数に関連します。
[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: Thu, 03 Oct 2024 Deviance: 967.46
Time: 15:45:07 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
=================================================================================
比較¶
上記の要約出力では、関連するWald推論を持つ`params`と`cov_params`がバージョン間で一致していることがわかりました。これを以下にまとめ、個々の結果属性をバージョン間で比較します。
パラメータ推定値`params`、パラメータの標準誤差`bse`、およびパラメータがゼロであるという検定の`pvalues`はすべて一致します。ただし、尤度と適合度統計量である`llf`、`deviance`、`pearson_chi2`は部分的にしか一致しません。具体的には、集計バージョンは元のデータを使用した結果と一致しません。
**警告**：`llf`、`deviance`、`pearson_chi2`の動作は、将来のバージョンで変更される可能性があります。
説明変数の一意の値に対する応答変数の合計と平均の両方に、適切な尤度解釈があります。ただし、この解釈はこれら3つの統計量には反映されていません。計算上、これは集計データが使用された場合の調整が不足していることが原因である可能性があります。ただし、理論的には、これらの場合、特に`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 (集計データの結果) | |
|---|---|---|---|---|
| 切片 | 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 (集計データの結果) | |
|---|---|---|---|---|
| 切片 | 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 (集計データの結果) | |
|---|---|---|---|---|
| 切片 | 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]:
(np.float64(52.91343161907935),
np.float64(45.726693322505525),
np.float64(-22.863346661251853))
[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.913431618650066),
np.float64(45.726693322507344),
np.float64(-22.863346661253672))
集計データ：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]:
(np.float64(-31.618527525103445),
np.float64(45.72669332250598),
np.float64(-22.863346661252763))
[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.61852752510788),
np.float64(45.72669332250621),
np.float64(-22.863346661252763))
ピアソンカイ二乗統計量の調査¶
まず、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]:
(np.float64(894.5809002315148),
np.float64(926.1994277566182),
np.float64(894.5809002315148),
np.float64(926.199427756618))
[31]:
res_e._results.resid_response.mean(), res_e.model.family.variance(res_e.mu)[
:5
], res_e.mu[:5]
[31]:
(np.float64(-2.815935518950797e-13),
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.1994277566182)
[33]:
res_e2._results.resid_response.mean(), res_e2.model.family.variance(res_e2.mu)[
:5
], res_e2.mu[:5]
[33]:
(np.float64(-2.361188168064333e-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.5809002315148)
[35]:
(res_e2._results.resid_response ** 2).sum(), (res_e._results.resid_response ** 2).sum()
[35]:
(np.float64(51204.85737832321), np.float64(47104.64779595939))
符号が正しくない理由の一つとして、異なる分母で割られた二次項を減算していることが考えられます。関連するケースでは、文献では共通の分母を使用することが推奨されています。完全モデルと縮小モデルで同じ分散の仮定を使用して、ピアソンカイ二乗統計量を比較できます。
この場合、すべてのバージョンで、縮小モデルと完全モデルの間でスケールされたピアソンカイ二乗の差分が同じになります。(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.43314175121899)
[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.43314175121923)
[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.43314175122017)
[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.43314175121962)
残りの部分¶
ノートブックの残りの部分は、追加のチェックのみを含んでおり、無視できます。
[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.664817945858474)
[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.5809002315161),
np.float64(926.199427756624),
np.float64(-42.34720713518717),
np.float64(-32.68238918932538))
[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.5809002315161)
[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.199427756624)
[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.1477584802967)
[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())
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, 03 Oct 2024 Deviance: 1013.2
Time: 15:45:08 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, 03 Oct 2024 Deviance: 967.46
Time: 15:45:08 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: Thu, 03 Oct 2024 Deviance: 15420.
Time: 15:45:08 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, 03 Oct 2024 Deviance: 15375.
Time: 15:45:08 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
=================================================================================
最終更新日：2024年10月3日