センタリング
太郎丸 博
2016年10月 1日
1 Cenetring in Conventional Regression Models: 通常の回帰分析でのセンタリング
マルチレベル分析でのセンタリングについて述べる前に、通常の最小二乗法 (OLS) での回帰分析におけるセンタリングを概説しておこう。
1.1 What is Centering: センタリングとは？
回帰分析におけるセンタリングとは、変数の値から何らかの値（ふつうはその変数の平均値）を引くことで、変数の分布を平行移動させることである。すなわち、\(X\)という変数の平均値を\(\bar X\)、\(X\) を平均値でセンタリングしたものを \(X_c\) とすると、 \[ X_c = X - \bar X \qquad (1) \] である。独立変数を平均値でセンタリングすると、計算が簡単になったり、交互作用を含むモデルの推定結果が解釈しやすくなる。以下に架空の例を示そう。
# 架空のデータを作る
n <- 10 # サンプル・サイズ
X <- 1 : n
Y <- -5 + 2 * X + rnorm(n, sd = 2) # rnorm(n) で n 個の正規分布する乱数を返す。sd= は標準偏差を指定
Xc <- X - mean(X) # X を平均でセンタリング
cbind(Y, X, Xc) #Y, X と Xc の対応を確認せよ## Y X Xc
## [1,] -1.1475736 1 -4.5
## [2,] 0.0279304 2 -3.5
## [3,] 5.0620017 3 -2.5
## [4,] 1.3020446 4 -1.5
## [5,] 5.8488709 5 -0.5
## [6,] 8.7211013 6 0.5
## [7,] 9.1699810 7 1.5
## [8,] 11.6380606 8 2.5
## [9,] 9.2680014 9 3.5
## [10,] 14.9436375 10 4.5# Y と X, Xc の散布図を作る
plot(Y ~ X, xlim = c(min(Xc), max(X)), xlab = "X or Xc") # xlim = 描画するx軸の最小値と最大値を指定
par(new = TRUE) # 重ねてプロットするための指定
plot(Y ~ Xc, xlim = c(min(Xc), max(X)), xlab = "", pch=2, col = "red")
par(new = FALSE) # 重ね書きを解除（やらなくていいがいちおう）
legend("topleft", pch = 1 : 2, col = 1 : 2, # 凡例の追加
legend = c("X と Y の散布図", "Xc 〃 "))
arrows(0, -2, 0, max(Y)) # Y軸を矢印で描画
lm1 <- lm(Y ~ X)
lm1c <- lm(Y ~ Xc)
abline(lm1) # 回帰直線の描画
abline(lm1c, col = "red", lty=2)
図1: X と Y の散布図と Xc（Xを平均値でセンタリングした変数）の散布図を重ねたもの
上の図の黒い丸がもともとの独立変数 X と従属変数 Y の散布図で、赤い三角が X を平均値でセンタリングして作った Xc という変数と Y の散布図である。X をセンタリングすると散布図がそのまま横に平行移動することがわかる。しかし、X と Y の関連の仕方は、Xc と Y の関連の仕方とまったく同じである。それゆえ、それぞれに関して回帰直線を OLS で推定すると、傾きは同じになり、切片（回帰直線が Y 軸と交わる点の座標）だけが異なることがわかる。社会学では回帰直線の傾きだけが問題であり切片はあまり重要でないので、独立変数を平均値でセンタリングしても実害ない場合が多い。またセンタリングした変数を使って推定した回帰直線から、もともとの回帰直線の切片を計算することもできる。具体的には、上の二つの回帰直線をそれぞれ \[ \hat Y = a + bX \qquad (2)\\ \hat Y = a_c + bX_c \qquad (3)\\ \] と表記する。(2) 式より、 \[ a = \hat Y - bX \qquad (4) \] である。(1) 式を (3) 式に代入すると \[ \hat Y = a_c + b(X - \bar X)\\ \hat Y = a_c + bX - b \bar X \\ \hat Y - b X = a_c - b \bar X \qquad (5) \] である。この (5) 式を (4) 式に代入すると、 \[ a = a_c - b \bar X \qquad (5) \] である。つまり、センタリングした変数を使った回帰直線の切片と傾き、もともとの X の平均値がわかれば、もともとの切片も計算できるというわけである。以下が計算例。
coef(lm1); coef(lm1c) # coef() で切片と傾きの推定値だけをベクトルとして返す## (Intercept) X
## -2.478519 1.629441## (Intercept) Xc
## 6.483406 1.629441# (5) 式での lm1 の切片の計算。一致していることを確認せよ
coef(lm1c)[1] - coef(lm1c)[2] * mean(X) # coef(lm1c)[1] で切片だけ返す## (Intercept)
## -2.4785191.2 Interaction Effects and Centering
独立変数をセンタリングせずに交互作用効果を推定しようとすると、独立変数の主効果の解釈が難しくなることがある。具体的には、もしも独立変数の最小値がゼロより大きいか、最大値がゼロよりも小さいならば（これを「レンジにゼロを含まない」と表現することにする）、主効果の解釈は非常に難しい。以下のような架空の例を考えてみよう。A 社と B 社における年齢と賃金の関係を図示すると図２のようになり、年齢をセンタリングせずに会社ダミーとの交互作用効果を回帰分析で推定すると、下のような結果になる。
# 架空のデータを作る。とりあえず理解できなくてもよい
n <- 10 # 1社あたりのサンプル・サイズ
data2 <- data.frame(firm = gl(2, n, labels = c("A社", "B社")
) # gl() は因子を作るための関数
) # data.frame(変数名 = ) でデータフレームを返す
data2$age <- rep(20:(19 + n), 2) # rep(A, B) で A を B回繰り返すベクトルを返す
set.seed(1986)
data2$wage <- 150 + 100 * (data2$firm == "A社") +
5 * (data2$age - 20) +
9 * (data2$firm == "A社") * (data2$age - 20) +
rnorm(n * 2, sd = 12)
data2$wage <- round(data2$wage, 0)
data2## firm age wage
## 1 A社 20 249
## 2 A社 21 267
## 3 A社 22 281
## 4 A社 23 280
## 5 A社 24 312
## 6 A社 25 312
## 7 A社 26 344
## 8 A社 27 339
## 9 A社 28 343
## 10 A社 29 362
## 11 B社 20 160
## 12 B社 21 153
## 13 B社 22 169
## 14 B社 23 173
## 15 B社 24 173
## 16 B社 25 168
## 17 B社 26 210
## 18 B社 27 201
## 19 B社 28 216
## 20 B社 29 184library(ggplot2)
g <- ggplot(data2, aes(x = age, y = wage, colour = firm))
g <- g + geom_point()
g
図2: 2社の年齢と賃金の架空データ
lm2.0 <- lm(wage ~ firm + age, data = data2)
lm2.1 <- update(lm2.0, ~. + firm:age)
library(texreg)
screenreg(list(lm2.0, lm2.1))##
## ====================================
## Model 1 Model 2
## ------------------------------------
## (Intercept) 89.59 ** 7.33
## (29.71) (31.52)
## firmB社 -128.20 *** 36.32
## (6.87) (44.58)
## age 8.95 *** 12.31 ***
## (1.20) (1.28)
## firmB社:age -6.72 **
## (1.81)
## ------------------------------------
## R^2 0.96 0.98
## Adj. R^2 0.95 0.97
## Num. obs. 20 20
## RMSE 15.37 11.61
## ====================================
## *** p < 0.001, ** p < 0.01, * p < 0.05図を見ると、A社のほうが平均賃金は高く、両社とも年齢に比例して賃金が上昇しているが、A社のほうが賃金の上昇が速く、20歳ぐらいの時には100万円程度であった両社の賃金格差は、29歳時には、150万円程度まで広がっている。回帰分析の結果を見ると、主効果のみの Model 1 では B社のダミー変数の係数は負の有意な値を示しており、グラフと整合する結果が得られている。ところが、交互作用効果を投入した Model 2 では、B社のダミー変数の係数は正の値を示しており、直感に反する結果である。確かに年齢の係数は正の値、B社ダミーと年齢の交互作用は負の値で、直感どおりの結果であるが、B社ダミーの主効果だけはおかしく感じる読者も多いだろう。なぜこうなるのか？
交互作用効果を含むモデルにおける主効果は、調整変数がゼロのときのその変数の傾きを指す。調整変数とは、例えば、以下のようなモデルを例に考えてみよう。 \[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{1i} X_{2i} + \epsilon_i. \qquad (6) \] このモデルでは \(X_{1i}\)と\(X_{2i}\) の交互作用が仮定されているので、\(X_{1i}\) にとっての調整変数は \(X_{2i}\) であり、逆に言えば \(X_{2i}\) にとっての調整変数は \(X_{1i}\) である。上の例で言えば、\(X_{1i}\) をB社ダミー、\(X_{2i}\) を年齢と考えればよい。つまり、交互作用効果を含むモデルとは、変数の傾きが調整変数の値によって異なると仮定するモデルであるが、そのようなモデルにおける主効果は、調整変数がゼロの時の変数の傾きのことなのである。以下、このことを数式で説明しよう。
もしも \(X_{1i} = 0\) ならば（つまりA社に限って年齢と賃金の関係を見ると）、(6) 式に \(X_{1i} = 0\) を代入して、 \[ Y_i = \beta_0 + \beta_1 \cdot 0 + \beta_2 X_{2i} + \beta_3 \cdot 0 \cdot X_{2i} + \epsilon_i \\ Y_i = \beta_0 + \beta_2 X_{2i} + \epsilon_i \qquad (7) \] である。つまり、\(\beta_2\) は調整変数（B社ダミー）がゼロの時の \(X_{2i}\) の傾きである。ちなみに調整変数であるB社ダミーが 1 のとき（つまりB社に関しては）、(6) 式に \(X_{1i} = 1\) を代入して、 \[ Y_i = \beta_0 + \beta_1 \cdot 1 + \beta_2 X_{2i} + \beta_3 \cdot 1 \cdot X_{2i} + \epsilon_i \\ Y_i = \beta_0 + \beta_1 + \beta_2 X_{2i} + \beta_3 X_{2i} + \epsilon_i \\ Y_i = (\beta_0 + \beta_1) + (\beta_2 + \beta_3) X_{2i} + \epsilon_i \qquad (8) \] となる。つまり、B社に関して回帰直線を引くと切片が \(\beta_0 + \beta_1\) で、傾きが \(\beta_2 + \beta_3\) になる。つまり調整変数のほうの主効果が切片の大きさの違いを示し、交互作用効果の係数が年齢の傾きの違いを示している。
同じようにして、(6) 式に \(X_{2i} = 0\) を代入して、 \[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 \cdot 0 + \beta_3 X_{1i} \cdot 0 + \epsilon_i \\ Y_i = \beta_0 + \beta_1 X_{1i} + \epsilon_i\qquad (9) \] である。つまり、\(\beta_1\) は調整変数（年齢）がゼロの時の \(X_{1i}\) の傾き（つまり2社の賃金の差）である。ゼロ歳時の賃金の差など考えても無意味であるが、主効果の係数とは調整変数がゼロの時の傾きを示しているということである。
交互作用効果に関する主効果の意味をもう少し詳しく考えるために、 上の例にもどって、Model 2 の推定結果を会社別の回帰直線で図示すると以下のようになる。
g <- ggplot(data2, aes(x = age, y = wage, colour = firm))
g <- g + geom_point()
g <- g + xlim(0, 30) + ylim(0, 400) # xlim(), ylim でプロットする範囲を指定
# (7) 式より、A社の切片はモデル全体の切片 coef(lm2.1)[1]、A社の年齢の傾きは年齢の主効果 coef(lm2.1)[3] なので
g <- g + geom_abline(intercept = coef(lm2.1)[1], # geom_abline() で図に直線を追加
slope = coef(lm2.1)[3],
colour = "red")
# (8) 式より同様に B 社の回帰直線を描画
g <- g + geom_abline(intercept = coef(lm2.1)[1] + coef(lm2.1)[2],
slope = coef(lm2.1)[3] + coef(lm2.1)[4],
colour = "blue")
g <- g + annotate("segment", x = 0, xend = 0, colour = "black", size = 1,
y = coef(lm2.1)[1], yend = coef(lm2.1)[1] + coef(lm2.1)[2],
arrow = arrow(ends = "both", length = unit(.2, "cm"))
)
g <- g + geom_text(x=2, y=25, label="B社ダミー\n の主効果", colour = "black")
g
図3: 会社別に回帰直線をゼロ歳時まで外挿
外挿（独立変数の最小値より低い値や最大値より大きい値まで回帰直線を伸ばし無理に結果を一般化すること）を行い、回帰直線をゼロ歳時まで延ばしていくと、2社の賃金は逆転し、A社よりも B社のほうが高くなる。このゼロ歳時の 2社の賃金の予測値の差が、B社ダミーの主効果の大きさである。
このような問題は、年齢を平均値でセンタリングしてやれば解消する。もともとの平均値（24.5歳）はセンタリング後にはゼロになるので、B社ダミーの主効果は、センタリング後には、24.5歳時の両社の賃金の差の予測値となる。センタリングして計算した結果は以下の通り。
data2$age.c <- data2$age - mean(data2$age)
lm2.2 <- lm(wage ~ firm * age.c, data = data2)
screenreg(list(lm2.1, lm2.2),
custom.coef.names = c(NA, NA, NA, NA, "age", "firmB社:age"))##
## ====================================
## Model 1 Model 2
## ------------------------------------
## (Intercept) 7.33 308.90 ***
## (31.52) (3.67)
## firmB社 36.32 -128.20 ***
## (44.58) (5.19)
## age 12.31 *** 12.31 ***
## (1.28) (1.28)
## firmB社:age -6.72 ** -6.72 **
## (1.81) (1.81)
## ------------------------------------
## R^2 0.98 0.98
## Adj. R^2 0.97 0.97
## Num. obs. 20 20
## RMSE 11.61 11.61
## ====================================
## *** p < 0.001, ** p < 0.01, * p < 0.05年齢をセンタリングしても年齢の主効果と交互作用効果は変わらないが、調整変数である B社ダミーと切片の値が変わっていることがわかる。24.5歳時には B社のほうが 128 だけ賃金が低いという結果であり、この数字のほうが、実際のデータの特質をわかりやすく伝えられるだろう。また、年齢の傾きの 12.31 は B 社ダミーがゼロのときの（つまり A社の）年齢の傾きになる。
ちなみに、センタリングの効用はもう一つあって、切片から意味のある情報を引き出すことができるようになる。切片とは、すべての独立変数がゼロの時の従属変数の予測値である。この場合、B社ダミーがゼロでセンタリングした年齢がゼロ（つまり、年齢が 24.5 歳）のときの賃金の予測値である。上のモデル 2 の切片は 308.9 と推定されているので、B社の 24.5 歳時の賃金は 308.9 万円と、このモデルからは推定できるということである。センタリングしていないと、切片はゼロ歳時の賃金の予測値になるので、解釈しても意味がないだろう。
センタリングは必ず平均値で行う必要はなく、最小値以上、最大値以下で何か意味のある数字であれば何でもよい。ダミー変数の場合、必ず最小値はゼロなのでセンタリングの必要は特にないが、してもよい。その場合、調整変数の主効果は全体の平均的な効果ということになる。
2 Centering in Random Slope Model: ランダム傾きモデルでのセンタリング
マルチレベル分析でも、交互作用効果を推定したり、ランダム傾きモデルを推定する場合、当該の変数がレンジにゼロを含んでいなければ、センタリングすべきである。レベル1の独立変数のグループ間効果とグループ内効果を識別する必要がある場合、やはりセンタリングが必要になるので、以下で概説する。
ランダム傾きモデルを推定する場合にどうしてセンタリングが必要か説明しよう。以下のようなランダム傾きモデルを考える。 \[ Y_{ij} = \beta_{0j} + \beta_{1j}X_{ij} + \epsilon_{ij} \\ \beta_{0j} = \gamma_{00} + \mu_{0j} \\ \beta_{0j} = \gamma_{10} + \mu_{1j} \\ \]
架空のデータで独立変数一つだけのランダム切片モデル (Model 1) とランダム傾きモデル (Model 2) を推定した結果が以下である。
head(data1)## id.firm id.worker firm.size educ wage educ.c size.c
## 1 a 1 10 11 242 -2.104 -24.0204
## 2 a 2 10 7 267 -6.104 -24.0204
## 3 a 3 10 11 449 -2.104 -24.0204
## 4 a 4 10 7 238 -6.104 -24.0204
## 5 a 5 10 13 345 -0.104 -24.0204
## 6 b 6 209 7 259 -6.104 -22.0304library(lme4)
r.intercept1 <- lmer(wage ~ educ + (1 | id.firm), data = data1)
r.slope1 <- lmer(wage ~ educ + (educ | id.firm), data = data1)
screenreg(list(r.intercept1, r.slope1))##
## =======================================================
## Model 1 Model 2
## -------------------------------------------------------
## (Intercept) 51.59 46.52
## (27.08) (26.85)
## educ 26.95 *** 26.56 ***
## (1.91) (2.47)
## -------------------------------------------------------
## AIC 1396.74 1380.71
## BIC 1408.05 1397.68
## Log Likelihood -694.37 -684.35
## Num. obs. 125 125
## Num. groups: id.firm 25 25
## Var: id.firm (Intercept) 2114.70 6590.24
## Var: Residual 3209.57 2646.89
## Var: id.firm educ 85.37
## Cov: id.firm (Intercept) educ -750.05
## =======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05注目してほしいのは、切片のランダム効果の分散 “Var: id.firm (Intercept)” と、残差（レベル1のランダム効果）の分散 “Var: Residual” である。Model 2 のほうが当てはまりがよいので、残差の分散は Model 2 のほうが小さくなる。しかし、切片のランダム効果の分散は、2114.7 から 6590.24 へと 3 倍以上に増えている。このような現象が起きるのは、切片の分散とは、\(X_{ij}=0\) 、この例では教育年数がゼロのときの会社別の賃金の平均値の分散を推定していることになるからである。このサンプルでは、教育年数の最小値は 9 なので、このような切片の分散を推定しても意味がない。そこで、教育年数をサンプル全体の平均値でセンタリングしてもう一度推定しなおしたのが、以下である。
data1$educ.c <- data1$educ - mean(data1$educ)
r.intercept1c <- lmer(wage ~ educ.c + (1 | id.firm), data = data1)
r.slope1c <- lmer(wage ~ educ.c + (educ.c | id.firm), data = data1)
screenreg(list(r.intercept1c, r.slope1c))##
## =========================================================
## Model 1 Model 2
## ---------------------------------------------------------
## (Intercept) 404.69 *** 394.61 ***
## (10.50) (9.41)
## educ.c 26.95 *** 26.56 ***
## (1.91) (2.47)
## ---------------------------------------------------------
## AIC 1396.74 1380.71
## BIC 1408.05 1397.68
## Log Likelihood -694.37 -684.35
## Num. obs. 125 125
## Num. groups: id.firm 25 25
## Var: id.firm (Intercept) 2114.70 1591.37
## Var: Residual 3209.57 2646.89
## Var: id.firm educ.c 85.36
## Cov: id.firm (Intercept) educ.c 368.57
## =========================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05Model 2 の切片のランダム効果の分散 “Var: id.firm (Intercept)” がセンタリングしなかった場合よりもずっと小さくなっている。これは平均的な教育年数の従業員の会社による賃金の違いとして解釈できるという点で、センタリングしなかった場合よりも利点がある。センタリングしていないとランダム傾きモデルにおいては切片のランダム効果の分散を解釈することがほとんど不可能であるのに対して、センタリングすればランダム切片モデルにおける切片のランダム効果の分散とほぼ同様に理解することが可能になる。このメリットは無視できないように思える。
ランダム切片モデルに関しては、交互作用効果を推定しないのならば、センタリングは必要ない。
3 Identifying Within and Between Group Effects: センタリングによるグループ間効果とグループ内効果の識別
ある変数の効果がグループ間とグループ内で異なる場合がある。以下のような架空のデータを見てみよう。
n <- 20 # 一社当たりの人数
M <- 7 # 会社数
data3 <- data.frame(firm = gl(M, n))
firm.n <- as.numeric(data3$firm)
set.seed(1986)
data3$edu <- 8 + firm.n + rnorm(M * n, sd = 2)
data3$edu <- round(data3$edu, 0)
data3$wage <- 100 + 3 * data3$edu + firm.n * 70 + rnorm(M * n, sd = 10)
data3$wage <- round(data3$wage, 0)
head(data3)## firm edu wage
## 1 1 9 190
## 2 1 10 196
## 3 1 10 219
## 4 1 7 179
## 5 1 10 194
## 6 1 8 209# 横軸に edu 、縦軸に wage をとって会社別に色分けしてプロット
# 会社別の回帰直線もつける
g <- ggplot(data3, aes(x = edu, y = wage, colour =firm))
g <- g + geom_point()
g <- g + geom_smooth(method = "lm", se = FALSE) # 会社別回帰直線の追加
g <- g + geom_smooth(mapping = aes(colour = NULL), method = "lm",
se = FALSE, color = "black")
# グループ平均をグループ平均に回帰させる
wage.firm <- tapply(data3$wage, data3$firm, mean)
edu.firm <- tapply(data3$edu, data3$firm, mean)
lm.firm <- lm(wage.firm ~ edu.firm)
g <- g + geom_point(aes(x = edu.firm, y = wage.firm, colour =NULL),
data = data.frame(wage.firm, edu.firm),
size = 3)
g <- g + geom_abline(intercept = coef(lm.firm)[1], slope = coef(lm.firm)[2],
linetype = 2)
g
図4: 全データで回帰直線を推定した場合（黒の実線）と、グループ別に回帰直線を推定した場合（黒以外の線）、賃金のグループ平均を教育年数のグループ平均に回帰させた場合（黒の点線）の傾きの違い （黒い大きな点は、各グループの賃金の平均値と教育年数の平均値を示す）
図4 の色のついた点は、個人の教育年数と賃金を示している。各点の色は会社ごとに異なる色をつけてある。黒い大きめの点は、各会社の教育年数の平均値と賃金の平均値を示している。会社別に OLS で回帰分析 (\(n=20\)) した結果が、7つの色付きの線であり、7つの黒い大きな点をデータとして、OLS で回帰直線を引いたものが、黒い破線である (\(n=7\)) 。これは各社の賃金の平均値を、教育年数の平均値に回帰させた結果である。さらに黒い実線は、すべての個人をデータとして、賃金を教育年数に回帰させた結果である (\(n=140\)) 。最初のグループごとの回帰直線の傾きをグループ内効果、二番目のグループ平均を使った回帰直線の傾きをグループ間効果、三番目のすべてのサンプルを使った回帰直線の傾きを総効果と呼ぶことにする。
総効果はグループ内効果とグループ間の中間にくることが知られている。グループ内効果とグループ間効果は一致することもあるし、図4 のように一致しないこともある。両者の相違が重要な意味を持つ場合がある。例えば、単純に教育年数＝従業員の労働生産性に比例して賃金が高まるのではなく、従業員の平均教育年数の高さが会社の信用を高め、それが会社の収益を高め、結局従業員の平均賃金も高める、という仮説が考えられる。もしもグループ内効果よりもグループ間効果のほうが大きければ、この例の仮説と整合的な結果である。このような仮説を検証したい場合、ランダム切片モデルに、各社の教育年数の平均値をグループレベルの独立変数として投入し、グループ平均センタリングした教育年数も個人レベルの独立変数として同時にモデルに投入すればよい。以下で具体的に説明しよう。
3.1 Group Mean Centering: グループ平均センタリング
一番、一般的なセンタリング法は、サンプルすべての平均をもとの変数の値から引くという方法で、全体平均センタリング (grand mean centering) と呼んでおく。 \(X_{ij}\) のサンプル全体の平均を \(\bar X_{\cdot \cdot}\) 、全体平均センタリングした変数を \(X_{ij}^{\mathrm{全}}\) とすると、 \[ X_{ij}^{\mathrm{全}} = X_{ij} - \bar X_{\cdot \cdot} \] である。 これはすでに交差レベル交互作用を推定する際に紹介した。もう一つのセンタリング法は、グループ平均センタリング (group mean centering) と呼ばれ、全体の平均ではなく、個々のグループの平均値をもとの変数の値から引く。すなわち、\(X_{ij}\) のグループ \(j\) の平均を \(\bar X_{\cdot j}\) 、グループ平均センタリングをと \(X_{ij}^{\mathrm{グ}}\) とすると、 \[ X_{ij}^{\mathrm{グ}} = X_{ij} - \bar X_{\cdot j} \qquad (2.1) \] である。
R でグループ平均センタリングした変数を計算するためには、まず各個人に対応するグループ平均を計算する必要がある。そのためには、
ave(グループ平均センタリングしたいベクトル, グループIDベクトル)
で、これらのベクトルと同じ長さのグループ平均を返す。以下が具体例。
data3$edu.m <- ave(data3$edu, data3$firm) # 各社の平均教育年数を個人に割り当てる
data3$edu.gc <- data3$edu - data3$edu.m # グループ平均センタリング
data3[1:4, ] # 一社目のデータを最初の4人分だけ見る## firm edu wage edu.m edu.gc
## 1 1 9 190 9.5 -0.5
## 2 1 10 196 9.5 0.5
## 3 1 10 219 9.5 0.5
## 4 1 7 179 9.5 -2.5data3[21:24, ] # 二社目のデータを最初の4人分だけ見る## firm edu wage edu.m edu.gc
## 21 2 9 278 10.3 -1.3
## 22 2 11 283 10.3 0.7
## 23 2 7 263 10.3 -3.3
## 24 2 9 280 10.3 -1.3教育年数が同じ 9年でも会社の平均によってグループ平均センタリングした値は違っている。
3.2 Estimating Within and Between Effects: グループ間／内効果の推定
これらのグループ平均とグループ平均偏差をランダム切片モデルに独立変数として投入すると、以下のようになる。 \[ \begin{array}{lll} Y_{ij} &= \ \beta_{0j} + \beta_{1} X_{ij}^{\mathrm{グ}} + R_{ij} & \\ \beta_{0j} &= \ \gamma_{00} + \gamma_{01} \bar X_{\cdot j} + U_j & \end{array} \] つまり、各社の平均的な賃金を、平均教育年数で予測し、各個人の賃金は、各社の平均よりもどの程度教育年数が乖離しているかで予測するというモデルになっている。グループ間効果は \(\gamma_{01}\) で、グループ内効果は \(\beta_1\)で表現される。以下が分析例である。
library(lmerTest)
r.itcpt1.0 <- lmer(wage ~ edu.gc + (1 | firm), data = data3) # 個人レベルの変数だけ投入したモデルもいちおう推定
r.itcpt1.1 <- lmer(wage ~ edu.gc + edu.m + (1 | firm), data = data3) # グループ内／間効果を推定したモデル
summary(r.itcpt1.1)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: wage ~ edu.gc + edu.m + (1 | firm)
## Data: data3
##
## REML criterion at convergence: 1060.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2802 -0.6467 -0.1460 0.7042 2.8858
##
## Random effects:
## Groups Name Variance Std.Dev.
## firm (Intercept) 257.4 16.04
## Residual 102.7 10.14
## Number of obs: 140, groups: firm, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -513.3729 39.5302 5.0000 -12.987 4.83e-05 ***
## edu.gc 3.4091 0.4356 132.0000 7.827 1.44e-12 ***
## edu.m 75.9430 3.1861 5.0000 23.835 2.42e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) edu.gc
## edu.gc 0.000
## edu.m -0.988 0.000fixed1.0 <- lm(wage ~ firm + edu, data = data3) # おまけで固定効果モデルを推定
screenreg(list(r.itcpt1.0, r.itcpt1.1, fixed1.0), omit.coef = "firm")##
## ============================================================
## Model 1 Model 2 Model 3
## ------------------------------------------------------------
## (Intercept) 417.47 *** -513.37 *** 165.31 ***
## (59.85) (39.53) (4.72)
## edu.gc 3.41 *** 3.41 ***
## (0.44) (0.44)
## edu.m 75.94 ***
## (3.19)
## edu 3.41 ***
## (0.44)
## ------------------------------------------------------------
## AIC 1101.07 1070.56
## BIC 1112.84 1085.27
## Log Likelihood -546.53 -530.28
## Num. obs. 140 140 140
## Num. groups: firm 7 7
## Var: firm (Intercept) 25072.94 257.40
## Var: Residual 102.74 102.74
## R^2 1.00
## Adj. R^2 1.00
## RMSE 10.14
## ============================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05このデータはグループ数が ７しかないので、ランダム効果の分散や p値の推定は過少になっている可能性がある（screenreg() は固定効果は正規分布すると仮定して p値を計算する）が、ここで興味があるのは固定効果の推定値なので、それらは気にしない。上の表の Model 2 の推定値がグループ内効果とグループ間効果を推定したものである。まったく異なる値になっていることがわかるだろう。グループ内効果とグループ間効果が一致しているかどうかを知るためのやや保守的だが簡単な方法は、両者の信頼区間が重なっているかどうかをチェックすることである。重なっていなければ、両者は有意に異なると考えられる。この例ではグループ数が少ないのでパラメトリック・ブートストラップで推定する。
confint(r.itcpt1.1, method="boot", parm = "beta_") # "parm = "beta_" で固定効果の信頼区間だけ表示## 2.5 % 97.5 %
## (Intercept) -591.278204 -428.78429
## edu.gc 2.608971 4.22142
## edu.m 69.039651 82.42817まったく異なる値だとわかるだろう。
グループ数が多く両者が等しいかどうか厳密に検定したいならば、Wald 検定するという方法もある。ただし、グループ数が少ないと p 値は不正確になる。Wald 検定には、aod パッケージの wald.test() 関数が使える。
wald.test(係数の分散共分散ベクトル, 係数ベクトル, L =cbind(0, 1, -1, 0…))
“L =” は行列で帰無仮説を表す引数であり、いろいろな帰無仮説を検定することができるのだが、二つの固定効果が等しいという帰無仮説を検定したい場合、cbind() で推定した固定効果と同じ長さの数値を指定し、等しいかどうかを検定したい係数のところで 1 と -1 を指定し、それ以外の係数のところではすべてゼロを指定する。1 と -1 はどちらにどちらを指定しても同じ結果が得られる。
library(aod)
wald.test(vcov(r.itcpt1.1), fixef(r.itcpt1.1), L = cbind(0, 1, -1))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 508.8, df = 1, P(> X2) = 0.03.3 Several Ways of Estimating Within Effects
グループ内効果とは、けっきょくグループによる違いを完全に統制したときの効果のことなので、固定効果モデルでも同じ推定値が得られる。これは上の表の Model 3 からもわかるだろう。また、グループ内効果だけが知りたいのならば、Model 1 のようにグループ平均センタリングした変数だけをランダム切片モデルに投入しても、Model 2 や Model 3 と同じ推定値が得られる。
4 Differences of Centering Methods
ちなみに、グループ内効果とグループ間効果の大きさが異なるかどうかを検定したい場合、個人レベルの変数はグループ平均センタリングではなく、全体平均センタリングした変数を使ったほうが、検定そのものは簡単である（あるいは主効果のみのランダム切片モデルならばセンタリングしていないもとの変数でもよい）。すなわち、 \[ \begin{array}{lll} Y_{ij} &= \ \beta_{0j} + \beta_{1} X_{ij} + R_{ij} & (2.2)\\ \beta_{0j} &= \ \gamma_{00} + \gamma_{01} \bar X_{\cdot j} + U_j & (2.3) \end{array} \] というモデルを推定すれば、\(\beta_{1}\) がグループ内効果、\(\beta_{1}+\gamma_{01}\) がグループ間効果を表すので、\(\gamma_{01}\) が統計的に有意ならばグループ間効果とグループ内効果は有意に異なると言える。 (2.2) 式に (2.3) 式を代入すると、 \[ \begin{array}{lll} Y_{ij} &= \ \gamma_{00} + \gamma_{01} \bar X_{\cdot j} + U_j + \beta_{1} X_{ij} + R_{ij} & \\ &= \ \gamma_{00} + \gamma_{01} \bar X_{\cdot j} + \beta_{1} X_{ij} + U_j + R_{ij} & (2.4) \end{array} \] (2.1) 式より、 \[ X_{ij} = X_{ij}^{\mathrm{グ}} + \bar X_{\cdot j} \] だからこれを (2.4) 式に代入して、 \[ \begin{array}{lll} Y_{ij} &= \ \gamma_{00} + \gamma_{01} \bar X_{\cdot j} + \beta_{1} (X_{ij}^{\mathrm{グ}} + \bar X_{\cdot j}) + U_j + R_{ij} & \\ &= \ \gamma_{00} + (\gamma_{01} + \beta_{1}) \bar X_{\cdot j} + \beta_{1} X_{ij}^{\mathrm{グ}} + U_j + R_{ij} & (2.5) \end{array} \] である。つまり、(2.5) 式のモデルにおいては、グループ間効果は \(\gamma_{01} + \beta_{1}\) で表され、グループ内効果は \(\beta_{1}\) で表される。それゆえ、\(\gamma_{01}\) はグループ内効果とグループ間効果の差を示す。\(\gamma_{01}\)がゼロなら両者に差はない。逆に言えば、\(\gamma_{01}\) が有意にゼロから離れていれば、グループ内効果とグループ間効果は異なると言える。以下は計算例。
r.itcpt1.2 <- lmer(wage ~ I(edu - mean(edu)) + edu.m + (1 | firm), data = data3)
screenreg(list(r.itcpt1.1, r.itcpt1.2), custom.model.names = c("Model 2", "Model 4"))##
## ===============================================
## Model 2 Model 4
## -----------------------------------------------
## (Intercept) -513.37 *** -471.59 ***
## (39.53) (39.89)
## edu.gc 3.41 ***
## (0.44)
## edu.m 75.94 *** 72.53 ***
## (3.19) (3.22)
## I(edu - mean(edu)) 3.41 ***
## (0.44)
## -----------------------------------------------
## AIC 1070.56 1070.56
## BIC 1085.27 1085.27
## Log Likelihood -530.28 -530.28
## Num. obs. 140 140
## Num. groups: firm 7 7
## Var: firm (Intercept) 257.40 257.40
## Var: Residual 102.74 102.74
## ===============================================
## *** p < 0.001, ** p < 0.01, * p < 0.05Model 2 はすでに推定してあったもので、グループ平均センタリングした教育年数 “edu.gc” を投入したモデルであり、Model 4 は全体平均センタリングした教育年数 “I(edu - mean(edu))” を投入してある。edu.gc と I(edu - mean(edu)) の係数がぴったり一致しているのがわかるだろう。“edu.m” は各社の平均教育年数だが、Model 2 ではこの係数がグループ間効果を表すのに対し、Model 4 ではグループ間／内効果の差を示す。\[75.94 - 3.41 = 72.53\] なのでぴったり一致していることがわかる。
4.1 Wchich Centering is Better? どの方法を用いるべきか
どのセンタリング法を用いても表現の仕方が違うだけなので、AICのような適合度指標は同じになる。それゆえ、問題はどうすれば係数を解釈しやすくなるのか、という点にある。これは結局、分析結果から何を言おうとしているのか、ということによって答えが変わってくる。上記のように、グループ間／内効果が異なることを主張したいならば、全体平均センタリングしてグループ平均の係数が有意かどうか検討するのがよかろうし、両者に差があるのは自明で、それぞれの効果の大きさがどれぐらいか示したいのならば、グループ平均センタリングしたほうがわかりやすかろう。
また、ランダム 傾き モデルでは、どちらのセンタリングでもあまり実害はないことが多いように思うが、もしもランダム傾きを仮定した個人レベルの変数（上の例では教育年数）の平均がグループ（上の例では会社）によって大きく異なるならば、グループ平均センタリングのほうが解釈しやすいかもしれない。なぜなら、グループ平均センタリングしたほうが、切片の分散の解釈が簡単だからである。グループ平均センタリングした場合、切片の分散は各グループの平均教育年数における予測賃金の分散と解釈できる。全体平均センタリングした場合、切片の分散は教育年数が全体平均に一致するとき（上の例では教育年数 = 12.26 のとき）の各グループの予測賃金の分散と解釈できる。どちらでもいいのだが、平均値より上の値の者ばかりのグループや逆に下の者ばかりのグループがあると、全体平均センタリングした場合の切片の分散は無意味であると考えられる。