コントラストの概要¶
[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`個の水準を持つカテゴリ変数のコントラストのセットは、ダミー変数の合計からも独立している、因子水準平均の関数的に独立した`k-1`個の線形結合のセットです。ダミーコーディング自体は間違っていません。すべての係数を捉えますが、ANOVAのようにモデルが係数の独立性を仮定する場合、問題が複雑になります。線形回帰モデルは係数の独立性を仮定しないため、ダミーコーディングはこの文脈で教えられている唯一のコーディングであることがよくあります。
Patsyのコントラスト行列を見てみましょう。UCLA ATSのデータを使用します。まず、データをロードしましょう。
サンプルデータ¶
[2]:
import pandas as pd
url = "https://stats.idre.ucla.edu/stat/data/hsb2.csv"
hsb2 = pd.read_table(url, delimiter=",")
[3]:
hsb2.head(10)
[3]:
| ID | 女性 | 人種 | 社会経済的地位 | 学校の種類 | プログラム | 読解力 | 作文力 | 数学 | 科学 | 社会 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 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 = 白人））の各水準の平均を見ることは有益です。
[4]:
hsb2.groupby("race")["write"].mean()
[4]:
race
1 46.458333
2 58.000000
3 48.200000
4 54.055172
Name: write, dtype: float64
トリートメント（ダミー）コーディング¶
ダミーコーディングは、おそらく最もよく知られているコーディングスキームです。カテゴリ変数の各水準を基準水準と比較します。基準水準は切片の値です。順序付けられていないカテゴリ因子に対するPatsyのデフォルトのコントラストです。人種に対するトリートメントコントラスト行列は次のようになります。
[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`変数をどのようにエンコードするかを見てみましょう。
[6]:
hsb2.race.head(10)
[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
[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.]]
[8]:
pd.get_dummies(hsb2.race.values, drop_first=False)
[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行×4列
これは少しトリックです。`race`カテゴリは、ゼロベースのインデックスに都合よくマッピングされているためです。そうでない場合、この変換は内部で行われるため、これは一般的には機能しませんが、アイデアを修正するための便利な演習です.以下の３つのコントラストを使った出力を示します。
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 15:50:53 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.
人種のコントラストを明示的に指定しました。ただし、トリートメントがデフォルトであるため、これを省略することもできます。
単純コーディング¶
トリートメントコーディングと同様に、単純コーディングは各水準を固定基準水準と比較します。ただし、単純コーディングでは、切片はすべての因子の水準の全体平均です。Patsyには単純コントラストが含まれていませんが、独自のコントラストを簡単に定義できます。そのためには、patsy.contrast.ContrastMatrixインスタンスを返すcode_with_interceptメソッドとcode_without_interceptメソッドを含むクラスを記述します。
[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]))
[11]:
hsb2.groupby("race")["write"].mean().mean()
[11]:
np.float64(51.67837643678162)
[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]]
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 15:50:53 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は他のすべてと比較されます。
[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.]]
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 15:50:53 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.
これは、すべての係数の合計がゼロになるように強制するパラメーター化に対応します。ここでの切片は、各水準ごとの従属変数の平均の平均である全体平均です。
[16]:
hsb2.groupby("race")["write"].mean().mean()
[16]:
np.float64(51.67837643678162)
逆差分コーディング¶
逆差分コーディングでは、ある水準の従属変数の平均が、前の水準の従属変数の平均と比較されます。このタイプのコーディングは、名義変数または順序変数に役立つ場合があります。
[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]]
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 15:50:53 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の係数は、水準1の平均と比較した水準2の`write`の平均です。つまり、
[19]:
res.params["C(race, Diff)[D.1]"]
hsb2.groupby("race").mean()["write"][2] - hsb2.groupby("race").mean()["write"][1]
[19]:
np.float64(11.541666666666664)
ヘルマートコーディング¶
ここで使用するヘルマートコーディングは、逆ヘルマートコーディングと呼ばれることもあります。ある水準の従属変数の平均は、それ以前のすべての水準にわたる従属変数の平均と比較されます。そのため、順方向ヘルマートコーディングと区別するために、「逆」という名前が付けられることがあります。この比較は、人種のような名義変数にはあまり意味がありませんが、ヘルマートコントラストは次のように使用します。
[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.]]
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.78e-05
Time: 15:50:53 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つの水準の従属変数の平均です。
[22]:
grouped = hsb2.groupby("race")
grouped.mean()["write"].loc[4] - grouped.mean()["write"].loc[:3].mean()
[22]:
np.float64(3.169061302681982)
[23]:
grouped.mean()["write"]
[23]:
race
1 46.458333
2 58.000000
3 48.200000
4 54.055172
Name: write, dtype: float64
ご覧のとおり、これらは定数までしか等しくありません。ヘルマートコントラストの他のバージョンは、平均の実際の差を示します。いずれにせよ、仮説検定は同じです。
[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`から順序付きカテゴリを作成する必要があります。
[25]:
hsb2["readcat"] = np.asarray(pd.cut(hsb2.read, bins=4))
hsb2["readcat"] = hsb2["readcat"].astype(object)
hsb2.groupby("readcat").mean()["write"]
[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
[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 ]]
[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: Thu, 03 Oct 2024 Prob (F-statistic): 5.95e-18
Time: 15:50:53 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`に有意な線形効果がありますが、二次効果と三次効果は有意ではありません。
最終更新日：2024年10月3日