Patsy：カテゴリ変数のためのコントラストコーディングシステム¶
注記
このドキュメントは、UCLAの優れたリソースに基づいています。
K個のカテゴリまたはレベルを持つカテゴリ変数は、通常、回帰分析ではK-1個のダミー変数のシーケンスとして入力されます。これは、レベルの平均値に関する線形仮説に相当します。つまり、これらの変数の各検定統計量は、そのレベルの平均値が基準カテゴリの平均値と統計的に有意に異なるかどうかを検定することに相当します。このダミーコーディングは、RではTreatmentコーディングと呼ばれ、この慣例に従います。しかし、異なる線形仮説に相当する異なるコーディング方法があります。
実際、ダミーコーディングは技術的にはコントラストコーディングではありません。これは、ダミー変数が1に足し合わされ、モデルの切片と機能的に独立していないためです。一方、k個のレベルを持つカテゴリ変数のコントラストの集合とは、因子レベルの平均値のk-1個の機能的に独立した線形結合の集合であり、ダミー変数の合計からも独立しています。ダミーコーディング自体は間違ってはいません。すべての係数を捉えますが、ANOVAなど、係数の独立性を仮定するモデルでは問題を複雑にします。線形回帰モデルは係数の独立性を仮定しないため、この文脈ではダミーコーディングがしばしば唯一教えられているコーディングです。
Patsyのコントラスト行列を確認するために、UCLA ATSのデータを使用します。まずデータをロードしましょう。
例データ¶
In [1]: import pandas
In [2]: url = 'https://stats.idre.ucla.edu/stat/data/hsb2.csv'
In [3]: hsb2 = pandas.read_csv(url)
各人種のレベル（(1 = ヒスパニック、2 = アジア系、3 = アフリカ系アメリカ人、4 = 白人)）について、従属変数writeの平均を見ることは有益でしょう。
In [4]: hsb2.groupby('race')['write'].mean()
Out[4]:
race
1 46.458333
2 58.000000
3 48.200000
4 54.055172
Name: write, dtype: float64
Treatment（ダミー）コーディング¶
ダミーコーディングはおそらく最もよく知られているコーディングスキームです。カテゴリ変数の各レベルを基準レベルと比較します。基準レベルは切片の値です。これは、順序付けられていないカテゴリ因子のPatsyにおけるデフォルトのコントラストです。人種に対するTreatmentコントラスト行列は次のようになります。
In [5]: from patsy.contrasts import Treatment
In [6]: levels = [1,2,3,4]
In [7]: contrast = Treatment(reference=0).code_without_intercept(levels)
In [8]: print(contrast.matrix)
[[0. 0. 0.]
[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
ここではreference=0を使用しました。これは、最初のレベルであるヒスパニックが、他のレベルの効果を測定するための基準カテゴリであることを意味します。前述のように、列の合計はゼロにならず、したがって切片とは独立していません。明確にするために、これがどのようにrace変数をエンコードするかを見てみましょう。
In [9]: contrast.matrix[hsb2.race-1, :][:20]
Out[9]:
array([[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 1., 0.],
[0., 0., 0.],
[0., 0., 1.],
[0., 1., 0.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 1., 0.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.],
[0., 0., 1.]])
これはちょっとしたトリックで、raceカテゴリが都合よく0ベースのインデックスにマッピングされているためです。そうでない場合は、この変換が内部で行われるため、一般的には機能しませんが、それでもアイデアを固定するための便利な練習です。以下は、上記3つのコントラストを使用した出力を示しています。
In [10]: from statsmodels.formula.api import ols
In [11]: mod = ols("write ~ C(race, Treatment)", data=hsb2)
In [12]: res = mod.fit()
In [13]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
Intercept 46.4583 1.842 25.218 0.000 42.825 50.091
C(race, Treatment)[T.2] 11.5417 3.286 3.512 0.001 5.061 18.022
C(race, Treatment)[T.3] 1.7417 2.732 0.637 0.525 -3.647 7.131
C(race, Treatment)[T.4] 7.5968 1.989 3.820 0.000 3.675 11.519
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 8.25
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
私たちは人種のコントラストを明示的に与えましたが、Treatmentがデフォルトであるため、これを省略することもできました。
Simpleコーディング¶
Treatmentコーディングと同様に、Simpleコーディングは各レベルを固定された基準レベルと比較します。しかし、Simpleコーディングでは、切片は因子のすべてのレベルの全体平均です。Simpleコントラストの実装方法については、ユーザー定義コーディングを参照してください。
In [14]: contrast = Simple().code_without_intercept(levels)
In [15]: print(contrast.matrix)
[[-0.25 -0.25 -0.25]
[ 0.75 -0.25 -0.25]
[-0.25 0.75 -0.25]
[-0.25 -0.25 0.75]]
In [16]: mod = ols("write ~ C(race, Simple)", data=hsb2)
In [17]: res = mod.fit()
In [18]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
Intercept 51.6784 0.982 52.619 0.000 49.741 53.615
C(race, Simple)[Simp.1] 11.5417 3.286 3.512 0.001 5.061 18.022
C(race, Simple)[Simp.2] 1.7417 2.732 0.637 0.525 -3.647 7.131
C(race, Simple)[Simp.3] 7.5968 1.989 3.820 0.000 3.675 11.519
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 7.03
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Sum（偏差）コーディング¶
Sumコーディングは、特定のレベルの従属変数の平均を、すべてのレベル全体での従属変数の全体平均と比較します。つまり、最初のk-1個のレベルとレベルkの間のコントラストを使用します。この例では、レベル1は他のすべてと比較され、レベル2は他のすべてと比較され、レベル3は他のすべてと比較されます。
In [19]: from patsy.contrasts import Sum
In [20]: contrast = Sum().code_without_intercept(levels)
In [21]: print(contrast.matrix)
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]
[-1. -1. -1.]]
In [22]: mod = ols("write ~ C(race, Sum)", data=hsb2)
In [23]: res = mod.fit()
In [24]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
Intercept 51.6784 0.982 52.619 0.000 49.741 53.615
C(race, Sum)[S.1] -5.2200 1.631 -3.200 0.002 -8.437 -2.003
C(race, Sum)[S.2] 6.3216 2.160 2.926 0.004 2.061 10.582
C(race, Sum)[S.3] -3.4784 1.732 -2.008 0.046 -6.895 -0.062
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 6.72
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
これは、すべての係数の合計をゼロにする制約を課すパラメータ化に対応します。ここで、切片は全体平均であり、全体平均は各レベルによる従属変数の平均の平均です。
In [25]: hsb2.groupby('race')['write'].mean().mean()
Out[25]: np.float64(51.67837643678162)
Backward Differenceコーディング¶
Backward Differenceコーディングでは、レベルの従属変数の平均を、前のレベルの従属変数の平均と比較します。このタイプのコーディングは、名義変数または順序変数に役立つ場合があります。
In [26]: from patsy.contrasts import Diff
In [27]: contrast = Diff().code_without_intercept(levels)
In [28]: print(contrast.matrix)
[[-0.75 -0.5 -0.25]
[ 0.25 -0.5 -0.25]
[ 0.25 0.5 -0.25]
[ 0.25 0.5 0.75]]
In [29]: mod = ols("write ~ C(race, Diff)", data=hsb2)
In [30]: res = mod.fit()
In [31]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
======================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------------
Intercept 51.6784 0.982 52.619 0.000 49.741 53.615
C(race, Diff)[D.1] 11.5417 3.286 3.512 0.001 5.061 18.022
C(race, Diff)[D.2] -9.8000 3.388 -2.893 0.004 -16.481 -3.119
C(race, Diff)[D.3] 5.8552 2.153 2.720 0.007 1.610 10.101
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 8.30
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
たとえば、ここではレベル1の係数は、レベル2におけるwriteの平均とレベル1における平均の差です。つまり、
In [32]: res.params["C(race, Diff)[D.1]"]
Out[32]: np.float64(11.541666666666575)
In [33]: hsb2.groupby('race').mean()["write"].loc[2] - \
....: hsb2.groupby('race').mean()["write"].loc[1]
....:
Out[33]: np.float64(11.541666666666664)
Helmertコーディング¶
私たちのHelmertコーディングのバージョンは、Reverse Helmertコーディングと呼ばれることもあります。レベルの従属変数の平均を、それ以前のすべてのレベルの従属変数の平均と比較します。したがって、「逆」という名前が、前方Helmertコーディングと区別するために使用されることがあります。この比較は、人種のような名義変数ではあまり意味がありませんが、Helmertコントラストは次のように使用します。
In [34]: from patsy.contrasts import Helmert
In [35]: contrast = Helmert().code_without_intercept(levels)
In [36]: print(contrast.matrix)
[[-1. -1. -1.]
[ 1. -1. -1.]
[ 0. 2. -1.]
[ 0. 0. 3.]]
In [37]: mod = ols("write ~ C(race, Helmert)", data=hsb2)
In [38]: res = mod.fit()
In [39]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
=========================================================================================
coef std err t P>|t| [0.025 0.975]
-----------------------------------------------------------------------------------------
Intercept 51.6784 0.982 52.619 0.000 49.741 53.615
C(race, Helmert)[H.2] 5.7708 1.643 3.512 0.001 2.530 9.011
C(race, Helmert)[H.3] -1.3431 0.867 -1.548 0.123 -3.054 0.368
C(race, Helmert)[H.4] 0.7923 0.372 2.130 0.034 0.059 1.526
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 7.26
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
例として、レベル4の比較は、レベル4の平均から取った前の3レベルの従属変数の平均です。
In [40]: grouped = hsb2.groupby('race')
In [41]: grouped.mean()["write"].loc[4] - grouped.mean()["write"].loc[:3].mean()
Out[41]: np.float64(3.169061302681982)
ご覧のとおり、これらは定数までしか等しくありません。Helmertコントラストの他のバージョンは、平均の実際の差を与えます。いずれにしても、仮説検定は同じです。
In [42]: k = 4
In [43]: 1./k * (grouped.mean()["write"].loc[k] - grouped.mean()["write"].loc[:k-1].mean())
Out[43]: np.float64(0.7922653256704955)
In [44]: k = 3
In [45]: 1./k * (grouped.mean()["write"].loc[k] - grouped.mean()["write"].loc[:k-1].mean())
Out[45]: np.float64(-1.3430555555555561)
直交多項式コーディング¶
k=4レベルの多項式コーディングで取られる係数は、カテゴリ変数の線形、2次、3次の傾向です。ここでのカテゴリ変数は、基礎となる等間隔の数値変数で表されていると仮定されています。したがって、このタイプのエンコーディングは、等間隔の順序付きカテゴリ変数にのみ使用されます。一般に、多項式コントラストはk-1次の多項式を生成します。raceは順序付き因子変数ではないため、例としてreadを使用しましょう。まず、readから順序付きカテゴリを作成する必要があります。
In [46]: _, bins = np.histogram(hsb2.read, 3)
In [47]: try: # requires numpy main
....: readcat = np.digitize(hsb2.read, bins, True)
....: except:
....: readcat = np.digitize(hsb2.read, bins)
....:
In [48]: hsb2['readcat'] = readcat
In [49]: hsb2.groupby('readcat').mean()['write']
Out[49]:
readcat
0 46.000000
1 44.980392
2 53.356436
3 60.127660
Name: write, dtype: float64
In [50]: from patsy.contrasts import Poly
In [51]: levels = hsb2.readcat.unique().tolist()
In [52]: contrast = Poly().code_without_intercept(levels)
In [53]: print(contrast.matrix)
[[-0.6708 0.5 -0.2236]
[-0.2236 -0.5 0.6708]
[ 0.2236 -0.5 -0.6708]
[ 0.6708 0.5 0.2236]]
In [54]: mod = ols("write ~ C(readcat, Poly)", data=hsb2)
In [55]: res = mod.fit()
In [56]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.320
Model: OLS Adj. R-squared: 0.309
Method: Least Squares F-statistic: 30.73
Date: Thu, 03 Oct 2024 Prob (F-statistic): 2.51e-16
Time: 16:08:41 Log-Likelihood: -694.54
No. Observations: 200 AIC: 1397.
Df Residuals: 196 BIC: 1410.
Df Model: 3
Covariance Type: nonrobust
==============================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------------------
Intercept 51.1161 2.018 25.324 0.000 47.135 55.097
C(readcat, Poly).Linear 11.3501 5.348 2.122 0.035 0.803 21.897
C(readcat, Poly).Quadratic 3.8954 4.037 0.965 0.336 -4.066 11.857
C(readcat, Poly).Cubic -2.4598 1.998 -1.231 0.220 -6.400 1.480
==============================================================================
Omnibus: 9.741 Durbin-Watson: 1.699
Prob(Omnibus): 0.008 Jarque-Bera (JB): 10.263
Skew: -0.535 Prob(JB): 0.00591
Kurtosis: 2.703 Cond. No. 13.7
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
ご覧のとおり、readcatは従属変数writeに有意な線形効果がありますが、有意な2次効果または3次効果はありません。
ユーザー定義コーディング¶
独自のコーディングを使用する場合は、patsy.contrast.ContrastMatrixインスタンスを返すcode_with_interceptメソッドとcode_without_interceptメソッドを含むコーディングクラスを作成する必要があります。
In [57]: from patsy.contrasts import ContrastMatrix
....:
....: def _name_levels(prefix, levels):
....: return ["[%s%s]" % (prefix, level) for level in levels]
....:
In [58]: class Simple(object):
....: def _simple_contrast(self, levels):
....: nlevels = len(levels)
....: contr = -1./nlevels * np.ones((nlevels, nlevels-1))
....: contr[1:][np.diag_indices(nlevels-1)] = (nlevels-1.)/nlevels
....: return contr
....:
....: def code_with_intercept(self, levels):
....: contrast = np.column_stack((np.ones(len(levels)),
....: self._simple_contrast(levels)))
....: return ContrastMatrix(contrast, _name_levels("Simp.", levels))
....:
....: def code_without_intercept(self, levels):
....: contrast = self._simple_contrast(levels)
....: return ContrastMatrix(contrast, _name_levels("Simp.", levels[:-1]))
....:
File <tokenize>:13
def code_without_intercept(self, levels):
^
IndentationError: unindent does not match any outer indentation level
In [60]: mod = ols("write ~ C(race, Simple)", data=hsb2)
....: res = mod.fit()
....: print(res.summary())
....:
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.107
Model: OLS Adj. R-squared: 0.093
Method: Least Squares F-statistic: 7.833
Date: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 16:08:41 Log-Likelihood: -721.77
No. Observations: 200 AIC: 1452.
Df Residuals: 196 BIC: 1465.
Df Model: 3
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
Intercept 51.6784 0.982 52.619 0.000 49.741 53.615
C(race, Simple)[Simp.1] 11.5417 3.286 3.512 0.001 5.061 18.022
C(race, Simple)[Simp.2] 1.7417 2.732 0.637 0.525 -3.647 7.131
C(race, Simple)[Simp.3] 7.5968 1.989 3.820 0.000 3.675 11.519
==============================================================================
Omnibus: 10.487 Durbin-Watson: 1.779
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.031
Skew: -0.551 Prob(JB): 0.00402
Kurtosis: 2.670 Cond. No. 7.03
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
最終更新日：2024年10月3日