【備忘】交互作用のparameterizationについて考えたこと

実務上の分析をするなかでふと思ったこと。
いちいち書き留めるほどのことでもないけど、そうしないと忘れちゃうので（老兵はつらい）。

何が問題か
じゃあどうすればいいか
- Zが離散変数：ダミー変数のcodingを工夫
- Zが連続変数の場合はどう考えるとよいか
Conclusion

何が問題か

たとえば（例はほんとになんでもいいんだけど）
テストの成績 $Y$ に対しての通塾 $X=0,1$ の効果が、朝ごはんを食べているかどうか $Z=0,1$ によって変わるか否か*1、ということに関して調べたいとする。

そういう関心をもっているとき、我々は往々にして以下のようなモデルを推定する。

$Y_i = \alpha X_i + \beta Z_i+ \gamma X_i Z_i +\epsilon_i$

このとき、 $\gamma$ の係数推定値（やその統計的有意性）をみて、前述の問いを検証することが多い。
$\gamma$ が正でかつ有意であれば、「朝ごはんを食べることにより、通塾の成績への効果が増幅される」といった風に解釈するわけだ。

ここまではまぁとても基礎的なハナシなんだけど、実はひとつ盲点を抱えている
それは $X$ の効果が $Z$ によって高められるのは、通塾勢の成績が上がることによってなのか、非通塾勢の成績が下がることによってなのかが、 $\gamma$ をみただけではわからない、ということである。
ことばを替えると、天井を引き上げる効果なのか床を押し下げる効果なのかがわからない。

図で表すと、以下の場合両方とも $X$ の効果が $Z=1$ のときに20増えているけれど、上の図は天井への効果によるもので、下の図は床への効果によるものである。

f:id:ronri_rukeichi:20201027220535p:plain — 交互作用1（天井への効果）

f:id:ronri_rukeichi:20201027220711p:plain — 交互作用2（床への効果）

しかし上の定式化であると、一様に $\gamma = 20$ ということになってしまう。
パラメタライゼーションを工夫することにより、天井への効果と床への効果を識別して推定する術を探す。

じゃあどうすればいいか

まずZが離散変数（というかダミー変数）の時を考えて、次に

Zが離散変数：ダミー変数のcodingを工夫

以下のような仮想データを生成する。天井への効果が+15,底への効果が-10で、 $\gamma=25$ となるような場合を想定している。

x0z0 <- rnorm(100 ,  30 ,10)
x0z1 <- rnorm( 100 , 20 ,10)
x1z0 <- rnorm(100 , 50 , 10)
x1z1 <- rnorm(100 , 65 , 10)
dta1 <- data.frame(  Y = c( x0z0, x0z1 , x1z0 , x1z1) , X = c(rep(0, 200), rep(1 , 200 )), Z = c(rep(0, 100), rep(1 , 100 ),rep(0, 100), rep(1 , 100 ) ))
lm_res <- lm( Y ~ X *Z , data = dta1)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  30.0541     0.9335  32.194  < 2e-16 ***
#   X            20.0450     1.3202  15.183  < 2e-16 ***
#   Z            -9.5118     1.3202  -7.205 2.96e-12 ***
#   X:Z        24.3167     1.8671  13.024  < 2e-16 ***

$\gamma = 24.3167$ である
われわれが目指すのは、 $\gamma_1\fallingdotseq15 , \gamma_2 \fallingdotseq -10$ でかつ $\gamma_1 - \gamma_2 = \gamma$ という推定値を与えるものだ

Xが2値 × Zが2値＝4通りがあるので、「Z=0のときのXの効果（ $\beta_0$ ）」のほかに
「Zの天井への効果(X=1のときの効果)」「Zの天井への効果(X=0のときの効果)」という三つを識別すると、自由度を消費しきってしまう。
なので、上述の通常モデルで把握できる「Z」の主効果（というか実はX=0のときの条件付き効果なのだがはあきらめるよりない。

結局、求める推定値 $\beta_0 , \gamma_1, \gamma_2$ と切片 $\alpha$ 使って、 $(X, Z)$ の四通りの状況を以下のようにあらわしたいと考えるのである。
$E( Y|X=0, Z =0) = \alpha$
$E( Y|X=1, Z =0) = \alpha + \beta_0$
$E( Y|X=0, Z =1) = \alpha + \gamma_2$
$E( Y|X=1, Z =1) = \alpha + \beta_0+ \gamma_1$

よって、 $(X=0, Z =1)$ のときは1,それ以外は0をとる変数 $U_1$ と、 $(X=1, Z =1)$ のときは1,それ以外の場合は0をとる変数 $U_2$ を作って、Xとともに回帰式に投入すればよいです。

dta1 <- dplyr::mutate( dta1 , gamma1 = if_else( X==1 & Z==1 , 1 ,0 ), gamma2 = if_else( X==0 & Z==1 , 1,0))
lm_res2 <- lm( Y ~ X+ gamma1 + gamma2 , data = dta1)

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  30.0541     0.9335  32.194  < 2e-16 ***
#   X            20.0450     1.3202  15.183  < 2e-16 ***
#   gamma1       14.8049     1.3202  11.214  < 2e-16 ***
#   gamma2       -9.5118     1.3202  -7.205 2.96e-12 ***

$14.8049 -(-9.5118 )= 24.3167$ なので、想定通り $\gamma_1 - \gamma_2 = \gamma$ になっている。

Zが連続変数の場合はどう考えるとよいか

状況をより複雑にしよう。さて、Zが連続変数の場合はどう考えるとよいだろうか。
さっきは $Z$ を朝ごはんをたべていたかどうか(二値変数)、としたが、今回は平日の自宅での平均勉強時間、にしよう（表記は $H$ とする）。

通塾の効果は+20点だが、相乗効果があって勉強時間が1hあがるごとに、効果は+3.5点されると考える。
勉強自体の効果は1hごとに+3.0点としよう。

仮想データをつくり、そして通常の交互作用モデルのOLSを推定する。

cpt <- rnorm( 200, mean=30 , sd = 0.5) #切片（30点とする)
juku <-  sapply( rnorm( 200) , function( x){
  return( if_else( x >=0 , 1 , 0))
}) #塾にいくか否か
std_h <-  rpois( 200 , 3)/2 #勉強時間（Poisson分布に従うと仮定）
errT <- rnorm( 200, sd= 2.5)#誤差項 
#点数
tensu <- cpt + 20 * juku +  3 *std_h +  3.5 * juku * std_h + errT
test_dta <- data.frame( tensu , juku , std_h )

test_res1 <- lm( tensu ~ juku + std_h + std_h:juku, data = test_dta)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  31.0673     0.5222  59.491  < 2e-16 ***
#   juku         19.3359     0.7091  27.270  < 2e-16 ***
#   std_h         2.7675     0.3151   8.784 7.90e-16 ***
#   juku:std_h    3.4089     0.4247   8.026 9.07e-14 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$\gamma = 3.4089$ である
さて、ここで第三の変数が連続変数である本ケースにおいても、「天井」と「床」それぞれの効果を識別したうえで推定値が得たい。

さきほどの離散変数の場合と同じく、
$Y_i = \beta_0 X_i + \gamma_1 X_i H_i + \gamma_2 (1-X_i) H_i$
といった形での定式化をする。

test_dta <- dplyr::mutate( test_dta , U0 = std_h *juku , U1 = std_h*(1- juku))
test_res2 <- lm( tensu ~ juku + U0 + U1, data = test_dta)

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  31.0673     0.5222  59.491  < 2e-16 ***
#   juku         19.3359     0.7091  27.270  < 2e-16 ***
#   U0            6.1765     0.2848  21.685  < 2e-16 ***
#   U1            2.7675     0.3151   8.784  7.9e-16 ***