連続変数を複数含む解析

　では、回帰分析の最後に連続変数を複数含んだ解析を紹介する。一般にはこのような解析を多重線形回帰ということが多いと思うが、これまでのページでは変数がカテゴリカルか連続かに問わず、複数の説明変数を含んだ線形回帰を多重線形回帰と紹介してきた。モデルの構造を見れば、そう考えても差し支えないことがすぐにわかるだろう。x1とx2を連続変数とすれば、

y ~ b1 + b2 × x1 + b3 × x2 + b4 (x1 × x2)

となり、連続×カテゴリカルの多重線形回帰と構造は大きく変わらない。このモデル式は次の2通りに変形ができる。

y ~ (b1 + b3 × x2) + (b2 + b4 × x2) x1……(1)

y ~ (b1 + b2 × x1) + (b3 + b4 × x1) x2……(2)

すなわち、(1)はx1-y平面で見たときにx2の値によって切片と傾きが変化する直線群、(2)はx2-y平面で見たときにx1の値によって切片と傾きが変化する直線群ととらえることが可能である。そして、x1やx2は連続変数なので、x1-x2-y空間内では徐々に切片や傾きが変化する直線が連続的につながった曲面を表している、と考えられる（下図）。しいて言うなら、この辺りがカテゴリカル変数を説明変数として持つ場合との違いであり、カテゴリカル変数が説明変数とした場合、各水準で独立な回帰直線を求めていた。

今回の帰無仮説は今までと同様に、残差が平均0の正規分布に従うと仮定し、「b1 = 0、b2 = 0、b3 = 0、b4 = 0」である。

係数の推定

　これまでの線形回帰と同様に最小二乗法を用いて係数を推定することができるが、いよいよこの辺りから、(係数) = (単純な数式)という表記が難しくなってくる。まず今回は、2つの連続変数とそれらの交互作用項を持つモデルに関して、最小二乗法で解く。そのあと、交互作用項がないときは2連続変数の係数がより簡単に表記できることを示す。そして、追記でより一般に複数の連続変数が存在し、交互作用がないときの係数の計算について紹介する。

このように推定する係数が3以上になると、一つ一つの偏微分の条件は簡単でも、解くべき連立方程式は複雑になり、簡単な式で示すことは困難になる。それゆえ、簡単のために最後は行列とベクトルを使って示した。推定する係数にかかるE(x)などの係数をまとめた行列Xについて、逆行列X^-1を求めることができれば、係数を推定することができる（ただし、ここではひたすら純朴に計算したが、行列Xの形を工夫することでもっと簡単な形まで落とし込める、追記にて記す）。Rで計算するうえでは、逆行列まで求める必要はなく、Xb = yの形から、推定値ベクトルbを計算してくれる。

　では、交互作用がない場合はどうなるだろうか。b4 = 0として考えればよいので、式(1)～(3)は以下の通りになり、さらに各係数まで計算する。

このように相関係数をもちいて比較的簡単な形にすることができる。上記のように、多重線形回帰によって求められた係数は偏回帰係数partial regression coefficientと呼ばれる。

多重線形回帰の検出力

　さて、今回もシミュレートしてみよう。連続変数×カテゴリカル変数の時とほとんど同じような感じになるのは想像に難くない。とはいえ、データの構造はやや異なるので、どんなデータが生成されているかをモデルから考えることは、統計の直感を養う上で重要だ。

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

library(plotn)

library(rgl) #作図に必要なので、ない場合はインストール

library(reshape2) #作図に必要なので、ない場合はインストール

#以下、作図前のおまじない

pp <- par3d(no.readonly = T)

pp$listeners <- 473

pp$userMatrix <- matrix(c(0.5437224, -0.8392414, 0.00631617, 0, 

                          0.2854379, 0.1919944, 0.93896931, 0, 

                          -0.7892347, -0.5087356, 0.34394285, 

                          0, 0, 0, 0, 1), nrow = 4,byrow = T)

pp$windowRect <- c(190, 213, 446, 469)

close3d()

#以上、作図前のおまじない

b1 <- 10

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 60

x1 <- runif(n, 0, 30)

x2 <- runif(n, 0, 30)

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 14.626422 10.628340  5.997423

## 2  1.703273  7.839096 15.706092

## 3 15.573617 20.087065  9.074815

## 4  9.704097 15.745849  5.306351

## 5 15.455977 11.738793  9.423308

## 6 21.048418 28.390182  9.771883

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

fit2 <- summary(fit1)

#以下、rglパッケージを使ったときの3D作図

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図1の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_1.png")

close3d()

#以上、rglパッケージを使ったときの3D作図

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

このように2連続変数によるデータ生成は3次元的広がりを見せる。この広がりを

y ~ b1 + b2 × x1 + b3 × x2 + b4 (x1 × x2)

とあらわされる曲面モデルで回帰しているのである。3次元空間を1方向から見て理解するのは難しい。ぜひとも、上のスクリプトでpar3d()を打たずに（これをした後は図が固定される）、いきなりpoints3d()～aspect3d()を打った後、ぐりぐりと図を動かしてみるとよい。rglパッケージで作図されたものは、コマンドライン以外にも直接図をドラッグなどして操作できる。ほかには、上記のような3次元図ではなく、2次元に要約した下記のような作図もあり得る。

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

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")#図2の作図

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")#図3の作図

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

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

図2と図3はそれぞれ、x1-y平面やx2-y平面から全データ点を見た図である。そして、x2やx1の値の大きさによってlow-medium-highと色分けした。全データのうち小さいものから33%までのデータがlow、次の33%がmedium、そして大きい方の33%がhighとなっている。また、図2を例に説明すると回帰直線は、x2のデータのうち小さい方から約16.65%（つまり、0%と33.3%の中央値）近辺の回帰直線をlow、50%（全体の中央値）近辺のものをmedium、83.3%（つまり、66.6%と100%の中央値）近辺のものをhighとしている。lowからhighにかけて、切片と傾きが連続的に変化しているが、そのなかのスナップショットを図に示している感じである。

　データの生成過程から、切片b1だけが0から異なっており、b2やb3のような傾きが存在しない。このようなデータ生成に沿ったものができていることが、図1～3からも見て取れるだろう。このような時の検出力を以下のように計算する。

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

b1 <- 10

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 60

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1 <- runif(n, 0, 30)

  x2 <- runif(n, 0, 30)

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

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

  fit1 <- lm(y ~ x1 * x2, 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 = "頻度")#図4の作図

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

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

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

sum(p1 < 0.05)

## [1] 9407

sum(p2 < 0.05)

## [1] 511

sum(p3 < 0.05)

## [1] 525

sum(p4 < 0.05)

## [1] 499

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

以上のように切片b1の検出力は94%で、残りの係数の危険率は5%付近に抑えられている。

　次は一方の連続変数に傾きを追加してみよう。

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

b1 <- 10

b2 <- 2

b3 <- 0

b4 <- 0

s <- 5

n <- 60

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1 <- runif(n, 0, 30)

  x2 <- runif(n, 0, 30)

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

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

  fit1 <- lm(y ~ x1 * x2, 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値

}

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図8の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_2.png")

close3d()

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

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

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")#図9の作図

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")#図10の作図

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

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

図8～11を見ると、x1軸方向には傾斜になっていて、x2軸方向はほとんど傾きがないことがわかる。交互作用がほとんどないので、x1軸方向の傾斜はx2による変化がほとんどない。

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

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

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

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

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

sum(p1 < 0.05)

## [1] 9431

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 524

sum(p4 < 0.05)

## [1] 537

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

以上のようにb1とb2の検出力は100%に近く、残りの係数の危険率は5%付近に抑えられている。

　次は両方の連続変数に傾きを追加してみよう。

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

b1 <- 10

b2 <- 2

b3 <- 5

b4 <- 0

s <- 5

n <- 60

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1 <- runif(n, 0, 30)

  x2 <- runif(n, 0, 30)

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

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

  fit1 <- lm(y ~ x1 * x2, 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値

}

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図15の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_3.png")

close3d()

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

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

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")#図16の作図

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")#図17の作図

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

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

図15～17を見ると、x1軸、x2軸の両方向に傾斜になっている。交互作用がほとんどないので、x1軸方向の傾斜もx2軸方向の傾斜も、他方の連続変数による変化がほとんどない。

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

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

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

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

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

sum(p1 < 0.05)

## [1] 9422

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 10000

sum(p4 < 0.05)

## [1] 469

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

以上のようにb1、b2、b3の検出力は100%に近く、b4の危険率は5%付近に抑えられている。

　最後に、b3 = 0として、交互作用を追加してみよう。

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

b1 <- 10

b2 <- 2

b3 <- 0

b4 <- -3

s <- 5

n <- 60

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1 <- runif(n, 0, 30)

  x2 <- runif(n, 0, 30)

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

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

  fit1 <- lm(y ~ x1 * x2, 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値

}

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図22の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_4.png")

close3d()

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

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

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")#図23の作図

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")#図24の作図

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

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

図22～24を見ると、x1やx2が大きな値を示すようになると、負の値の傾斜が強くなるさまが見て取れる。図22で言うならば、手前側は平地となっているが、奥側はくぼんでいる、という感じである。これはデータの生成過程でb4に負の値を入れているためである。

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

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

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

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

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

sum(p1 < 0.05)

## [1] 9360

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 487

sum(p4 < 0.05)

## [1] 10000

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

以上のようにb1、b2、b4の検出力は100%に近く、b3の危険率は5%付近に抑えられている。

多重線形回帰の検出力

　では、すべての係数が0のとき、危険率を調査してみよう。

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

b1 <- 0

b2 <- 0

b3 <- 0

b4 <- 0

s <- 5

n <- 60

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:10000){

  x1 <- runif(n, 0, 30)

  x2 <- runif(n, 0, 30)

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

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

  fit1 <- lm(y ~ x1 * x2, 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値

}

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図29の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_5.png")

close3d()

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

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

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")#図30の作図

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")#図31の作図

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

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

図29～31をみると、y = 0近辺のほぼ水平な面が推定されていることがわかる。

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

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] 9360

sum(p2 < 0.05)

## [1] 10000

sum(p3 < 0.05)

## [1] 487

sum(p4 < 0.05)

## [1] 10000

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

そして、予想通り各係数の危険率は5%付近に抑えられている。

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

　では今回も実際の解析例を示そう。下記のようなデータを解析する。まずは作図してみよう。

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

x1 <- c(1, 18, 5, 16, 10, 0, 13, 1, 6, 5, 

        19, -10, 1, 1, -9, 7, 2, 6, 11, -10, 

        11, 15, 12, 3, 8, -6, 18, 15, 16, -1,

        10, 15, 6, 4, 4, 12, 2, -6, -7, -5, 14, 

        -3, 3, -8, 5, 12, -7, 10, 0, 12, 18, -4,

        0, 2, -7, 18, -8, 3, 17, -4)

x2 <- c(-8, 8, 5, -8, 7, -4, 2, -6, 3, -5, 12, 9, 

        -10, 8, 2, 9, 2, 2, 7, -7, -6, 18, 19, -9,

        6, -2, 18, 9, -4, -4, 8, 9, 8, 11, 14, 10,

        2, -1, -5, 4, 17, -9, 17, -2, 16, 0, 0, 17, 

        15, 7, 3, 17, -5, 15, -5, -7, -1, -6, -8, 5)

y <- c(4, 930, 187, -418, 442, 10, 259, -4, 149, -52, 1328, 

       -606, -7, 40, -202, 381, 32, 115, 468, 280, -147, 1467, 

       1170, -56, 303, -10, 1790, 805, -98, 7, 488, 828, 303, 

       232, 268, 700, 39, -53, 114, -193, 1301, 134, 281, -4, 

       406, 146, -100, 868, -43, 535, 475, -437, 8, 120, 96, 

       -384, -53, -35, -439, -148)

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

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black")#図36の作図

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_6.png")

close3d()

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

さすがに今回に関しては図36だけを見てデータを予測するのは難しいので、ぜひ、各自でrglの図をぐりぐり動かしてほしい。するとわかるのが、x1とx2が大きくなるとyが大きくなり、x1とx2が中間的だとyの変動がほぼなくなり、x1とx2が小さくなるとyが大きくなる、というような関係にあることがわかる。図36で言えば、奥と手前が上にせり出しており、左端と右端は下にせり出している、という感じだ。では、このようなデータを多重線形回帰してみよう。

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

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

fit2 <- summary(fit1)

fit1

## 

## Call:

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

## 

## Coefficients:

## (Intercept)           x1           x2        x1:x2  

##      0.4766      12.1310      -3.3761       4.9573

fit2

## 

## Call:

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

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -42.156  -7.526  -1.626   7.197  48.704 

## 

## Coefficients:

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

## (Intercept)  0.47659    2.42166   0.197    0.845    

## x1          12.13100    0.26168  46.359  < 2e-16 ***

## x2          -3.37611    0.31736 -10.638 4.63e-15 ***

## x1:x2        4.95725    0.03097 160.089  < 2e-16 ***

## ---

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

## 

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

## Multiple R-squared:  0.999,  Adjusted R-squared:  0.9989 

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

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

特に変わったことをすることなく、今まで通り、lm関数を使って解析可能だ。Coefficients項を見ると、x1およびx2を変数と見たときの回帰直線は次のようになる。

y = (0.48 - 3.38 x2) + (12.13 + 4.96 x2) x1……(1)

y = (0.48 + 12.13 x1) + (-3.38 + 4.96 x1) x2……(2)

(1)はx2が増加すると傾きが増加する直線、(2)はx1が増加すると傾きが増加する直線を示している。つまり、x1とx2が正に大きくなれば

y = (負の値) + (正の値) x1……(1)

y = (正の値)  + (正の値) x2……(2)

となり、　yが正に大きな値を示すことが予測される。一方で、負に大きくなれば、

y = (正の値) + (負の値) x1……(1)

y = (負の値)  + (負の値) x2……(2)

となり、x1とx2は負の値なのだから、やはりyの値が正に大きくなることが予測される。したがって、実際のデータに合致した予測が得られているだろう。回帰曲面を作図すると以下のようになる。

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

d$pred <- predict(fit1)

x1_pred <- seq(min(d$x1),max(d$x1),len=100)

x2_pred <- seq(min(d$x2),max(d$x2),len=100)

plot_pred <- expand.grid(x1 = x1_pred, x2 = x2_pred)

plot_pred$y <- predict(fit1, newdata = plot_pred)

y_pred <- dcast(plot_pred, x1 ~ x2, value.var="y")[-1]

par3d(pp)

par3d(windowRect = c(20, 30, 800, 600))

points3d(d$x1, d$x2, d$y, col="black") #図37の作図

surface3d(x1_pred, x2_pred, as.matrix(y_pred),col="gray",alpha = 0.2)

axes3d()

title3d(xlab="x1", ylab="x2", zlab="y")

aspect3d(1,1,0.6)

rgl.snapshot(fmt="png", "plot3d_7.png")

close3d()

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

そして、

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

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") #図38の作図

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") #図39の作図

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

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

となり、予測通りになることが確認できるだろう。では、各項の計算をしてみよう。残差Resiguals項の四分位数は、以下のように計算できる。

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

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

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

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

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

ye2 <- c(b4e * d$x1 * d$x2 + b3e * d$x2 + b2e * d$x1 + b1e) #yの推定値

quantile(d$y - ye2)

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

## -42.156279  -7.525634  -1.625622   7.197243  48.704296

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

次にCoefficients項の各係数Estimateだが、以下のように計算できる。非常にややこしい式になっている。自分で定義したFunc関数は共分散公式を利用して、積の平均を共分散と平均の積の和から計算する関数である。また、solve関数は行列XとベクトルYを引数とし、solve(X, Y)とすれば、Y = Xbを満たすベクトルbを計算してくれる。

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

Scov <- function(x, y) cov(x, y) * (length(x)-1)/length(x) #標本共分散の計算

Func <- function(x, y = NULL) {

  if(is.null(y)) y <- x

  Scov(x, y) + mean(x) * mean(y)} #積の平均の計算

X <- matrix(c(1, mean(d$x1), mean(d$x2), Func(d$x1, d$x2),

              mean(d$x1), Func(d$x1), Func(d$x1, d$x2), Func(d$x1^2, d$x2),

              mean(d$x2), Func(d$x1, d$x2), Func(d$x2), Func(d$x1, d$x2^2),

              Func(d$x1, d$x2), Func(d$x1^2, d$x2), Func(d$x1, d$x2^2), Func(d$x1^2, d$x2^2)),

  4, 4)

Y <- c(mean(d$y), Func(d$x1,d$y), Func(d$x2, d$y), Func(d$x1 * d$x2, d$y))

b <- solve(X, Y)

b

## [1]  0.4765948 12.1309962 -3.3761055  4.9572532

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

そして標準誤差Std.Errorを計算しよう。ここで計算方法の詳細を紹介することはしないが、追記に書かれた方法から、行列Xにサンプルサイズを掛けたnXの逆行列に残差不偏分散を掛けたものがパラメータの分散共分散行列となることを利用して、以下のように解ける。なお、solve関数は行列だけを引数とすると、逆行列を計算してくれる。

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

se <- sum((d$y - ye2)^2)/(length(d$x1) - 4)

#自由度はサンプルサイズ60で、推定するパラメータが4つなので、60 - 4 = 56

bS <- sqrt(se * solve(nrow(d)*X))

## Warning in sqrt(se * solve(nrow(d) * X)): 計算結果が NaN になりました

bse <- c(bS[1,1],bS[2,2],bS[3,3],bS[4,4])

bse

## [1] 2.42165900 0.26167583 0.31735800 0.03096553

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

係数と標準誤差が求まったのでt値は、

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

t <- b/bse

t

## [1]   0.1968051  46.3588708 -10.6381610 160.0894010

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

そして、このt値に対応するP値は、

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

pt(q = t[1], df = nrow(d)-4, lower.tail = sign(t[1]) < 0) * 2

## [1] 0.8446926

pt(q = t[2], df = nrow(d)-4, lower.tail = sign(t[2]) < 0) * 2

## [1] 2.313786e-46

pt(q = t[3], df = nrow(d)-4, lower.tail = sign(t[3]) < 0) * 2

## [1] 4.631783e-15

pt(q = t[4], df = nrow(d)-4, lower.tail = sign(t[4]) < 0) * 2

## [1] 3.1973e-76

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

残差の不偏分散の平方根Residual standard errorは、

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

sqrt(se) #Residual standard error

## [1] 15.77042

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

R二乗値は、

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

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

## [1] 0.998986

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

F値とその時のP値は

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

sr <- sum((ye2 - mean(ye2))^2)/3

sr/se #F-value

## [1] 18389.84

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

## [1] 8.931333e-84

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

以上のように、一通りの数値を計算できた。これにて、回帰分析の根本はほとんど理解できた、と考えてもよいだろう。とはいえ、実際のデータ解析にはまだまだ、悩ましい問題が山積だ。上記の例を見ても、たった2つ連続変数において、交互作用を考慮すると非常にややこしくなってしまった。ちなみに、例えば3変数の交互作用を考慮したモデルは下記のようになる。

y ~ b1 + b2 x1 + b3 x2 + b4 x3 + b5 x1 x2 + b6 x2 x3 + b7 x3 x1 + b8 x1 x2 x3

全ての交互作用を考慮する場合、なんと、推定するべき係数の数は2^(変数の数)と、爆発的に増える。こんなに交互作用があるとデータの挙動が複雑すぎて、人間の理解の限度を軽く超えてしまう。また、推定するべき係数が増えるということは、その分だけデータを増強する必要があるが、指数増殖するパラメータを高い精度で予測できるデータをそろえられることはまずないだろう。交互作用を考えることでデータの解釈の幅は広がるものの、逆にモデルの複雑化を招く厄介者でもあるのだ。もし、もっと多数の変数を扱ったら……、と思うとぞっとするだろう。それゆえか、多数の変数を扱う場合は交互作用を考慮しないことも多い。それにモデルが複雑になればなるほど、予測精度が落ちることが一般に知られている。もちろん、交互作用に限らず、単純主効果もあればあるだけいいモデルになるわけではない。もし、主効果間で極めて高い相関が存在する場合、解析が失敗することがある。どの変数をモデルに組み込むか、は実のところかなり難しい問題で、かつデータを見定めるセンスがいる部分なのだ……。とりあえず、私なりの1つの指針は、多変数を扱うことの意義（メタ解析や気象データなど）がない場合、特に実験操作による応答を見る場合は、可能な限り最大でも2変数の解析に落とし込むのが良いと思う。3変数存在するときは、どれかカテゴリカル変数ならそれを基準に分けて残りの2変数の応答を見るし、カテゴリカル変数がない場合も便宜的にカテゴリカルにして残りの2変数の応答を見る。これは一つの方法に過ぎず、正解というわけではない。ともかくも解析は手段であって、データの見せ方はより分かりやすくできるように工夫してみるのが良いだろう。

追記1：交互作用がないときの係数の推定

　以下のように、行列を用いて比較的簡単な表記が存在することを示す。

なお、補足であるが、推定したいパラメータ群^bがうまく推定されるためには、式の様子から、X^TXに逆行列が存在しなければならない。逆行列が存在しない条件は、変数間に極めて強い相関が存在する=多重共線性がある場合である。

　また、分散は次のように表記できる。最終的に求まるものは分散共分散行列であり、その行列の対角成分が対応するパラメータの分散である。