コントラストの概要¶
In [1]:
import numpy as np
import statsmodels.api as sm
このドキュメントは、UCLAの優れたリソース（http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm）を元にしています。
Kカテゴリまたはレベルのカテゴリカル変数は、通常、K-1個のダミー変数のシーケンスとして回帰に組み込まれます。これは、レベルの平均に関する線形仮説に相当します。つまり、これらの変数に対する各テスト統計量は、そのレベルの平均が基準カテゴリの平均と統計的に有意に異なるかどうかをテストすることになります。このダミー符号化は R 用語で「処置符号化」と呼ばれ、この慣例に従います。ただし、異なる符号化方法があり、それぞれ異なる線形仮説のセットに相当します。
実際、ダミー符号化は技術的にはコントラスト符号化ではありません。これは、ダミー変数が 1 に加算され、モデルの切片から機能的に独立していないためです。一方、k レベルのカテゴリカル変数に対するコントラストのセットは、ファクターのレベル平均の線形結合であり、ダミー変数の合計からも独立しています。ダミー符号化は本質的に間違っているわけではありません。それはすべての係数を捕えますが、モデルが係数の独立性を前提とする場合（例：ANOVA）には問題を複雑にします。線形回帰モデルは係数の独立性を前提としないため、ダミー符号化はこの文脈で教えられる唯一の符号化方法であることが多いです。
Patsyのコントラスト行列を確認するために、UCLA ATSのデータを使用します。まず、データをロードしましょう。
例題データ¶
In [2]:
import pandas as pd
url = "https://stats.idre.ucla.edu/stat/data/hsb2.csv"
hsb2 = pd.read_table(url, delimiter=",")
In [3]:
hsb2.head(10)
Out[3]:
| id | female | race | ses | schtyp | prog | read | write | math | science | socst | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 70 | 0 | 4 | 1 | 1 | 1 | 57 | 52 | 41 | 47 | 57 |
| 1 | 121 | 1 | 4 | 2 | 1 | 3 | 68 | 59 | 53 | 63 | 61 |
| 2 | 86 | 0 | 4 | 3 | 1 | 1 | 44 | 33 | 54 | 58 | 31 |
| 3 | 141 | 0 | 4 | 3 | 1 | 3 | 63 | 44 | 47 | 53 | 56 |
| 4 | 172 | 0 | 4 | 2 | 1 | 2 | 47 | 52 | 57 | 53 | 61 |
| 5 | 113 | 0 | 4 | 2 | 1 | 2 | 44 | 52 | 51 | 63 | 61 |
| 6 | 50 | 0 | 3 | 2 | 1 | 1 | 50 | 59 | 42 | 53 | 61 |
| 7 | 11 | 0 | 1 | 2 | 1 | 2 | 34 | 46 | 45 | 39 | 36 |
| 8 | 84 | 0 | 4 | 2 | 1 | 1 | 63 | 57 | 54 | 58 | 51 |
| 9 | 48 | 0 | 3 | 2 | 1 | 2 | 57 | 55 | 52 | 50 | 51 |
それぞれの人種のレベルごとに（1 = ヒスパニック、2 = アジア系、3 = アフリカ系アメリカ人、4 = 白人）、従属変数の平均を調べることは有益です。
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
処置（ダミー）符号化¶
ダミー符号化は最もよく知られた符号化方式の1つです。これはカテゴリカル変数の各レベルを基準となる参照レベルと比較します。基準となる参照レベルは切片の値です。これは順序のないカテゴリカル因子に対するPatsyのデフォルトのコントラストです。人種に対する処置コントラスト行列は次のようになります。
In [5]:
from patsy.contrasts import Treatment
levels = [1, 2, 3, 4]
contrast = Treatment(reference=0).code_without_intercept(levels)
print(contrast.matrix)
[[0. 0. 0.] [1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]
ここでは reference=0 を使用しましたが、これは最初のレベルであるヒスパニックが参照カテゴリであり、他のレベルの効果がこのカテゴリに対して測定されることを意味します。上記のように、列はゼロに合計されず、したがって切片とは独立していません。明示的に言うと、これが race 変数をどのようにエンコードするかを見てみましょう。
In [6]:
hsb2.race.head(10)
Out[6]:
0 4 1 4 2 4 3 4 4 4 5 4 6 3 7 1 8 4 9 3 Name: race, dtype: int64
In [7]:
print(contrast.matrix[hsb2.race - 1, :][:20])
[[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.]]
In [8]:
pd.get_dummies(hsb2.race.values, drop_first=False)
Out[8]:
| 1 | 2 | 3 | 4 | |
|---|---|---|---|---|
| 0 | False | False | False | True |
| 1 | False | False | False | True |
| 2 | False | False | False | True |
| 3 | False | False | False | True |
| 4 | False | False | False | True |
| ... | ... | ... | ... | ... |
| 195 | False | True | False | False |
| 196 | False | False | False | True |
| 197 | False | False | False | True |
| 198 | False | False | False | True |
| 199 | False | False | False | True |
200 rows × 4 columns
これは少しコツが必要ですが、race カテゴリは便利にゼロベースのインデックスにマッピングされます。もしそうでない場合、この変換は内部で行われるため、一般的にはこの方法は機能しませんが、それでもアイデアを定着させるためには有用な練習です。以下は、上記の3つのコントラストを使用した出力を示しています。
In [9]:
from statsmodels.formula.api import ols
mod = ols("write ~ C(race, Treatment)", 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: Wed, 28 Aug 2024 Prob (F-statistic): 5.78e-05
Time: 10:25:31 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）がデフォルトであるため、これを省略することもできました。
シンプル符号化¶
処置符号化と同様に、シンプル符号化は各レベルを固定の参照レベルと比較します。ただし、シンプル符号化では、切片が要因のすべてのレベルの全体平均となります。Patsy には単純コントラストが含まれていませんが、自分自身でコントラストを定義することは簡単です。そのためには、code_with_interceptとcode_without_interceptメソッドを含むクラスを作成し、patsy.contrast.ContrastMatrixインスタンスを返すようにします。
In [10]:
from patsy.contrasts import ContrastMatrix
def _name_levels(prefix, levels):
return ["[%s%s]" % (prefix, level) for level in levels]
class Simple(object):
def _simple_contrast(self, levels):
nlevels = len(levels)
contr = -1.0 / nlevels * np.ones((nlevels, nlevels - 1))
contr[1:][np.diag_indices(nlevels - 1)] = (nlevels - 1.0) / 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]))
In [11]:
hsb2.groupby("race")["write"].mean().mean()
Out[11]:
51.67837643678162
In [12]:
contrast = Simple().code_without_intercept(levels)
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 [13]:
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: Wed, 28 Aug 2024 Prob (F-statistic): 5.78e-05
Time: 10:25:31 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.
合計（偏差）符号化¶
合計符号化は、特定のレベルにおける従属変数の平均と、すべてのレベルにおける従属変数の全体平均を比較します。つまり、最初の k-1 レベルと k レベルの間でのコントラストを使用します。この例では、レベル1は他のすべてのレベルと比較され、レベル2は他のすべてのレベルと比較され、レベル3も他のすべてのレベルと比較されます。
In [14]:
from patsy.contrasts import Sum
contrast = Sum().code_without_intercept(levels)
print(contrast.matrix)
[[ 1. 0. 0.] [ 0. 1. 0.] [ 0. 0. 1.] [-1. -1. -1.]]
In [15]:
mod = ols("write ~ C(race, Sum)", 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: Wed, 28 Aug 2024 Prob (F-statistic): 5.78e-05
Time: 10:25:31 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 [16]:
hsb2.groupby("race")["write"].mean().mean()
Out[16]:
51.67837643678162
後方差分符号化¶
後方差分符号化では、あるレベルの従属変数の平均と、前のレベルの従属変数の平均が比較されます。このタイプの符号化は、名義尺度または順序尺度の変数に役立つ場合があります。
In [17]:
from patsy.contrasts import Diff
contrast = Diff().code_without_intercept(levels)
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 [18]:
mod = ols("write ~ C(race, Diff)", 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: Wed, 28 Aug 2024 Prob (F-statistic): 5.78e-05
Time: 10:25:31 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 [19]:
res.params["C(race, Diff)[D.1]"]
hsb2.groupby("race").mean()["write"][2] - hsb2.groupby("race").mean()["write"][1]
Out[19]:
11.541666666666664
ヘルマート符号化¶
私たちのヘルマート符号化のバージョンは、時々「逆ヘルマート符号化（Reverse Helmert Coding）」と呼ばれます。この方法では、あるレベルの従属変数の平均を、それ以前のすべてのレベルにおける従属変数の平均と比較します。したがって、前方ヘルマート符号化との違いを示すために「逆」という名前が付けられています。この比較は、人種のような名義尺度の変数にはあまり意味がありませんが、私たちはヘルマートコントラストを次のように使用します。
In [20]:
from patsy.contrasts import Helmert
contrast = Helmert().code_without_intercept(levels)
print(contrast.matrix)
[[-1. -1. -1.] [ 1. -1. -1.] [ 0. 2. -1.] [ 0. 0. 3.]]
In [21]:
mod = ols("write ~ C(race, Helmert)", 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: Wed, 28 Aug 2024 Prob (F-statistic): 5.78e-05
Time: 10:25:31 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 [22]:
grouped = hsb2.groupby("race")
grouped.mean()["write"].loc[4] - grouped.mean()["write"].loc[:3].mean()
Out[22]:
3.169061302681982
In [23]:
grouped.mean()["write"]
Out[23]:
race 1 46.458333 2 58.000000 3 48.200000 4 54.055172 Name: write, dtype: float64
ご覧のように、これらは定数を除いて等しいだけです。ヘルマートコントラストの他のバージョンでは、実際の平均の差が得られます。それでも、仮説検定は同じです。
In [24]:
k = 4
c4 = 1.0 / k * (grouped.mean()["write"].loc[k] - grouped.mean()["write"].loc[: k - 1].mean())
k = 3
c3 = 1.0 / k * (grouped.mean()["write"].loc[k] - grouped.mean()["write"].loc[: k - 1].mean())
print(c4, c3)
0.7922653256704955 -1.3430555555555561
直交多項式符号化¶
k=4 のレベルでの多項式符号化で得られる係数は、カテゴリ変数における線形、二次、三次のトレンドです。ここでのカテゴリ変数は、基になる等間隔の数値変数で表されると仮定されています。したがって、このタイプのエンコーディングは、等間隔の順序付きカテゴリ変数にのみ使用されます。一般的に、多項式コントラストは k-1 次の多項式を生成します。race は順序付き因子変数ではないので、例として read を使用します。まず、read から順序付きカテゴリを作成する必要があります。
In [25]:
hsb2["readcat"] = np.asarray(pd.cut(hsb2.read, bins=4))
hsb2["readcat"] = hsb2["readcat"].astype(object)
hsb2.groupby("readcat").mean()["write"]
Out[25]:
readcat (27.952, 40.0] 42.772727 (40.0, 52.0] 49.978495 (52.0, 64.0] 56.563636 (64.0, 76.0] 61.833333 Name: write, dtype: float64
In [26]:
from patsy.contrasts import Poly
levels = hsb2.readcat.unique()
contrast = Poly().code_without_intercept(levels)
print(contrast.matrix)
[[-0.67082039 0.5 -0.2236068 ] [-0.2236068 -0.5 0.67082039] [ 0.2236068 -0.5 -0.67082039] [ 0.67082039 0.5 0.2236068 ]]
In [27]:
mod = ols("write ~ C(readcat, Poly)", data=hsb2)
res = mod.fit()
print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: write R-squared: 0.346
Model: OLS Adj. R-squared: 0.336
Method: Least Squares F-statistic: 34.51
Date: Wed, 28 Aug 2024 Prob (F-statistic): 5.95e-18
Time: 10:25:31 Log-Likelihood: -690.69
No. Observations: 200 AIC: 1389.
Df Residuals: 196 BIC: 1403.
Df Model: 3
Covariance Type: nonrobust
==============================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------------------
Intercept 52.7870 0.634 83.268 0.000 51.537 54.037
C(readcat, Poly).Linear 14.2587 1.484 9.607 0.000 11.332 17.186
C(readcat, Poly).Quadratic -0.9680 1.268 -0.764 0.446 -3.468 1.532
C(readcat, Poly).Cubic -0.1554 1.006 -0.154 0.877 -2.140 1.829
==============================================================================
Omnibus: 4.467 Durbin-Watson: 1.768
Prob(Omnibus): 0.107 Jarque-Bera (JB): 4.289
Skew: -0.307 Prob(JB): 0.117
Kurtosis: 2.628 Cond. No. 3.01
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
ご覧の通り、readcat は従属変数 write に対して有意な線形効果を持っていますが、二次効果および三次効果は有意ではありません。