交互作用を

考慮しない重回帰？

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

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

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

・交差が定義域内で起きるときに交互作用を考慮

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

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

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"))

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

図1 x1-y平面上からみたデータの様子と特定のx2に対する交互作用を考慮した回帰直線

図2 x2-y平面上からみたデータの様子と特定のx1に対する交互作用を考慮した回帰直線

・交差が定義域内で起きるときに交互作用を考慮しない

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

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

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"))

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

図3 x1-y平面上からみたデータの様子と特定のx2に対する交互作用を考慮しない回帰直線

図4 x2-y平面上からみたデータの様子と特定のx1に対する交互作用を考慮しない回帰直線

当然、真のモデルは異なっているので、推定としては正しくないのは当然なのだが、傾向として読み取れるのは、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"))

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

図5 x1-y平面上からみたデータの様子と特定のx2に対する交互作用を考慮した回帰直線

図6 x2-y平面上からみたデータの様子と特定のx1に対する交互作用を考慮した回帰直線

・交差が定義域外で起きるときに交互作用を考慮しない

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

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

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"))

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

図7 x1-y平面上からみたデータの様子と特定のx2に対する交互作用を考慮しない回帰直線

図8 x2-y平面上からみたデータの様子と特定のx1に対する交互作用を考慮しない回帰直線

このようにずれはあるが、図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の作図

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

図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の作図

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

図10 交互作用を考慮しないときのb1のP値

図11 交互作用を考慮しないときのb2のP値

図12 交互作用を考慮しないときのb3のP値

図13 交互作用を考慮しないときのb4のP値

図14 交互作用を考慮しないときのb5のP値

図15 交互作用を考慮しないときのb6のP値

交互作用を考慮した時の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の作図

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

図16 交互作用を考慮したときのb1のP値

図17 交互作用を考慮したときのb2のP値

図18 交互作用を考慮したときのb3のP値

図19 交互作用を考慮したときのb4のP値

図20 交互作用を考慮したときのb5のP値

図21 交互作用を考慮したときのb6のP値

また、各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の作図

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

図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の作図

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

図23 交互作用を考慮しないときのb1のP値

図24 交互作用を考慮しないときのb2のP値

図25 交互作用を考慮しないときのb3のP値

図26 交互作用を考慮しないときのb4のP値

図27 交互作用を考慮しないときのb5のP値

図28 交互作用を考慮しないときのb6のP値

交互作用を考慮した時の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の作図

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

図29 交互作用を考慮したときのb1のP値

図30 交互作用を考慮したときのb2のP値

図31 交互作用を考慮したときのb3のP値

図32 交互作用を考慮したときのb4のP値

図33 交互作用を考慮したときのb5のP値

図34 交互作用を考慮したときのb6のP値

また、各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の基準などでは交互作用項を含むモデルが選択される傾向があるが、それは真の交互作用を示しているか不明である点、場所差などを考慮した一般化線形混合モデルを活用することでその交互作用の効果が消える点などを理由として挙げている。

交互作用を

考慮しない重回帰？

目次