交互作用を考慮しなくてもよいのか？

　多重線形回帰の項で、交互作用をすべて考慮するのは難しく、しばしば交互作用を考慮しない解析も行われることを紹介した。実際のところ、それはどれくらい妥当な解析になるのかをシミュレーションしてみよう。

交互作用があるときの交互作用を考慮しない推定

　以下のような設定でデータを生成して、交互作用を考慮して推定をしてみる。

------------------------------------------------------

library(plotn)

b1 <- 10

b2 <- -4

b3 <- 7

b4 <- 1

s <- 5

n <- 60

x1 <- runif(n, -20, 20)

x2 <- runif(n, -20, 20)

y <- b4 * x1 * x2 + b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

d <- data.frame(x1 = x1, x2 = x2, y = y)

head(d)

##            x1          x2          y

## 1  12.5349397   0.8149608  -19.50287

## 2  -2.1099159 -16.9715837  -57.39466

## 3   0.5116144   6.6575488   62.99118

## 4  -9.2932711   5.0021768   41.31684

## 5 -18.7211646 -15.9722696  264.85309

## 6   7.5937164 -13.2066298 -211.76806

fit1 <- lm(y ~ x1 * x2, data = d) #線形回帰

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 * x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -13.2058  -2.9061  -0.2171   4.0197   9.6583 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  9.665456   0.674810   14.32   <2e-16 ***

## x1          -4.131569   0.064722  -63.84   <2e-16 ***

## x2           6.909432   0.059591  115.95   <2e-16 ***

## x1:x2        0.991207   0.005845  169.58   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5.171 on 56 degrees of freedom

## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9987 

## F-statistic: 1.528e+04 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

順当に解析ができている。

------------------------------------------------------

coefs <- coef(fit2)[,1]

x1s <- quantile(d$x1, c(0.33/2, 0.5, 0.5+0.66/2))

x2s <- quantile(d$x2, c(0.33/2, 0.5, 0.5+0.66/2))

d <- d[order(d$x1), ]

d$x1_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x1_lv <- factor(d$x1_lv, levels = c("low","medium","high"))

d <- d[order(d$x2), ]

d$x2_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x2_lv <- factor(d$x2_lv, levels = c("low","medium","high"))

plotn(y ~ x1, data = d, group = "x2_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x2") #図1の描画

overdraw(abline(a =  coefs[1] + coefs[3] * x2s[1], b = coefs[2] + coefs[4] * x2s[1], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[3] * x2s[2], b = coefs[2] + coefs[4] * x2s[2]),

         abline(a =  coefs[1] + coefs[3] * x2s[3], b = coefs[2] + coefs[4] * x2s[3], 

                col = "red"))

plotn(y ~ x2, data = d, group = "x1_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x1") #図2の描画

overdraw(abline(a =  coefs[1] + coefs[2] * x1s[1], b = coefs[3] + coefs[4] * x1s[1], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[2] * x1s[2], b = coefs[3] + coefs[4] * x1s[2]),

         abline(a =  coefs[1] + coefs[2] * x1s[3], b = coefs[3] + coefs[4] * x1s[3], 

                col = "red"))

------------------------------------------------------

さて、これを交互作用を考慮せず解析してみよう。

------------------------------------------------------

fit3 <- lm(y ~ x1 + x2, data = d) #交互作用を考慮しない回帰

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1 + x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -266.648  -64.712    2.293   65.634  235.802 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)    4.959     15.159   0.327 0.744755    

## x1            -5.864      1.437  -4.081 0.000141 ***

## x2             5.365      1.324   4.052 0.000156 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 116.3 on 57 degrees of freedom

## Multiple R-squared:  0.372,  Adjusted R-squared:   0.35 

## F-statistic: 16.88 on 2 and 57 DF,  p-value: 1.746e-06

------------------------------------------------------

当然のことながら係数が変わる。これがどのような推定しているかを作図して確認しよう。

------------------------------------------------------

coefs <- coef(fit4)[,1]

x1s <- quantile(d$x1, c(0.33/2, 0.5, 0.5+0.66/2))

x2s <- quantile(d$x2, c(0.33/2, 0.5, 0.5+0.66/2))

d <- d[order(d$x1), ]

d$x1_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x1_lv <- factor(d$x1_lv, levels = c("low","medium","high"))

d <- d[order(d$x2), ]

d$x2_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x2_lv <- factor(d$x2_lv, levels = c("low","medium","high"))

plotn(y ~ x1, data = d, group = "x2_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x2") #図3の描画

overdraw(abline(a =  coefs[1] + coefs[3] * x2s[1], b = coefs[2], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[3] * x2s[2], b = coefs[2]),

         abline(a =  coefs[1] + coefs[3] * x2s[3], b = coefs[2], 

                col = "red"))

plotn(y ~ x2, data = d, group = "x1_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x1") #図4の描画

overdraw(abline(a =  coefs[1] + coefs[2] * x1s[1], b = coefs[3], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[2] * x1s[2], b = coefs[3]),

         abline(a =  coefs[1] + coefs[2] * x1s[3], b = coefs[3], 

                col = "red"))

------------------------------------------------------

当然、真のモデルは異なっているので、推定としては正しくないのは当然なのだが、傾向として読み取れるのは、y ~ x1の回帰直線はx2の中央値周りで精度が良く、中央値から離れるほど精度が悪くなるということ、y ~ x2ではx1の中央値周りで精度が良くなる、ということである。言い換えれば、重回帰で求まったx1の係数はx2の中央値周り、x2の係数のx1の中央値周りの推定値となっているということである。交互作用を考慮すればそれはデータ生成過程と一致するので、説明変数外の予測も可能だが、考慮しない場合は得られている説明変数周りでの近似として推定されている、ともいえるだろう。

　また、図から読み取れることとして、変数の定義域内に回帰直線の交点があるとき、あるいはデータが蝶ネクタイ型に分布しているとき、解析は容易に失敗するだろうことが予測できる。だから、設定するパラメータは同じでも、x1やx2の定義域を交点の外側に設定すると、推定はおおむねうまくいく。まず、交互作用を考慮した場合、

------------------------------------------------------

b1 <- 10

b2 <- -4

b3 <- 7

b4 <- 1

s <- 5

n <- 60

x1 <- runif(n, 10, 50)

x2 <- runif(n, 10, 50)

y <- b4 * x1 * x2 + b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

d <- data.frame(x1 = x1, x2 = x2, y = y)

head(d)

##         x1       x2         y

## 1 19.24801 17.15594  382.9880

## 2 36.32044 33.37036 1319.7946

## 3 31.59725 32.55700 1127.8097

## 4 37.51360 19.18257  710.6238

## 5 18.75059 41.10368  993.1381

## 6 17.20444 44.67140 1028.1602

fit1 <- lm(y ~ x1 * x2, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 * x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -13.3935  -3.0150  -0.4684   2.2259  15.2094 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept) 19.218942   6.779806   2.835  0.00637 ** 

## x1          -4.298360   0.201689 -21.312  < 2e-16 ***

## x2           6.846338   0.192156  35.629  < 2e-16 ***

## x1:x2        1.006020   0.005665 177.590  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 4.89 on 56 degrees of freedom

## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 

## F-statistic: 2.523e+05 on 3 and 56 DF,  p-value: < 2.2e-16

coefs <- coef(fit2)[,1]

x1s <- quantile(d$x1, c(0.33/2, 0.5, 0.5+0.66/2))

x2s <- quantile(d$x2, c(0.33/2, 0.5, 0.5+0.66/2))

d <- d[order(d$x1), ]

d$x1_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x1_lv <- factor(d$x1_lv, levels = c("low","medium","high"))

d <- d[order(d$x2), ]

d$x2_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x2_lv <- factor(d$x2_lv, levels = c("low","medium","high"))

plotn(y ~ x1, data = d, group = "x2_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x2") #図5の作図

overdraw(abline(a =  coefs[1] + coefs[3] * x2s[1], b = coefs[2] + coefs[4] * x2s[1], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[3] * x2s[2], b = coefs[2] + coefs[4] * x2s[2]),

         abline(a =  coefs[1] + coefs[3] * x2s[3], b = coefs[2] + coefs[4] * x2s[3], 

                col = "red"))

plotn(y ~ x2, data = d, group = "x1_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x1") #図6の作図

overdraw(abline(a =  coefs[1] + coefs[2] * x1s[1], b = coefs[3] + coefs[4] * x1s[1], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[2] * x1s[2], b = coefs[3] + coefs[4] * x1s[2]),

         abline(a =  coefs[1] + coefs[2] * x1s[3], b = coefs[3] + coefs[4] * x1s[3], 

                col = "red"))

------------------------------------------------------

交互作用を考慮しない場合、

------------------------------------------------------

fit3 <- lm(y ~ x1 + x2, data = d) #線形回帰

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1 + x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -215.635  -78.849   -7.802   61.545  302.743 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -1082.944     64.254  -16.85   <2e-16 ***

## x1             29.847      1.434   20.82   <2e-16 ***

## x2             39.395      1.359   28.99   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 115.1 on 57 degrees of freedom

## Multiple R-squared:  0.9583, Adjusted R-squared:  0.9568 

## F-statistic: 654.5 on 2 and 57 DF,  p-value: < 2.2e-16

coefs <- coef(fit4)[,1]

x1s <- quantile(d$x1, c(0.33/2, 0.5, 0.5+0.66/2))

x2s <- quantile(d$x2, c(0.33/2, 0.5, 0.5+0.66/2))

d <- d[order(d$x1), ]

d$x1_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x1_lv <- factor(d$x1_lv, levels = c("low","medium","high"))

d <- d[order(d$x2), ]

d$x2_lv <- c(rep("low",20),rep("medium",20),rep("high",20))

d$x2_lv <- factor(d$x2_lv, levels = c("low","medium","high"))

plotn(y ~ x1, data = d, group = "x2_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x2") #図7の作図

overdraw(abline(a =  coefs[1] + coefs[3] * x2s[1], b = coefs[2], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[3] * x2s[2], b = coefs[2]),

         abline(a =  coefs[1] + coefs[3] * x2s[3], b = coefs[2], 

                col = "red"))

plotn(y ~ x2, data = d, group = "x1_lv", mode = "m", pch = 15:17,

      col.dot = c("blue","black","red"), legend = T, lty.leg = 1, leg.sp = 5,

      leg.title = "x1") #図8の作図

overdraw(abline(a =  coefs[1] + coefs[2] * x1s[1], b = coefs[3], 

                col = "blue"),

         abline(a =  coefs[1] + coefs[2] * x1s[2], b = coefs[3]),

         abline(a =  coefs[1] + coefs[2] * x1s[3], b = coefs[3], 

                col = "red"))

------------------------------------------------------

このようにずれはあるが、図3、4ほどひどい状況ではない。ただし、データ生成過程の時のパラメータが推定されるわけではないことはよく覚えておくべきだ。例えば、x1の係数の設定は-4だったのに、交互作用を考慮しない場合は29.8となっている。これも、得られている説明変数周りでの推定だから、ずれていて当然である。感覚的な話であるが、データの分布を見て、明らかな交互作用、蝶ネクタイ型のデータ分布が見られれば、それらの変数間の交互作用を警戒して、そうでなければ交互作用を考慮しない解析を選択するだろうか。後述するように、データの分布をみることはもっと多数の変数を扱うときに重要になる。

交互作用は交互作用のない変数の推定にも影響する

　次は交互作用のある変数だけでなく、交互作用のない変数も混ぜてみよう。x3とx4の間だけ交互作用を仮定し、残りには仮定しない。データを生成したのち、総当たりプロットを描いてみる。

------------------------------------------------------

b1 <- 10

b2 <- -4

b3 <- 2

b4 <- 1

b5 <- -1

b6 <- 0

b3_4 <- 0.1

s <- 5

n <- 60

x1 <- runif(n, -20, 20)

x2 <- runif(n, -20, 20)

x3 <- runif(n, -20, 20)

x4 <- runif(n, -20, 20)

x5 <- runif(n, -20, 20)

y <- b3_4 * x3 * x4 + b6 * x5 +  b5 * x4 + b4 * x3 +  b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

d <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, x5 = x5, y = y)

head(d)

##           x1         x2         x3         x4          x5          y

## 1  18.491446  0.5388135 -16.572097  10.703666  10.8071188 -108.78102

## 2  -5.449471  6.5626611  11.672179   5.756849 -10.6390651   55.61920

## 3  -1.298966 -4.2627077   9.495639  17.607186   6.9212634   17.48456

## 4 -11.850411 18.6771198  14.676023  -1.121954   0.9132696  113.54449

## 5  -1.808212 15.8076049 -15.083052 -14.086600  -2.5212866   70.08076

## 6 -16.311126 -8.8528322   9.922507 -18.340546  -1.9552636   65.80573

fit1 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -26.117  -8.439   1.398   8.316  28.973 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  9.81174    1.87603   5.230 2.82e-06 ***

## x1          -4.15840    0.15473 -26.875  < 2e-16 ***

## x2           1.87870    0.17811  10.548 1.00e-14 ***

## x3           0.98611    0.15540   6.346 4.73e-08 ***

## x4          -0.48685    0.15419  -3.158   0.0026 ** 

## x5           0.01355    0.16130   0.084   0.9334    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 13.16 on 54 degrees of freedom

## Multiple R-squared:  0.9486, Adjusted R-squared:  0.9439 

## F-statistic: 199.5 on 5 and 54 DF,  p-value: < 2.2e-16

fit3 <- lm(y ~ x1 + x2 + x3*x4 + x5, data = d)

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1 + x2 + x3 * x4 + x5, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -14.530  -2.830  -0.382   3.062  12.416 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  9.450343   0.703691  13.430   <2e-16 ***

## x1          -4.036521   0.058401 -69.118   <2e-16 ***

## x2           1.899593   0.066791  28.441   <2e-16 ***

## x3           0.964033   0.058278  16.542   <2e-16 ***

## x4          -0.943582   0.063025 -14.971   <2e-16 ***

## x5          -0.051986   0.060586  -0.858    0.395    

## x3:x4        0.092135   0.005063  18.196   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 4.936 on 53 degrees of freedom

## Multiple R-squared:  0.9929, Adjusted R-squared:  0.9921 

## F-statistic:  1238 on 6 and 53 DF,  p-value: < 2.2e-16

plotn(d) #図9の作図

------------------------------------------------------

総当たりプロットのうち、最下段あるいは最右列に各説明変数と被説明変数の関係が描かれている。相対的に交互作用が弱いため、x3-yのプロットやx4-yのプロットを確認しても、蝶ネクタイ型の分布は見られない。こんな時、交互作用を含まないモデルと含むモデルで検出力がどれくらい変わるか確認してみよう。

------------------------------------------------------

#simulate

p1_1 <- NULL

p2_1 <- NULL

p3_1 <- NULL

p4_1 <- NULL

p5_1 <- NULL

p6_1 <- NULL

b1_1 <- NULL

b2_1 <- NULL

b3_1 <- NULL

b4_1 <- NULL

b5_1 <- NULL

b6_1 <- NULL

p1_2 <- NULL

p2_2 <- NULL

p3_2 <- NULL

p4_2 <- NULL

p5_2 <- NULL

p6_2 <- NULL

b1_2 <- NULL

b2_2 <- NULL

b3_2 <- NULL

b4_2 <- NULL

b5_2 <- NULL

b6_2 <- NULL

for(i in 1:10000){

  x1 <- runif(n, -20, 20)

  x2 <- runif(n, -20, 20)

  x3 <- runif(n, -20, 20)

  x4 <- runif(n, -20, 20)

  x5 <- runif(n, -20, 20)

  y <- b3_4 * x3 * x4 + b6 * x5 +  b5 * x4 + b4 * x3 +  b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

  d <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, x5 = x5, y = y)

  fit1 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = d) #交互作用を含まない

  fit2 <- summary(fit1)

  fit3 <- lm(y ~ x1 + x2 + x3*x4 + x5, data = d) #交互作用を含む

  fit4 <- summary(fit3)

  p1_1 <- c(p1_1, coef(fit2)[1,4]) #交互作用を含まないb1のP値

  p2_1 <- c(p2_1, coef(fit2)[2,4]) #交互作用を含まないb2のP値

  p3_1 <- c(p3_1, coef(fit2)[3,4]) #交互作用を含まないb3のP値

  p4_1 <- c(p4_1, coef(fit2)[4,4]) #交互作用を含まないb4のP値

  p5_1 <- c(p5_1, coef(fit2)[5,4]) #交互作用を含まないb5のP値

  p6_1 <- c(p6_1, coef(fit2)[6,4]) #交互作用を含まないb6のP値

  p1_2 <- c(p1_2, coef(fit4)[1,4]) #交互作用を含むb1のP値

  p2_2 <- c(p2_2, coef(fit4)[2,4]) #交互作用を含むb2のP値

  p3_2 <- c(p3_2, coef(fit4)[3,4]) #交互作用を含むb3のP値

  p4_2 <- c(p4_2, coef(fit4)[4,4]) #交互作用を含むb4のP値

  p5_2 <- c(p5_2, coef(fit4)[5,4]) #交互作用を含むb5のP値

  p6_2 <- c(p6_2, coef(fit4)[6,4]) #交互作用を含むb6のP値

  b1_1 <- c(b1_1, coef(fit2)[1,1]) #交互作用を含まないb1値

  b2_1 <- c(b2_1, coef(fit2)[2,1]) #交互作用を含まないb2値

  b3_1 <- c(b3_1, coef(fit2)[3,1]) #交互作用を含まないb3値

  b4_1 <- c(b4_1, coef(fit2)[4,1]) #交互作用を含まないb4値

  b5_1 <- c(b5_1, coef(fit2)[5,1]) #交互作用を含まないb5値

  b6_1 <- c(b6_1, coef(fit2)[6,1]) #交互作用を含まないb6値

  b1_2 <- c(b1_2, coef(fit4)[1,1]) #交互作用を含むb1値

  b2_2 <- c(b2_2, coef(fit4)[2,1]) #交互作用を含むb2値

  b3_2 <- c(b3_2, coef(fit4)[3,1]) #交互作用を含むb3値

  b4_2 <- c(b4_2, coef(fit4)[4,1]) #交互作用を含むb4値

  b5_2 <- c(b5_2, coef(fit4)[5,1]) #交互作用を含むb5値

  b6_2 <- c(b6_2, coef(fit4)[6,1]) #交互作用を含むb6値

}

histn(p1_1, xlab = "P value", ylab = "頻度") #図10の作図

histn(p2_1, xlab = "P value", ylab = "頻度") #図11の作図

histn(p3_1, xlab = "P value", ylab = "頻度") #図12の作図

histn(p4_1, xlab = "P value", ylab = "頻度") #図13の作図

histn(p5_1, xlab = "P value", ylab = "頻度") #図14の作図

histn(p6_1, xlab = "P value", ylab = "頻度") #図15の作図

------------------------------------------------------

交互作用を考慮した時のP値のヒストグラムは、

------------------------------------------------------

histn(p1_2, xlab = "P value", ylab = "頻度") #図16の作図

histn(p2_2, xlab = "P value", ylab = "頻度") #図17の作図

histn(p3_2, xlab = "P value", ylab = "頻度") #図18の作図

histn(p4_2, xlab = "P value", ylab = "頻度") #図19の作図

histn(p5_2, xlab = "P value", ylab = "頻度") #図20の作図

histn(p6_2, xlab = "P value", ylab = "頻度") #図21の作図

------------------------------------------------------

また、各P値を比較するが、併せてデータ生成過程のパラメータの正負と一致しているかも、同時に考慮しよう。数値そのものはデータ生成過程の時から大きくずれていても仕方ないが、せめて、その要因が正の効果を持つか負の効果を持つかくらいは一致していてほしいだろう。

------------------------------------------------------

#交互作用を考慮しなかったとき

sum((p1_1 < 0.05) * (b1_1 > 0))

## [1] 9995

sum((p2_1 < 0.05) * (b2_1 < 0))

## [1] 10000

sum((p3_1 < 0.05) * (b3_1 > 0))

## [1] 10000

sum((p4_1 < 0.05) * (b4_1 > 0))

## [1] 9988

sum((p5_1 < 0.05) * (b5_1 < 0))

## [1] 9995

sum(p6_1 < 0.05)

## [1] 508

#交互作用を考慮したとき

sum((p1_2 < 0.05) * (b1_2 > 0))

## [1] 10000

sum((p2_2 < 0.05) * (b2_2 < 0))

## [1] 10000

sum((p3_2 < 0.05) * (b3_2 > 0))

## [1] 10000

sum((p4_2 < 0.05) * (b4_2 > 0))

## [1] 10000

sum((p5_2 < 0.05) * (b5_2 < 0))

## [1] 10000

sum(p6_2 < 0.05)

## [1] 494

------------------------------------------------------

確認をすると、切片(b1)、x1（b2）およびx2（b3）に関して、検出力は100%近いが、交互作用を考慮しないとわずかに切片での検出力が落ちていることがわかる。この後も紹介するが、交互作用がさらに強くなると、交互作用を考慮しない場合、検出力がさらに落ちる。x5はb6 = 0として設定しているので、有意差が出ないことが望ましいのだが、危険率は5%で変わらないことがわかる。交互作用しているx3（b4）とx4（b5）はそもそもモデルが間違っているし、交互作用がある場合にはもとのパラメータ設定に関わらず、交互作用するもう一方の変数の大きさや定義域によって効果は正にも負にもなるので、あまり考えても仕方ない面がある。参考程度に、というところだろうか。

　では、下記のように設定して交互作用を強くしよう。

------------------------------------------------------

b1 <- 10

b2 <- -4

b3 <- 2

b4 <- 1

b5 <- -1

b6 <- 0

b3_4 <- 1

s <- 5

n <- 60

x1 <- runif(n, -20, 20)

x2 <- runif(n, -20, 20)

x3 <- runif(n, -20, 20)

x4 <- runif(n, -20, 20)

x5 <- runif(n, -20, 20)

y <- b3_4 * x3 * x4 + b6 * x5 +  b5 * x4 + b4 * x3 +  b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

d <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, x5 = x5, y = y)

head(d)

##           x1         x2        x3         x4         x5          y

## 1 -15.738834 -10.822455 -7.174249 -12.249903 -10.470701 138.892863

## 2  14.198674  -8.706631 10.394956  14.274739  16.780552  80.566917

## 3  -3.030901  -7.655370  2.375297 -18.788921   2.958674  -5.773333

## 4 -12.967925 -12.288806  5.099475   9.964350  -5.028253  79.127102

## 5 -14.905394 -18.194765 10.846586  -7.822847  -9.839541 -30.861562

## 6 -18.187971  -8.370323 -3.569154  12.473597  19.545878   9.924296

fit1 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -211.178  -47.550    0.128   61.603  295.207 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)  

## (Intercept)  35.7308    14.4386   2.475   0.0165 *

## x1           -2.6964     1.3249  -2.035   0.0467 *

## x2            1.1774     1.1672   1.009   0.3176  

## x3           -1.0911     1.4505  -0.752   0.4552  

## x4            2.5462     1.2321   2.067   0.0436 *

## x5            0.8177     1.1447   0.714   0.4781  

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 107.1 on 54 degrees of freedom

## Multiple R-squared:  0.1706, Adjusted R-squared:  0.09379 

## F-statistic: 2.221 on 5 and 54 DF,  p-value: 0.06532

fit3 <- lm(y ~ x1 + x2 + x3*x4 + x5, data = d)

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1 + x2 + x3 * x4 + x5, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -12.0322  -3.0128   0.2809   3.1530  12.8173 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept) 10.473483   0.678071  15.446   <2e-16 ***

## x1          -4.012373   0.061076 -65.694   <2e-16 ***

## x2           2.024699   0.053584  37.786   <2e-16 ***

## x3           0.917416   0.067432  13.605   <2e-16 ***

## x4          -0.863810   0.060152 -14.360   <2e-16 ***

## x5           0.084040   0.052492   1.601    0.115    

## x3:x4        0.994487   0.006189 160.692   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 4.895 on 53 degrees of freedom

## Multiple R-squared:  0.9983, Adjusted R-squared:  0.9981 

## F-statistic:  5191 on 6 and 53 DF,  p-value: < 2.2e-16

plotn(d) #図22の作図

------------------------------------------------------

相対的に交互作用が強く、x3-yのプロットやx4-yに蝶ネクタイ型の分布が見られる。こんな時の検出力は以下の通り。

------------------------------------------------------

#simulate

p1_1 <- NULL

p2_1 <- NULL

p3_1 <- NULL

p4_1 <- NULL

p5_1 <- NULL

p6_1 <- NULL

b1_1 <- NULL

b2_1 <- NULL

b3_1 <- NULL

b4_1 <- NULL

b5_1 <- NULL

b6_1 <- NULL

p1_2 <- NULL

p2_2 <- NULL

p3_2 <- NULL

p4_2 <- NULL

p5_2 <- NULL

p6_2 <- NULL

b1_2 <- NULL

b2_2 <- NULL

b3_2 <- NULL

b4_2 <- NULL

b5_2 <- NULL

b6_2 <- NULL

for(i in 1:10000){

  x1 <- runif(n, -20, 20)

  x2 <- runif(n, -20, 20)

  x3 <- runif(n, -20, 20)

  x4 <- runif(n, -20, 20)

  x5 <- runif(n, -20, 20)

  y <- b3_4 * x3 * x4 + b6 * x5 +  b5 * x4 + b4 * x3 +  b3 * x2 + b2 * x1 + rnorm(n = n, mean = b1, sd = s)

  d <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, x5 = x5, y = y)

  fit1 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = d) #交互作用を含まない

  fit2 <- summary(fit1)

  fit3 <- lm(y ~ x1 + x2 + x3*x4 + x5, data = d) #交互作用を含む

  fit4 <- summary(fit3)

  p1_1 <- c(p1_1, coef(fit2)[1,4]) #交互作用を含まないb1のP値

  p2_1 <- c(p2_1, coef(fit2)[2,4]) #交互作用を含まないb2のP値

  p3_1 <- c(p3_1, coef(fit2)[3,4]) #交互作用を含まないb3のP値

  p4_1 <- c(p4_1, coef(fit2)[4,4]) #交互作用を含まないb4のP値

  p5_1 <- c(p5_1, coef(fit2)[5,4]) #交互作用を含まないb5のP値

  p6_1 <- c(p6_1, coef(fit2)[6,4]) #交互作用を含まないb6のP値

  p1_2 <- c(p1_2, coef(fit4)[1,4]) #交互作用を含むb1のP値

  p2_2 <- c(p2_2, coef(fit4)[2,4]) #交互作用を含むb2のP値

  p3_2 <- c(p3_2, coef(fit4)[3,4]) #交互作用を含むb3のP値

  p4_2 <- c(p4_2, coef(fit4)[4,4]) #交互作用を含むb4のP値

  p5_2 <- c(p5_2, coef(fit4)[5,4]) #交互作用を含むb5のP値

  p6_2 <- c(p6_2, coef(fit4)[6,4]) #交互作用を含むb6のP値

  b1_1 <- c(b1_1, coef(fit2)[1,1]) #交互作用を含まないb1値

  b2_1 <- c(b2_1, coef(fit2)[2,1]) #交互作用を含まないb2値

  b3_1 <- c(b3_1, coef(fit2)[3,1]) #交互作用を含まないb3値

  b4_1 <- c(b4_1, coef(fit2)[4,1]) #交互作用を含まないb4値

  b5_1 <- c(b5_1, coef(fit2)[5,1]) #交互作用を含まないb5値

  b6_1 <- c(b6_1, coef(fit2)[6,1]) #交互作用を含まないb6値

  b1_2 <- c(b1_2, coef(fit4)[1,1]) #交互作用を含むb1値

  b2_2 <- c(b2_2, coef(fit4)[2,1]) #交互作用を含むb2値

  b3_2 <- c(b3_2, coef(fit4)[3,1]) #交互作用を含むb3値

  b4_2 <- c(b4_2, coef(fit4)[4,1]) #交互作用を含むb4値

  b5_2 <- c(b5_2, coef(fit4)[5,1]) #交互作用を含むb5値

  b6_2 <- c(b6_2, coef(fit4)[6,1]) #交互作用を含むb6値

}

histn(p1_1, xlab = "P value", ylab = "頻度") #図23の作図

histn(p2_1, xlab = "P value", ylab = "頻度") #図24の作図

histn(p3_1, xlab = "P value", ylab = "頻度") #図25の作図

histn(p4_1, xlab = "P value", ylab = "頻度") #図26の作図

histn(p5_1, xlab = "P value", ylab = "頻度") #図27の作図

histn(p6_1, xlab = "P value", ylab = "頻度") #図28の作図

------------------------------------------------------

交互作用を考慮した時のP値のヒストグラムは、

------------------------------------------------------

histn(p1_2, xlab = "P value", ylab = "頻度") #図29の作図

histn(p2_2, xlab = "P value", ylab = "頻度") #図30の作図

histn(p3_2, xlab = "P value", ylab = "頻度") #図31の作図

histn(p4_2, xlab = "P value", ylab = "頻度") #図32の作図

histn(p5_2, xlab = "P value", ylab = "頻度") #図33の作図

histn(p6_2, xlab = "P value", ylab = "頻度") #図34の作図

------------------------------------------------------

また、各P値を比較するが、併せてデータ生成過程のパラメータの正負と一致しているかも、同時に考慮しよう。数値そのものはデータ生成過程の時から大きくずれていても仕方ないが、せめて、その要因が正の効果を持つか負の効果を持つかくらいは一致していてほしいだろう。

------------------------------------------------------

#交互作用を考慮しなかったとき

sum((p1_1 < 0.05) * (b1_1 > 0))

## [1] 908

sum((p2_1 < 0.05) * (b2_1 < 0))

## [1] 7177

sum((p3_1 < 0.05) * (b3_1 > 0))

## [1] 2524

sum((p4_1 < 0.05) * (b4_1 > 0))

## [1] 1593

sum((p5_1 < 0.05) * (b5_1 < 0))

## [1] 1672

sum(p6_1 < 0.05)

## [1] 518

#交互作用を考慮したとき

sum((p1_2 < 0.05) * (b1_2 > 0))

## [1] 10000

sum((p2_2 < 0.05) * (b2_2 < 0))

## [1] 10000

sum((p3_2 < 0.05) * (b3_2 > 0))

## [1] 10000

sum((p4_2 < 0.05) * (b4_2 > 0))

## [1] 10000

sum((p5_2 < 0.05) * (b5_2 < 0))

## [1] 10000

sum(p6_2 < 0.05)

## [1] 501

------------------------------------------------------

確認をすると、切片(b1)、x1（b2）およびx2（b3）に関して交互作用を考慮しないことで、極端に検出力が落ちることがわかる。これくらいの交互作用の強さになると警戒をするべき、ということだろう。x5の危険率はあまり影響を受けない。

　ただ、いろいろシミュレートしてみるとわかるが、交互作用項の動態は極めて複雑だ。交互作用項の絶対値の大きさではなく、交互作用項と単純主効果の相対的な大きさによって、単純主効果の検出力は大きく変わる。解析前に総当たりプロットを確認して、変数間の様子を確認することは必須だろう。これは多重共線性の確認にも役立つ。多数の変数を組み込んで解析をするとき、その解析は検出力が落ちていることを念頭において解釈をするべきであろう。

　久保先生の著書『データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC』(pp.127-130)によれば、交互作用項をむやみに入れることは避けることが推奨されている。AICの基準などでは交互作用項を含むモデルが選択される傾向があるが、それは真の交互作用を示しているか不明である点、場所差などを考慮した一般化線形混合モデルを活用することでその交互作用の効果が消える点などを理由として挙げている。