差の差の分析の話というか Python の statsmodels.formula.api.glm は R 言語の glm と違って被説明変数が集約されているときはつかえないのだろうかという話だが、集約しなければいいので何も本質的ではないし、集約されていてもつかえるのに見落としているだけかもしれない。
参考文献
- hypothesis testing - How do I perform a statistical test for a Difference-in-Differences analysis? - Cross Validated（2022年4月25日参照）．
- 「diference in diference statistical test」のような検索語で検索すると出てきたので参考にした。1位の回答者が R 言語を、2位の回答者が Stata 言語を使用している。
- http://www.omori.e.u-tokyo.ac.jp/coretext/（2022年4月25日参照）．
- [1] の手法がある文献を他に探したがこの「差の差の分析 (第5章）」が同じであった。この文献は Stata 言語である。この文献での例題は「ごみ焼却施設の建設が住宅価格を下落させたか」であり、分析結果は「下落させている（有意）」である。
- Weighted Generalized Linear Models — statsmodels（2022年4月25日参照）．
問題
以下の問題を考える。
※ [1] の例と同じ数値だがシチュエーションを勝手に改変している。なお、改変した意味はこれといってない。
この世界にはえいけんという検定試験があり、春試験と秋試験がある。えいけんの予備校であるほげ予備校には東口校と西口校がある。ある年の春試験の後、ほげ予備校は秋試験に向けて東口校でだけ新指導方法を試行した。結果、春試験と秋試験の結果を並べると以下のようになった（ほげ予備校の生徒は全員えいけんを受けるものとする）。新指導方法に切り替えたことによる合格率向上への効果は有意だろうか。
| 春試験 | 秋試験 | 合格率の差 | |||||
|---|---|---|---|---|---|---|---|
| 生徒数 | 合格者数 | 合格率 | 生徒数 | 合格者数 | 合格率 | ||
| 東口校 | 130 | 64 | 0.49 | 118 | 82 | 0.69 | 0.2 |
| 西口校 | 100 | 44 | 0.44 | 100 | 60 | 0.6 | 0.16 |
| 合格率の差の差 | 0.04 | ||||||
「新指導方法に切り替えたことによる合格率向上への効果」は、新指導方法に切り替えた東口校の「秋 - 春」差分から、西口校の「秋 - 春」差分を差し引けばよく、0.04 という正の値である。いま、この値が有意であるかに興味がある。以下、[1] の1位の回答と同様に、ロジスティック回帰の交互作用項の係数を検定する。
なお、「『新指導方法を受けたいから東口校を選ぶ』などということはできない」「『西口校のほうが秋に優秀な生徒が集まる』などということはない」などは仮定されているものとする。
回答
R による回答例
[1] のベストアンサーの方が示しているコードと同じである。ただ [1] のベストアンサーの方は数字を1箇所タイポしているので係数の推定値は [1] と揃っていない。何にせよ、交互作用項の係数 Branch1:Term1（東口校かつ秋試験）は正の値 0.2073 ではあるが、「交互作用項の係数がゼロである」の棄却域にはない。よって新指導方法の効果は有意であるとはいえない。N <- c(130, 118, 100, 100) NPassed <- c(64, 82, 44, 60) Branch <- factor(c("1", "1", "0", "0")) # 東口校:0, 西口校:1 Term <- factor(c("0", "1", "0", "1")) # 春試験:0, 秋試験:1 PassRate <- NPassed / N fit <- glm(PassRate ~ Branch * Term, family=binomial, weights=N) print(summary(fit))
Call:
glm(formula = PassRate ~ Branch * Term, family = binomial, weights = N)
Deviance Residuals:
[1] 0 0 0 0
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2412 0.2015 -1.197 0.2313
Branch1 0.2104 0.2671 0.788 0.4309
Term1 0.6466 0.2868 2.255 0.0242 *
Branch1:Term1 0.2073 0.3912 0.530 0.5961
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.7868e+01 on 3 degrees of freedom
Residual deviance: -1.7764e-15 on 0 degrees of freedom
AIC: 28.454
Number of Fisher Scoring iterations: 3
Python による失敗例
同じことを Python でやろうとして以下のようにかくとエラーになる。メッセージから察するに、「東口校の生徒 130 人のうち 49% が春試験に合格した」と解釈してくれず、「東口校の生徒 130 人は全員スコア 0.49 を取ったのだ」という扱いになっている。なので「完全に分離できてしまいます」などといってくるがそうではないのである。import pandas as pd import statsmodels.formula.api as smf import statsmodels.api as sm df = pd.DataFrame({ 'n': [130, 118, 100, 100], 'n_passed': [64, 82, 44, 60], 'branch': ['1', '1', '0', '0'], # 東口校:0, 西口校:1 'term': ['0', '1', '0', '1'], # 春試験:0, 秋試験:1 }) df['pass_rate'] = df['n_passed'] / df['n'] family = sm.families.Binomial() result = smf.glm( formula='pass_rate ~ branch * term', data=df, family=sm.families.Binomial(), var_weights=np.asarray(df['n']), # この状況で var_weights をつかうのは誤りである ).fit() print(result.summary())
PerfectSeparationError: Perfect separation detected, results not available
Python による愚直な成功例
「東口校の生徒の中に春試験に合格した子も合格しなかった子もいるんやぞ」ということを否応なくわからせるために、各行が生徒1人に対応するように愚直に整形した例が以下である。R での実行例と全く同じ結果が返ってきて安心はする。無論、このような実装は推奨できない。branch = ['1'] * 130 + ['1'] * 118 + ['0'] * 100 + ['0'] * 100 # 東口校:0, 西口校:1 term = ['0'] * 130 + ['1'] * 118 + ['0'] * 100 + ['1'] * 100 # 春試験:0, 秋試験:1 passed = [(1 if (i < 64) else 0) for i in range(130)] + \ [(1 if (i < 82) else 0) for i in range(118)] + \ [(1 if (i < 44) else 0) for i in range(100)] + \ [(1 if (i < 60) else 0) for i in range(100)] df = pd.DataFrame({'branch': branch, 'term': term, 'passed': passed}) print(df.head(3)) print(df.tail(3)) result = smf.glm(formula='passed ~ branch * term', data=df, family=sm.families.Binomial()).fit() print(result.summary())
branch term passed
0 1 0 1
1 1 0 1
2 1 0 1
branch term passed
445 0 1 0
446 0 1 0
447 0 1 0
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: passed No. Observations: 448
Model: GLM Df Residuals: 444
Model Family: Binomial Df Model: 3
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -298.57
Date: Mon, 25 Apr 2022 Deviance: 597.14
Time: 22:03:23 Pearson chi2: 448.
No. Iterations: 4 Pseudo R-squ. (CS): 0.03910
Covariance Type: nonrobust
=========================================================================================
coef std err z P>|z| [0.025 0.975]
-----------------------------------------------------------------------------------------
Intercept -0.2412 0.201 -1.197 0.231 -0.636 0.154
branch[T.1] 0.2104 0.267 0.788 0.431 -0.313 0.734
term[T.1] 0.6466 0.287 2.255 0.024 0.085 1.209
branch[T.1]:term[T.1] 0.2073 0.391 0.530 0.596 -0.559 0.974
=========================================================================================
Python による改善した成功例
先ほどはつい各行を生徒1人にしてしまったが、「合格した子たち」「合格しなかった子たち」はそれぞれ同一視してもよさそうである。というのが以下の例である。R と同じ結果になる。df = pd.DataFrame({
'n': [64, 130 - 64, 82, 118 - 82, 44, 100 - 44, 60, 100 - 60],
'passed': [1, 0, 1, 0, 1, 0, 1, 0],
'branch': ['1', '1', '1', '1', '0', '0', '0', '0'], # 東口校:0, 西口校:1
'term': ['0', '0', '1', '1', '0', '0', '1', '1'], # 春試験:0, 秋試験:1
})
family = sm.families.Binomial()
result = smf.glm(
formula='passed ~ branch * term',
data=df,
family=sm.families.Binomial(),
freq_weights=np.asarray(df['n']), # 同一の観測をまとめているとき freq_weights のはずである
).fit()
print(result.summary())
Generalized Linear Model Regression Results ============================================================================== Dep. Variable: passed No. Observations: 8 Model: GLM Df Residuals: 444 Model Family: Binomial Df Model: 3 Link Function: Logit Scale: 1.0000 Method: IRLS Log-Likelihood: -298.57 Date: Mon, 25 Apr 2022 Deviance: 597.14 Time: 22:26:35 Pearson chi2: 448. No. Iterations: 4 Pseudo R-squ. (CS): 0.8929 Covariance Type: nonrobust ========================================================================================= coef std err z P>|z| [0.025 0.975] ----------------------------------------------------------------------------------------- Intercept -0.2412 0.201 -1.197 0.231 -0.636 0.154 branch[T.1] 0.2104 0.267 0.788 0.431 -0.313 0.734 term[T.1] 0.6466 0.287 2.255 0.024 0.085 1.209 branch[T.1]:term[T.1] 0.2073 0.391 0.530 0.596 -0.559 0.974 =========================================================================================
