目次

・線形モデルの前提の再考

・分散が説明変数に依存している

・Box-Cox変換

　・λの決定方法―最尤法―

　・Box-Cox変換後の回帰分析

・補足：Box-Cox変換の連続性

・補足：Yeo-Johnson変換

線形モデルの前提の再考

　これまで様々なパターンの線形モデルを取り扱ってきた。今ここで、検定の時と同じように、線形モデルを使える前提について考えてみよう。線形モデルは、（1）残差の分布が正規分布に従う、（2）残差の分散は説明変数の値によらず一定、（3）説明変数と被説明変数が直線関係（線形）であることを仮定している。これらの仮定が重要であることは、線形回帰にて最小二乗法を使って係数の推定を行っていれば、おのずと理解できる。ゆえに本来は、これらのいずれかの仮定が満たされないときに単純な線形モデルに当てはめて解析するのは問題がある。（3）に関しては以前、非線形最小二乗法でも取り扱った。あるいは、今後、一般化線形モデルでも解決する術を紹介することになるだろう。今回は、（1）や（2）も絡んでくる状況を考えよう。

分散が説明変数に依存している

　重み付き最小二乗法の時は、ほとんどのデータは1直線状なのに、外れ値があることで回帰結果が大きくゆがむことを紹介した。これも残差の分散が一定ではない例であるが、今回は分散が明らかに説明変数に依存している状況を紹介しよう。

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

library(plotn)

library(car)

x1 <- c(10.1, 9.43, 2.34, 5.28, 17.32, 15.13, 12.68, 5.12, 28.99, 18.44, 

        19.11, 28.44, 18.43, 16.87, 23.41, 2.61, 13.98, 6.4, 8.7, 19.95, 

        10.08, 26.24, 22.27, 0.44, 20.61, 10.97, 15.66, 17.43, 14.97, 12.84,

        19.74, 24.69, 12.98, 20.48, 28.01, 13.69, 10.3, 0.56, 7.33, 23.6, 

        17.48, 0.15, 7.74, 29.8, 1.3, 27.84, 9.56, 7.15, 24.71, 22.13, 

        0.34, 4.58, 24.01, 11.33, 24.97, 17.72, 21.25, 14.61, 4.07, 17.79, 

        23.39, 21.75, 3.37, 5.96, 2.34, 24.74, 16.73, 17.63, 29.01, 8.46, 

        14.68, 25.35, 15.58, 5.89, 7.61, 20.34, 17.32, 18.85, 19.1, 15.87,

        29.98, 19.49, 14.78, 0.72, 27.62, 9.79, 0.28, 25.55, 25.38, 16.17,

        9.92, 6.06, 26.76, 21, 17.59, 8.69, 7.26, 18.79, 21.8, 6.86)

y <- c(11.92, 6.98, 0.79, 5.88, 7.18, 4.75, 8.22, 4.72, 17.88, 6.9, 16.99, 

       15.91, 7.1, 4.83, 8.92, 1.08, 5.12, 4.26, 4.52, 13.26, 7.39, 15.56, 

       15.84, 0.58, 16.05, 5.1, 6.22, 10.4, 4.33, 8.65, 5.68, 14.35, 14.26,

       12.98, 11.46, 7.19, 4.19, 1.06, 3.38, 19.45, 6.38, 0.82, 3.05, 13.4,

       2.73, 8.92, 7.77, 6.72, 17.08, 11.69, 1.23, 4.8, 12.34, 7.87, 19.82,

       10.75, 6.41, 12.21, 5.81, 11.38, 11.06, 9.51, 3.95, 4.05, 1.82, 8.06,

       11.45, 7.52, 28.43, 3.46, 16.03, 6.63, 4.77, 3.29, 5.81, 7.12, 12.08,

       18.49, 13.86, 8.23, 18.59, 9.66, 4.59, 0.87, 14.55, 5.91, 1.11, 8.23, 

       12.72, 8.75, 5.7, 4.06, 10.04, 12.98, 7.55, 4.43, 3.02, 3.69, 9.69, 5.56)

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

head(d)

##      x1     y

## 1 10.10 11.92

## 2  9.43  6.98

## 3  2.34  0.79

## 4  5.28  5.88

## 5 17.32  7.18

## 6 15.13  4.75

plotn(y ~ x1, data = d)#図1の描画

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

上記のデータは明らかにx1が大きくなるほど、分散が大きくなっている。こんなときも、一般的な線形回帰の仮定にはそぐわない。これによってどういう問題が起こりえるかは、難しい問題な気がするが、あり得るのはx1が大きくなると、外れ値の出現が起こりやすくなるだろう。そうなると、外れ値によってロバストな解析が阻害されると予測される。

Box-Cox変換

　今回は、データを変換することで線形回帰の仮定を満たすような状況に落とし込むことを目指す。今回のデータはx1が大きくなるにつれて分散が大きくなっている。このような状況の時、経験的には対数をとると分散が均一になることが知られている。

　より一般に、yi > 0において、以下のような変換をすることで、残差を分散が均等な正規分布に近づけられる。この変換をBox-Cox変換あるいは冪正規変換と呼ぶ（負の値は変換できないので、すべてのデータが正になるように調整するか、補足のYeo-Johnson変換を用いる）。Box-Cox変換は、累乗変換、平方根変換、対数変換を統合した変換である。

変換のされ方はλに依存する。後述するように、データを最も正規分布に近づけられるλを最尤法で求める（最尤法は当該ページ参照）。この変換はλに応じて、以下のような曲線を描く。

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

BoxCox <- function(x, lambda){

  mapply(function(tmp1, tmp2){

    if(tmp2 == 0){

      log(tmp1)

    } else {

      (tmp1^tmp2 - 1)/tmp2

    }

  }, x, lambda)

}#Box-Cox変換

x <- seq(0.0001, 2, length = 200)

d <- data.frame(x = x, 

                y1 = BoxCox(x, lambda = -1),

                y2 = BoxCox(x, lambda = -0.5),

                y3 = BoxCox(x, lambda = 0),

                y4 = BoxCox(x, lambda = 0.5),

                y5 = BoxCox(x, lambda = 1),

                y6 = BoxCox(x, lambda = 2),

                y7 = BoxCox(x, lambda = 3))

plotn(d[,1], d[,2:8], mode = "m", type = "l", ylim = c(-2,2),

      xlab = "y", ylab = "transformed y",

      legend = T, leg.title = "lambda", leg.lab = c(-1,-0.5,0,0.5,1,2,3),

      pch.leg = NA, line = T, leg.sp = 4)#図2の描画

overdraw(abline(v = 1, lty = 2),

         abline(h = 0, lty = 2))

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

λ > 1のとき変換は下に凸な関数となり、λ < 1のとき上に凸な関数となる。このことから、より具体的に変換の様子を確認すると以下のような感じである。

変換後の値が均一になるように線を引き、それが変換前にどんな値だったかを見ると、変換が下に凸であれば右の値の密度が高く、上に凸だと左の値の密度が高くなる。つまり、Box-Cox変換は、λ > 1であれば左に尾を引く分布を、λ < 1であれば右に尾を引く分布を左右対称にする効果を持つ。

　さらに、もし変換後の値が正規分布に従っているなら、同様に変換前は λ > 1であれば左に尾を引く分布を、λ < 1であれば右に尾を引く分布となることがシミュレーションからわかる。言い換えれば、Box-Cox変換において適切にλを指定することで、左右対称ではない分布を正規分布に近づけることが可能である。

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

Inv_BoxCox <- function(x, lambda){

  mapply(function(tmp1, tmp2){

    if(tmp2 == 0){

      exp(tmp1)

    } else {

      (tmp2 * tmp1 + 1)^(1/tmp2)

    }

  }, x, lambda)

}#Box-Cox変換の逆変換

histn(Inv_BoxCox(rnorm(10000, mean = 2, sd = 1), lambda = 0.3))#図3の描画

histn(Inv_BoxCox(rnorm(10000, mean = 2, sd = 1), lambda = 4))#図4の描画

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

・λの決定方法―最尤法―

　では、具体的にどのようにしてλを求めるかというと、最尤法という方法を用いる。様々なλで変換したデータが正規分布から得られる「確率（本当は確率ではないがそれに相当するもの、と考えればよい）」を計算し、その「確率」が最も大きいλを採用する方法である。ここでいう「確率」は尤度likelihoodと呼ばれる。

　Box-Cox変換時の尤度関数は以下の通りである。

yiの変換後をyi^(λ)とし、その平均を̄yi^(λ)、分散をσ^2とすれば、{}は変換後のデータは正規分布に従うと仮定した時の尤度となる。それだけでなく、データを変換したときに、データ間の距離は変換によって変化するので、その変化分を考慮しなければならない。それがヤコビ行列式、ヤコビアンJacobian determinantと呼ばれる係数であり、偏微分係数を用いる。ちょうど、積分で変数変換したときに登場する調整項みたいなものである。

　上記を対数をとって整理すると以下の通り。

この対数尤度を最大にするλが、最もデータを正規分布に近づけるパラメータであると判断するわけである。

　実際にRで解析してみよう。carパッケージのpowerTransform関数を用いると、適切なλを自動で計算してくれる。ここで求めたλが尤度を最大化する値であることを、自分で尤度関数を定義することで確認してみよう。

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

powerTransform(d$y)#最尤推定されたλ

## Estimated transformation parameter 

##       d$y 

## 0.4755349

lL <- NULL

for(i in seq(-2, 2, length = 200)){

  y_t <- BoxCox(y, lambda = i)

  n_t <- length(y_t)

  sy_t <- sqrt(var(y_t) * (n_t - 1)/n_t)

  lL_t <- sum(log(y^(i-1)/sy_t))

  lL <- c(lL, lL_t)

}

plotn(seq(-2, 2, length = 200), lL, type = "l")#図5の描画

overdraw(abline(v = powerTransform(d$y)$lambda, col = "red"))

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

確かに、推定されたλは尤度関数を最大化している。このλを使って、bcPower関数によってBox-Cox変換してみよう（自分で定義したBoxCox関数でもよいが）。

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

y_t <- bcPower(d$y, powerTransform(d$y)$lambda)#Box-Cox変換

plotn(y_t ~ x1)#図6の描画

histn(y)#図7の描画

histn(y_t)#図8の描画

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

すると、変換後はxによらずデータの分散が均質化されている。また、データ全体の分布は右に尾を引く形から、左右対称な分布に変わっている。実際に、正規分布に近づいているかどうかを累積分布関数を使って確認してみよう。

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

cdf <- function(p, v){

  sapply(v, function(x) sum(p < x))/length(p)

}

m <- mean(y)

s <- sd(y)

ry <- seq(min(y)-5, max(y), length = 200)

cy <- cdf(p = y, v = ry)

plotn(ry, cy, type = "l", xlab = "y", ylab = "CDF")#図9の描画

overdraw(points(ry, pnorm(ry, mean = m, sd = s), type = "l", col = "red"))

m <- mean(y_t)

s <- sd(y_t)

ry <- seq(min(y_t)-2, max(y_t), length = 200)

cy <- cdf(p = y_t, v = ry)

plotn(ry, cy, type = "l", xlab = "y", ylab = "CDF")#図10の描画

overdraw(points(ry, pnorm(ry, mean = m, sd = s), type = "l", col = "red"))

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

すると、データの歪みが変換によって改善し、正規分布に近づいていることがわかる。

　実はこの変換は説明変数が正規分布に従っていることが重要である。この仮定を聞いて、私は過去に抱いたBox-Cox変換に対する不信感を払拭できた。というのも、以下のような疑問を抱いたからだ。xによって分散が変わるということは、x各点におけるyの分布が違っている可能性があるということだ。それならば、xに依存して変化するλが必要なのではないか、そうしないとx各点におけるyの分布が等分散の正規分布になるはずがないと思ったのだ。

実際には、これまで示したようにxの各点を考慮するのではなくy全体の分布を1つのλで一括変換するだけでよい。

さてこの図を見ると、この変換がうまくいくかは、結局、Box-Cox変換後のデータ全体が正規分布になれるか否が重要であることがわかる。

読み飛ばし可******************************************

　ちょっと難しい話をする。分布の混合がわかればすぐにわかることだが、イメージを大切に説明を尽くす。平均値が連続的に変化する正規分布の混合分布を考えてみよう。例えば、分散が一定で、平均-1の正規分布の確率密度を考える。0.01刻みとして、次に-0.99、-0.98……、1まですべての正規分布を考えて、全ての確率密度を平均する。この刻みをもっと小さくして、さらに-∞から+∞まで考えることとする。さらに、この平均値自体が正規分布に従って生成される状況を考える。このとき、この混合分布は正規分布になることが知られる。上の右図であれば、xが正規分布に従っているとして、各xの値に正規分布が対応して、yi^(λ)が生成されてるとする。そのすべての正規分布の確率密度の平均は、全体で大きな正規分布＝yi^(λ)全体の分布が正規分布になるということなのだ。全体が大きな正規分布になるためには、xが正規分布に従っている必要がある。

*********************************************************

xが正規分布に従っていれば、yi^(λ)全体の分布が正規分布になることができる。そして、それに合わせたλを1つだけ決定するだけでよいわけである。そして、yi^(λ)全体の分布が正規分布なら、x各点のyi^(λ)も正規分布に従っている。正規分布に従ってる各点のyi^(λ)を逆変換したものが元のデータだ。

・Box-Cox変換後の回帰分析

　さて、Box-Cox変換をすることで、線形回帰を利用できる形になった。こうなればシンプルで、変換後データに対してlm関数を使って線形回帰を行う。

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

fit1 <- lm(y_t ~ x1)#変換後データを使った回帰分析

plotn(y_t ~ x1)#図11の描画

xx <- seq(min(x1), max(x1), length = 200)

yy <- coef(fit1)[1] + coef(fit1)[2] * xx

overdraw(points(xx, yy, type = "l", col = "red"))

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

では、逆変換によって、元のデータに戻してみよう。併せて、通常の回帰分析（OLS）、そして、前回行った重み付き最小二乗法（WLS）の結果も示す。

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

yy_inv <- (powerTransform(d$y)$lambda*yy + 1)^(1/powerTransform(d$y)$lambda)

plotn(y ~ x1)#図12の描画

overdraw(points(xx, yy_inv, type = "l", col = "red"))

fit2 <- lm(y ~ x1)#通常の回帰分析

yy <- coef(fit2)[1] + coef(fit2)[2] * xx

overdraw(points(xx, yy_inv, type = "l", col = "red"),

         points(xx, yy, type = "l", col = "green"))

w <- abs(y - (coef(fit2)[1] + coef(fit2)[2] * x1))

fit3 <- lm(y ~ x1, weights = (1/w)^2)#重み付き最小二乗法

yy <- coef(fit3)[1] + coef(fit3)[2] * xx

overdraw(points(xx, yy, type = "l", col = "blue"))

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

今回の結果では、OLSもWLSも結果に大きな違いはない。というのも、回帰直線の上下でデータが均等にばらついている場合は、OLSもWLSも大差のない結果が出力される。前回紹介したように、WLSは外れ値による分散のアンバランスに効果を発揮する。Box-Cox変換の逆変換による回帰線は直線となる。何を重要視しているかが異なるため、Box-Cox変換とWLSは単純比較はできない。いずれも、ある意味で正しい分析である。

　もうひとつ、以下の例の解析をしてみよう。

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

x1 <- c(15.4, 22.45, 2.91, 25.76, 25.14, 11.78, 18.28, 14.57, 18.17, 27.16, 

        3.44, 16.31, 1.63, 22.56, 7.66, 27.55, 12.47, 6.55, 4.72, 4.46, 

        20.11, 19.96, 16.61, 19.25, 16.87, 18.56, 18.12, 5.76, 10.02, 3.43,

        13.79, 11.73, 4.08, 24.7, 19.69, 22.05, 25.81, 6.22, 5.89, 21.69, 

        17.97, 26.74, 13.77, 19.3, 26.92, 20.21, 18.14, 25.61, 20.57, 9.2, 

        1.65, 24.35, 20.63, 25.47, 2.19, 11.5, 17.91, 23.29, 29.99, 5.53, 

        18.19, 10, 14.87, 23.38, 16.25, 5.79, 26.43, 15.24, 14.62, 10.88, 

        2.2, 19.73, 10.33, 2.53, 9.33, 1.19, 13.31, 3.45, 12.32, 22.16, 

        25.67, 14.67, 15.34, 19.03, 10.53, 23.21, 9.84, 28.64, 13.75, 0.19, 

        8.37, 0.56, 7.89, 16.23, 5.31, 1.72, 6.89, 19.32, 5.85, 3.19)

y <- c(5.98, 6.78, 1.79, 18.21, 8.94, 1.12, 3.51, 5.77, 6.23, 9.25, 1.09, 

       8.16, 0.95, 4.96, 2.78, 8.95, 4.65, 0.91, 0.58, 1.63, 4.79, 4.6, 

       5.65, 7.81, 8.99, 10.15, 5.61, 1.38, 4.91, 2.1, 1.74, 2.32, 1.05, 

       7.91, 11.39, 14.77, 4.6, 2.55, 2.53, 6.85, 6.31, 10.39, 2.1, 7.2, 

       9.76, 5.24, 7.97, 11.59, 6.58, 2.87, 1.27, 14.66, 8.32, 14.03, 2.13, 

       2.79, 6.06, 8.03, 24.4, 0.71, 7.74, 1.8, 5.43, 13.17, 1.87, 0.59, 

       6.88, 4.82, 3.43, 1.62, 2.18, 8.39, 2.57, 1.86, 2.37, 0.56, 3.58, 

       3.88, 4.08, 5.8, 13.6, 2.34, 1.93, 6.31, 1.49, 11.39, 1.77, 26.22, 

       5.22, 0.96, 1.48, 1.56, 2.44, 2.61, 1.59, 1.1, 4.14, 5.71, 2.06, 2.2)

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

head(d)

##      x1     y

## 1 15.40  5.98

## 2 22.45  6.78

## 3  2.91  1.79

## 4 25.76 18.21

## 5 25.14  8.94

## 6 11.78  1.12

plotn(y ~ x1, data = d)#図13の描画

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

今度のデータは下に凸の曲線のようなデータであり、この場合はそもそもが直線に当てはめるのに無理がある。

　このデータもBox-Cox変換して、回帰分析してみよう。

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

powerTransform(d$y)

## Estimated transformation parameter 

##        d$y 

## 0.07636666

y_t <- bcPower(d$y, powerTransform(d$y)$lambda)

plotn(y_t ~ x1)#図14の描画

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

すると、歪みが減り、線形回帰ができそうなデータの分布となった。データの分布を確認してみよう。

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

histn(y)#図15の描画

histn(y_t)#図16の描画

m <- mean(y)

s <- sd(y)

ry <- seq(min(y)-5, max(y), length = 200)

cy <- cdf(p = y, v = ry)

plotn(ry, cy, type = "l", xlab = "y", ylab = "CDF")#図17の描画

overdraw(points(ry, pnorm(ry, mean = m, sd = s), type = "l", col = "red"))

m <- mean(y_t)

s <- sd(y_t)

ry <- seq(min(y_t)-2, max(y_t)+1, length = 200)

cy <- cdf(p = y_t, v = ry)

plotn(ry, cy, type = "l", xlab = "y", ylab = "CDF")#図18の描画

overdraw(points(ry, pnorm(ry, mean = m, sd = s), type = "l", col = "red"))

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

明らかに、変換後は正規分布に近づいている。では、回帰分析を行う。

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

fit1 <- lm(y_t ~ x1)

plotn(y_t ~ x1)#図19の描画

xx <- seq(min(x1), max(x1), length = 200)

yy <- coef(fit1)[1] + coef(fit1)[2] * xx

overdraw(points(xx, yy, type = "l", col = "red"))

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

では、逆変換によって、元のデータに戻してみよう。併せて、ガンマ分布、対数リンクを仮定した一般化線形モデル（GLM）の結果も示す（実はこのデータはガンマ分布に従うように生成した）。

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

yy_inv <- (powerTransform(d$y)$lambda*yy + 1)^(1/powerTransform(d$y)$lambda)

plotn(y ~ x1)#図20の描画

overdraw(points(xx, yy_inv, type = "l", col = "red"))

fit2 <- glm(y ~ x1, family = Gamma(link = "log"))

plotn(y ~ x1)

yy <- exp(coef(fit2)[1] + coef(fit2)[2] * xx)

overdraw(points(xx, yy_inv, type = "l", col = "red"),

         points(xx, yy, type = "l", col = "blue"))

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

回帰結果は類似しているが、微妙に異なっている。これも、どちらが正しいというものではない。しかし、覚えておくとよいのは、GLMの青線は変換前のデータの平均値を通っている。一方で、Box-Cox変換の回帰の赤線は、変換後のデータの平均値を通っている。確率分布の変数変換をf(x)と表せば、変換した後の確率分布の期待値E(f(x))は、変換前の確率分布の期待値を変数変換したものf(E(x))とは一致しない。それの最たる例が不偏分散の平方根が不偏標準偏差にならない問題である。 油断していると勘違いするので気を付けよう。

補足：Box-Cox変換の連続性

　Box-Coxはλが0か否かで、関数が異なっている。ちゃんとこれらが連続的であることを以下に示す。

補足：Yeo-Johnson変換

　Box-Cox変換はyi < 0のときに変換ができない。そこで、実数全体で変換できるようにした変換がYeo-Johnson変換である。その変換は以下の通り。

結構ややこしいうえに、ネットの日本語文献ではなぜか間違った式になっているので注意が必要だ。また、上記の表記は結構混乱を産むので、以下のようにまとめ直した。

あるいは、以下のようにλ≠0, 2のときに一つの式にまとめることも可能である。sign関数は、入力した値が正なら+1、負なら-1を返す。

変換の様子は以下の通り。

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

YeoJohn <- function(x, lambda){

  mapply(function(tmp1, tmp2){

    if(!tmp2 == 0 & !tmp2 == 2){

      ((abs(tmp1) + 1)^(sign(tmp1)*(tmp2-1)+1) - 1)/(tmp2-1+sign(tmp1))

    } else {

      if(tmp2 == 0) {

        if(tmp1 >= 0) {

          log(tmp1 + 1)

        } else {

          -((-tmp1+1)^(2-tmp2) - 1)/(2-tmp2)

        }

      } else {

        if(tmp1 >= 0) {

          ((tmp1+1)^tmp2 - 1)/tmp2

        } else {

          -log(-tmp1 + 1)

        }

      }

    }

  }, x, lambda)

}

x <- seq(-2, 2, length = 200)

d <- data.frame(x = x, 

                y1 = YeoJohn(x, lambda = -2),

                y2 = YeoJohn(x, lambda = -1),

                y3 = YeoJohn(x, lambda = 0),

                y4 = YeoJohn(x, lambda = 0.5),

                y5 = YeoJohn(x, lambda = 1),

                y6 = YeoJohn(x, lambda = 2),

                y7 = YeoJohn(x, lambda = 3))

plotn(d[,1], d[,2:8], mode = "m", type = "l", ylim = c(-2,2),

      xlab = "y", ylab = "transformed y",

      legend = T, leg.title = "lambda", leg.lab = c(-2,-1,0,0.5,1,2,3),

      pch.leg = NA, line = T, leg.sp = 4)#図21の描画

overdraw(abline(v = 0, lty = 2),

         abline(h = 0, lty = 2))

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