分散分析と線形回帰の組み合わせ技

　これまでカテゴリカル変数の解析は分散分析、連続変数の解析は線形回帰、と区別して紹介してきた。実はこれらの変数は、一緒に説明変数に組み込むことが可能である。つまり、説明変数にカテゴリカル変数と連続変数が混在していてもよいし、連続変数が複数含まれていてもよい。このような、連続変数×カテゴリカル変数、連続変数×連続変数のような解析を多重線形回帰multiple linear regression、あるいは単に重回帰とよぶ。本項では連続変数とカテゴリカル変数が説明変数の重回帰を紹介し、連続変数が複数含まれる解析は次回に紹介する。

　では、連続変数xと2水準（A、B）を持つカテゴリカル変数gのモデルはどうなるだろうか。線形モデルの表記にならって下記のようになる。カテゴリカル変数は2水準だけなので下記の一本の式で十分だ（Aのときはg = 0、Ｂのときはg = 1）。

y ~ b1 + b2  x + b3  g + b4 (x × g)

この式はどんなデータの状況を想定していることになるだろうか。以下のように変形するとわかりやすくなる。

y ~ (b1 + b3 g) + (b2 + b4 g) x

カテゴリカル変数がAのとき、g = 0を代入すれば、

y ~ b1 + b2  x

カテゴリカル変数がBのとき、g = 1を代入すれば、

y ~ (b1 + b3) + (b2 + b4) x

このように、カテゴリカル変数によって傾きと切片が異なる直線を推定してることになる。もし、交互作用を考慮しない（b4 = 0と仮定）ならば、カテゴリカル変数がAのとき、g = 0を代入すれば、

y ~ b1 + b2  x

カテゴリカル変数がBのとき、g = 1を代入すれば、

y ~ (b1 + b3) + b2 x

となり、2本の直線は切片が異なり、傾きは等しい、平行な直線を推定することになり、以前紹介したような交互作用がなければ平行、あれば平行ではない、という関係が強調されていることになる。このような、交互作用を考慮しない（ないと仮定できる）、連続変数×カテゴリカル変数の多重線形回帰を共分散分析（Analysis of covariance、ANCOVA）と呼ぶこともある。帰無仮説は今までと同様に、残差が平均0の正規分布に従うと仮定し、「b1 = 0、b2 = 0、b3 = 0、b4 = 0」である。

係数の推定

　係数の推定は線形回帰で紹介したように最小二乗法を用いる。線形回帰の時に明らかになったことを使いまわすことで比較的簡単に示すことができる。

これの意味するところは、まず、g = 0のときその水準に含まれる標本を単独で線形回帰することで、b1とb2の推定値を得ることができる。次に、g = 1のとき、その水準に含まれる標本だけで線形回帰をするとb1 + b3を切片、b2 + b4を傾きとする推定値が得られる。ここから、g = 0のときの推定値b1、b2を差し引けばb3、b4が求まるという寸法だ。それほど難しいことではないだろう。

多重線形回帰の検出力

　さて、今回もシミュレートしてみよう。次のようにb1 ≠ 0、b2、b3、b4 = 0のデータを生成して、線形回帰を行う。

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

library(plotn)

cols <- col_genelator(number = 2, alpha = 1)

fils <- col_genelator(number = 2, alpha = 0.5)

b1 <- 10

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 30

x1_1 <- runif(n, 0, 30)

x1_2 <- runif(n, 0, 30)

y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)#水準Aのデータ生成

y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)#水準Bのデータ生成

d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

head(d)

##           x          y g

## 1 14.626422  8.1310194 A

## 2  1.703273 12.1936144 A

## 3 15.573617  8.6202482 A

## 4  9.704097 18.4047000 A

## 5 15.455977  0.7768798 A

## 6 21.048418 17.2785586 A

fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

#連続変数xとカテゴリカル変数gが説明変数になっていることに注目

fit2 <- summary(fit1)

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図1の描画

b1e <- coef(fit2)[1,1] #b1の推定値

b2e <- coef(fit2)[2,1] #b2の推定値

b3e <- coef(fit2)[3,1] #b3の推定値

b4e <- coef(fit2)[4,1] #b4の推定値

overdraw(abline(a = b1e, b = b2e, col = cols[1]),  #水準Aの直線

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2)) #水準Bの直線

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

さて、b1だけが0ではなく、b2 = b3 = b4 = 0を仮定しているので、2水準の直線は切片が同じで、傾きはともに0であることを想定している。実際、図1を見るとそのようなデータ生成になっている。ではこの時の検出力を確かめよう。

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

b1 <- 10

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #b1のP値

  p2 <- c(p2, coef(fit2)[2,4]) #b2のP値

  p3 <- c(p3, coef(fit2)[3,4]) #b3のP値

  p4 <- c(p4, coef(fit2)[4,4]) #b4のP値

}

histn(p1, xlab = "P value", ylab = "頻度") #図2の描画

histn(p2, xlab = "P value", ylab = "頻度") #図3の描画

histn(p3, xlab = "P value", ylab = "頻度") #図4の描画

histn(p4, xlab = "P value", ylab = "頻度") #図5の描画

sum(p1 < 0.05)

## [1] 9985

sum(p2 < 0.05)

## [1] 536

sum(p3 < 0.05)

## [1] 506

sum(p4 < 0.05)

## [1] 553

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

見ると仮定通り、b1の検出力は100%近いのに対して、残りの係数の危険率は5%に抑えられている。次は傾きも足してみよう。ただし、水準間の違いはない（b3 = b4 = 0）。

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

b1 <- 10

b2 <- 2

b3 <- 0

b4 <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #b1のP値

  p2 <- c(p2, coef(fit2)[2,4]) #b2のP値

  p3 <- c(p3, coef(fit2)[3,4]) #b3のP値

  p4 <- c(p4, coef(fit2)[4,4]) #b4のP値

}

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図6の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

histn(p1, xlab = "P value", ylab = "頻度") #図7の描画

histn(p2, xlab = "P value", ylab = "頻度") #図8の描画

histn(p3, xlab = "P value", ylab = "頻度") #図9の描画

histn(p4, xlab = "P value", ylab = "頻度") #図10の描画

sum(p1 < 0.05)

## [1] 9984

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 520

sum(p4 < 0.05)

## [1] 487

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

b1およびb2の検出力は100%近くなった。一方で、b3、b4は危険率5%程度に抑えられている。では、今度は水準間で切片も変えてみよう。ただし、傾きは一緒だ。

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

b1 <- 10

b2 <- 2

b3 <- 5

b4 <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #b1のP値

  p2 <- c(p2, coef(fit2)[2,4]) #b2のP値

  p3 <- c(p3, coef(fit2)[3,4]) #b3のP値

  p4 <- c(p4, coef(fit2)[4,4]) #b4のP値

}

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図11の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

histn(p1, xlab = "P value", ylab = "頻度") #図12の描画

histn(p2, xlab = "P value", ylab = "頻度") #図13の描画

histn(p3, xlab = "P value", ylab = "頻度") #図14の描画

histn(p4, xlab = "P value", ylab = "頻度") #図15の描画

sum(p1 < 0.05)

## [1] 9973

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 4694

sum(p4 < 0.05)

## [1] 514

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

b1およびb2の検出力は100%近くなった。そして、この条件のもとではb3の検出力は47%程度である。そして、b4の危険率は5%程度である。では、今度は水準間で傾きだけを変えてみよう。

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

b1 <- 10

b2 <- 2

b3 <- 0

b4 <- -3

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4]) #b1のP値

  p2 <- c(p2, coef(fit2)[2,4]) #b2のP値

  p3 <- c(p3, coef(fit2)[3,4]) #b3のP値

  p4 <- c(p4, coef(fit2)[4,4]) #b4のP値

}

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図16の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

histn(p1, xlab = "P value", ylab = "頻度") #図17の描画

histn(p2, xlab = "P value", ylab = "頻度") #図18の描画

histn(p3, xlab = "P value", ylab = "頻度") #図19の描画

histn(p4, xlab = "P value", ylab = "頻度") #図20の描画

sum(p1 < 0.05)

## [1] 9980

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 516

sum(p4 < 0.05)

## [1] 10000

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

b1、b2、b4の検出力は100%近くなった。一方、b3の危険率は5%程度である。では、最後に変わり種の状態として、交互作用項b4だけが0ではない状況を考える。なお、以下のシミュレーションではあることを確認するためにちょっとした計算を追加している。

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

b1 <- 0

b2 <- 0

b3 <- 0

b4 <- -10

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

p5 <- NULL

p6 <- NULL

p7 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  fit3 <- anova(fit1) #分散分析表を出力

  p1 <- c(p1, coef(fit2)[1,4]) #b1のP値

  p2 <- c(p2, coef(fit2)[2,4]) #b2のP値

  p3 <- c(p3, coef(fit2)[3,4]) #b3のP値

  p4 <- c(p4, coef(fit2)[4,4]) #b4のP値

p5 <- c(p5, fit3$`Pr(>F)`[1]) #連続変数のP値

p6 <- c(p6, fit3$`Pr(>F)`[2]) #カテゴリカル変数のP値

p7 <- c(p7, fit3$`Pr(>F)`[3]) #交互作用のP値

}

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図21の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

histn(p1, xlab = "P value", ylab = "頻度") #図22の描画

histn(p2, xlab = "P value", ylab = "頻度") #図23の描画

histn(p3, xlab = "P value", ylab = "頻度") #図24の描画

histn(p4, xlab = "P value", ylab = "頻度") #図25の描画

sum(p1 < 0.05)

## [1] 522

sum(p2 < 0.05)

## [1] 514

sum(p3 < 0.05)

## [1] 487

sum(p4 < 0.05)

## [1] 10000

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

このとき、b4の検出力が100%となり、残りの係数の危険率は5%程度に抑えられた。さて、この状況は二元配置分散分析で、特定の要因組み合わせの標本だけの平均値だけが異なり、各要因単独では効果がなかった状況と類似していることは想像つくだろうか。こんな時、分散分析では交互作用だけでなく、単純主効果も有意になってしまっていた。上のシミュレーションのように、実は線形回帰でも分散分析表をanova関数に線形回帰の解析結果を入力することで計算することができる（名前に反して、anova関数は分散分析を行う関数ではない（！）、分散分析を行う関数はaov関数だ……）。fit2とfit3の中身を確認してみよう。

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

fit2

## 

## Call:

## lm(formula = y ~ x * g, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -11.5924  -3.9156  -0.6812   3.0212  13.1318 

## 

## Coefficients:

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

## (Intercept)  -3.1629     2.3181  -1.364    0.178    

## x             0.1613     0.1212   1.331    0.189    

## gB            5.1716     3.4705   1.490    0.142    

## x:gB        -10.2728     0.1846 -55.662   <2e-16 ***

## ---

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

## 

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

## Multiple R-squared:  0.9969, Adjusted R-squared:  0.9967 

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

fit3

## Analysis of Variance Table

## 

## Response: y

##           Df Sum Sq Mean Sq F value    Pr(>F)    

## x          1  72530   72530  2129.3 < 2.2e-16 ***

## g          1 427601  427601 12553.5 < 2.2e-16 ***

## x:g        1 105533  105533  3098.2 < 2.2e-16 ***

## Residuals 56   1907      34                      

## ---

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

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

fit2では唯一、交互作用項の係数だけが有意になっているが、fit3の分散分析表を見ると交互作用項だけでなく連続変数とカテゴリカル変数が有意になっている。これは矛盾しているわけではなく、計算の性質によるものだ。重回帰モデルに立ち戻ってみると、以下のような二つのパターンに式変形できる。

y ~ (b1 + b3 g) + (b2 + b4 g) x

y ~ (b1 + b2 x) + (b3 + b4 x) g

分散分析にとって「xの効果がある」というのは、xによってyの平均が変化することを意味するのだから、(b2 + b4 g) ≠ 0であればよい。つまり、例え、b2 = 0であっても、b4 ≠ 0であれば分散分析においては「xの効果がある」として検出されうる、ということだ。同様の議論でb4 ≠ 0であれば、分散分析においては「gの効果がある」と推定されうる。逆に、この考察からb2 ≠ 0かつb4 ≠ 0なのに、分散分析において「xに効果がない」とさせるようなデータを作らせることができる。つまり、(b2 + b4 g) = 0となればよいのだから、b2 = - b4 gとなるような組み合わせならよい。gは変数だが、gの変動に対して平均的に見て(b2 + b4 g) = 0となるように設定できれば良いから、平均をEで表して、b2 = - b4 E(g)を満たすデータになればよい。分散分析ならg = 0 or 1なので、E(g) = nB/(nA + nB)とする（nAは水準A、nBは水準Bのサンプルサイズ）。もちろん同様の議論で、b3 ≠ 0なのに分散分析で「gに効果がない」データも作らせることができる。

　では、逆に2つのカテゴリカル変数を説明変数とするデータに線形回帰モデルを当てはめて解析をしてみよう。二元配置分散分析の時に扱った、主効果がなく、交互作用だけが有意になりそうなデータを対象にしてみる。

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

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 70

meAY <- 50

meBX <- 50

meBY <- 50

s <- 15

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  mAX <- rnorm(nAX, mean = meAX, sd = s)

  mAY <- rnorm(nAY, mean = meAY, sd = s)

  mBX <- rnorm(nBX, mean = meBX, sd = s)

  mBY <- rnorm(nBY, mean = meBY, sd = s)

  d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )

  d$g1 <- factor(d$g1, levels = c("B","A")) #水準Bを切片にするために水準間の順位を変更

  d$g2 <- factor(d$g2, levels = c("Y","X")) #水準Yを切片にするために水準間の順位を変更

  fit1 <- lm(m ~ g1 * g2, data = d) #多重線形回帰

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  p3 <- c(p3, coef(fit2)[3,4])

  p4 <- c(p4, coef(fit2)[4,4])

}

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図26の描画

fit2

## 

## Call:

## lm(formula = m ~ g1 * g2, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -41.997  -8.578   0.662   8.999  35.066 

## 

## Coefficients:

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

## (Intercept)   46.470      2.703  17.192  < 2e-16 ***

## g1A            6.406      3.823   1.676  0.09647 .  

## g2X            2.602      3.823   0.681  0.49738    

## g1A:g2X       17.043      5.406   3.153  0.00206 ** 

## ---

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

## 

## Residual standard error: 14.8 on 116 degrees of freedom

## Multiple R-squared:  0.3309, Adjusted R-squared:  0.3136 

## F-statistic: 19.12 on 3 and 116 DF,  p-value: 3.815e-10

histn(p1, xlab = "P value", ylab = "頻度") #図27の描画

histn(p2, xlab = "P value", ylab = "頻度") #図28の描画

histn(p3, xlab = "P value", ylab = "頻度") #図29の描画

histn(p4, xlab = "P value", ylab = "頻度") #図30の描画

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 500

sum(p3 < 0.05)

## [1] 498

sum(p4 < 0.05)

## [1] 9474

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

　確かに上のシミュレーションでは、b4の検出力が95%程度で、b2、b3の危険率が5%程度となっており、予測通りの結果となった。以上の議論をまとめて判断すると、今までさんざん、aov関数を使って分散分析を行ってきたが、特別な事情がない限りは、線形回帰のlm関数を使って係数を確認したほうが、データを正確に紐解くことができる、ということだろうか。もちろん、分散分析表は一つの要約として便利なのだが、いきなり要約を確認するのではなく、その中身を見るという意味ではlm関数のほうが優秀だ。また、係数をすぐ見える形で算出してくれているので、データの予測の観点からも便利だろう。

　唯一、注意が必要な点は切片にあたるデータによっては出力結果が変わって見えることだろうか。実際は、本質は変わっていないのだが、異なる結果のように見える、ということである。例えば、上記の例では切片を標本BYにすることで、要因g1の単独効果b2 = Aに変わることの効果、要因g2の単独効果b3 = Xに変わることの効果、が無いような設定にした。つまり、E(標本BY) = E(標本AY)、E(標本BY) = E(標本BX)となっている。そして、交互作用効果b4 = AとXに代わることの効果は存在する、つまり、E(標本BY) ≠ E(標本AX)である。

　一方で、もし、切片を標本AXにしたらどうだろう。すると、b2はBに変わる効果なので、E(標本BX) - E(標本AX)として推定される。つまり、有意な効果が検出されるだろう。同様にb3も有意になる。そして、交互作用b4 = E(標本BY) - E(標本AX) - b2 - b3として推定される。ところで、データの仮定からE(標本BY) - E(標本AX) = b2 = b3なので、b4 = -b2 = -b3 ≠ 0となり、やはり有意になる。というように、切片によって結果の見え方は異なる。

多重線形回帰の危険率

　さて、危険率もシミュレートしてみよう。次のようにb1 、b2、b3、b4 = 0のデータを生成して、線形回帰を行う。

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

b1 <- 0

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1_1 <- runif(n, 0, 30)

  x1_2 <- runif(n, 0, 30)

  y1 <- b2 * x1_1 + rnorm(n = n, mean = b1, sd = s)

  y2 <- (b2+b4) * x1_2 + rnorm(n = n, mean = b1+b3, sd = s)

  d <- data.frame(x = c(x1_1, x1_2), y = c(y1, y2), g = rep(c("A","B"), each = n))

  fit1 <- lm(y ~ x * g, data = d)

  fit2 <- summary(fit1)

  fit3 <- anova(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

  p3 <- c(p3, coef(fit2)[3,4])

  p4 <- c(p4, coef(fit2)[4,4])

}

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図31の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

histn(p1, xlab = "P value", ylab = "頻度") #図32の描画

histn(p2, xlab = "P value", ylab = "頻度") #図33の描画

histn(p3, xlab = "P value", ylab = "頻度") #図34の描画

histn(p4, xlab = "P value", ylab = "頻度") #図35の描画

sum(p1 < 0.05)

## [1] 490

sum(p2 < 0.05)

## [1] 507

sum(p3 < 0.05)

## [1] 504

sum(p4 < 0.05)

## [1] 495

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

　以上のようにいずれの係数も危険率5%程度となった。

Rで行う多重線形回帰ーその1

　では、Rで多重線形回帰を行ってみよう。今回は下記のようなデータだ。さっそく、図にしてみよう。

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

x <- c(4, -9, 5, -1, 5, 11, 9, -5, 1, -8, 

       -10, 9, -4, -6, 6, 14, -3, 4, -9, 2,

       -10, -2, -4, 14, 16, -8, -8, 17, 18,

       17, 15, -10, 16, 14, -2, 7, -9, 7, -10,

       5, 9, -1, 4, 4, 7, 0, 14, 4, 0, 5, 14, 

       -9, 16, -4, 9, 4, 18, 9, -4, 11)

y <- c(6, 14, 7, 26, -9, 24, 14, 14, 46, 10, 

       7, 18, -6, 19, 9, 10, 15, 18, 2, 14, 

       -2, 4, 31, -7, 16, 29, 8, 12, 33, -4, 

       -57, 88, -53, -51, 38, -7, 75, -3, 87,

       1, -15, 14, 10, 6, -22, 27, -45, -15, 21,

       -3, -28, 78, -63, 40, -14, 14, -65, -40, 26, -17)

g <- rep(c("A","B"), each = 30)

d <- data.frame(x = x, y = y, g = g)

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T) #図36の描画

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

みると、水準Aはほとんど傾きが0で、水準Bは負の傾きを持っていそうだ。ゆえに交互作用の検出が期待される。水準Aのほうをg = 0としたときのモデル式は、

水準A：y ~ b1 + b2  x

水準B：y ~ (b1 + b3) + (b2 + b4) x

となることを考慮して係数の推定を見てみよう。

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

fit1 <- lm(y ~ x * g, data = d) #多重線形回帰

fit2 <- summary(fit1)

fit1

## 

## Call:

## lm(formula = y ~ x * g, data = d)

## 

## Coefficients:

## (Intercept)            x           gB         x:gB  

##    12.56325      0.01696     12.77297     -5.14344

fit2

## 

## Call:

## lm(formula = y ~ x * g, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -21.648  -5.697   1.419   5.983  33.420 

## 

## Coefficients:

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

## (Intercept) 12.56325    2.16176   5.812 3.07e-07 ***

## x            0.01696    0.23145   0.073 0.941842    

## gB          12.77297    3.25706   3.922 0.000242 ***

## x:gB        -5.14344    0.34668 -14.836  < 2e-16 ***

## ---

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

## 

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

## Multiple R-squared:  0.8798, Adjusted R-squared:  0.8734 

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

plotn(y ~ x, data = d, mode = "m", group = "g", pch = c(16, 17), legend = T,

      lty.leg = c(1,2)) #図37の描画

b1e <- coef(fit2)[1,1]

b2e <- coef(fit2)[2,1]

b3e <- coef(fit2)[3,1]

b4e <- coef(fit2)[4,1]

overdraw(abline(a = b1e, b = b2e, col = cols[1]), 

         abline(a = b1e+b3e, b = b2e+b4e, col = cols[2], lty = 2))

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

b1 = 12.56、b2 = 0.017、b3 = 12.77、b4 = -5.14で、b2以外は有意なようだ。つまり、水準Aの回帰直線はy = 12.56 + 0.017x、水準Bはy = 12.56 + 12.77 + (0.017 - 5.14)xになるということで、確かにデータと一致している。では、ここからfit2内に格納されている各データを計算していこう。まず、残差Residualsの分布は、

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

ye <- c(b2e * d[d$g == "A", "x"] + b1e, 

        (b2e+b4e) * d[d$g == "B", "x"] + (b1e+b3e)

        ) 

#yの推定値(水準AとBで分けて計算している)

quantile(d$y - ye) #yの実測値と推定値の差

##         0%        25%        50%        75%       100% 

## -21.648058  -5.696577   1.418660   5.982939  33.419789

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

係数Coefficientsの推定値Estimateの項は、

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

b2e_h <- cov(d[d$g == "A", "x"], d[d$g == "A", "y"])/var(d[d$g == "A", "x"])

#xの係数b2

b1e_h <- mean(d[d$g == "A", "y"]) - mean(d[d$g == "A", "x"]) * b2e_h

#(Intercept)のb1

b4e_h <- cov(d[d$g == "B", "x"], d[d$g == "B", "y"])/var(d[d$g == "B", "x"]) - b2e_h

#gBの係数b3

b3e_h <- mean(d[d$g == "B", "y"]) - mean(d[d$g == "B", "x"]) * (b4e_h + b2e_h) - b1e_h

#x:gBの交互作用係数b4

b1e_h

## [1] 12.56325

b2e_h

## [1] 0.0169617

b3e_h

## [1] 12.77297

b4e_h

## [1] -5.143441

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

係数Coefficientsの標準誤差Std.Errorの項は、

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

se <- sum((d$y - ye)^2)/(length(x) - 4)

#残差の不偏分散、自由度はサンプルサイズ- 推定するパラメータ数なので、60 - 4 = 56

sb1 <- sqrt(

  se * (1/length(d[d$g == "A", "x"]) + 

          (mean(d[d$g == "A", "x"])^2)/(var(d[d$g == "A", "x"]) * (length(d[d$g == "A", "x"])-1)

                       )

        )

  ) #b1 standard error

sb2 <- sqrt(se/(var(d[d$g == "A", "x"]) * (length(d[d$g == "A", "x"])-1))) #b2 standard error

sb3 <- sqrt(

  se * (1/length(d[d$g == "B", "x"]) + 

          (mean(d[d$g == "B", "x"])^2)/(var(d[d$g == "B", "x"]) * (length(d[d$g == "B", "x"])-1)

                       )

        )

  + sb1^2) #b3 standard error

sb4 <- sqrt(se/(var(d[d$g == "B", "x"]) * (length(d[d$g == "B", "x"])-1)) + sb2^2) #b4 standard error

sb1

## [1] 2.161762

sb2

## [1] 0.231455

sb3

## [1] 3.257065

sb4

## [1] 0.3466769

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

係数Coefficientsのt valueの項は、

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

b1t <- b1e_h/sb1 #b1 t-value

b2t <- b2e_h/sb2 #b2 t-value

b3t <- b3e_h/sb3 #b3 t-value

b4t <- b4e_h/sb4 #b4 t-value

b1t

## [1] 5.81158

b2t

## [1] 0.07328293

b3t

## [1] 3.92162

b4t

## [1] -14.83641

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

係数CoefficientsのPr(>|t|)の項は、

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

pt(q = b1t, df = length(x)-4, lower.tail = sign(b1t) < 0) * 2

## [1] 3.072685e-07

pt(q = b2t, df = length(x)-4, lower.tail = sign(b2t) < 0) * 2

## [1] 0.9418421

pt(q = b3t, df = length(x)-4, lower.tail = sign(b3t) < 0) * 2

## [1] 0.0002423744

pt(q = b4t, df = length(x)-4, lower.tail = sign(b4t) < 0) * 2

## [1] 4.696814e-21

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

残差の標準誤差Residual standard errorは

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

sqrt(se) #Residual standard error

## [1] 11.51745

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

R二乗値は

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

var(ye)/var(y) #Multiple R-squared

## [1] 0.8798313

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

F値とその時のP値は

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

sr <- sum((ye - mean(ye))^2)/3 #パラメータの自由度で割る

sr/se #F-value

## [1] 136.6705

pf(df1 = 3, df2 = length(x) - 4, q = sr/se, lower.tail = F)

## [1] 9.753464e-26

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

というように計算される。以上のように、2水準に含まれる連続変数に関して、回帰直線を推定することができた。これで条件によって異なる連続変数の応答を推定することが可能となった。