ほくそ笑む

R言語と統計解析について

log 変換する?しない?AICでモデル比較するときの注意点

データを分析にかける前に、出力変数を log 変換する、というのはよくあることだと思います。
次のデータを見て下さい。

このデータ、線形モデルに当てはめる前に log 変換したほうがよさそうだなーというのが見てとれます。
それもそのはず、このデータは次のように作っています。

N <- 100
x <- runif(N, min = 1, max = 2)
y <- exp(x + rnorm(N, sd = 0.3))
data <- data.frame(x, y)

それでは、log 変換しないバージョンと、するバージョンでモデルを作成して、AIC を比較してみましょう。

model <- lm(y ~ x, data)
model.log <- lm(log(y) ~ x, data)
aic <- AIC(model, model.log)
print(aic)
##           df    AIC
## model      3 364.09
## model.log  3  38.94

このように、log 変換するバージョンの方が AIC が小さいので、log 変換したほうが良いモデルだということがわかります。

こんどは、次のデータを見てみましょう。

このデータは、次のように生成されているので、log 変換はしないほうが良いはずです。

N <- 100
x2 <- runif(N, min = 1, max = 2)
y2 <- x2 + rnorm(N, sd = 0.4)
data2 <- data.frame(x2, y2)

それでは、log 変換しないバージョンと、するバージョンでモデルを作成して、AIC を比較してみましょう。

model2 <- lm(y2 ~ x2, data2)
model.log2 <- lm(log(y2) ~ x2, data2)
aic2 <- AIC(model2, model.log2)
print(aic2)
##            df   AIC
## model2      3 96.48
## model.log2  3 37.99

なんと、log 変換するバージョンのモデルの方が、AIC が小さくなってしまいました。
これはなぜでしょう?

変数変換のワナ

実は、変数変換を行う前後のモデルは AIC でそのまま比較できません。
なぜなら、AIC の骨格を成す尤度が、変数変換前後で尺度が変わってしまうためです。
ここでは、変数変換後の AIC を変数変換前の尺度にするための補正項を求めてみます。

今、出力変数を変換する関数を  w(y) とします。
このとき、変換前の確率密度関数  f(y) と変換後の確率密度関数  g(w) の間には、次のような関係が成り立ちます。
  f(y) = g(w(y))\frac{dw}{dy}
したがって、変換前の尺度での尤度  Lik_y^* は、変換後の尺度での尤度  Lik_w を用いて次のように表すことができます。
  Lik_y^* = \prod_{i=1}^N f(y_i) = \prod_{i=1}^N g(w(y_i))\frac{dw}{dy}(y_i) = Lik_w \prod_{i=1}^N \frac{dw}{dy}(y_i)
よって、対数尤度は次のようになります。
  logLik_y^* = logLik_w + \sum_{i=1}^N \log\big(\frac{dw}{dy}(y_i)\big)
これにより、変換前の尺度での AIC を求めることができます。
  AIC_y^* = -2logLik_y^* + 2p = -2logLik_w - 2\sum_{i=1}^N \log\big(\frac{dw}{dy}(y_i)\big) + 2p
  AIC_y^* = AIC_w - 2\sum_{i=1}^N \log\big(\frac{dw}{dy}(y_i)\big)
今、 w(y) = \log(y) とすると、 w^\prime = \frac{1}{y} なので、対数変換後のモデルの変換前の尺度での AIC
  AIC_y^* = AIC_w -2\sum_{i=1}^N \log\big(\frac{1}{y_i}\big)
となります。
すなわち、対数変換前後のモデルを AIC で比較したい場合は、対数変換後の AIC に補正項を加えてやればいいだけです。

Rでやってみる

上記の例を用いて、R で実際にやってみましょう。

まずは、最初の例です。

adj <- -2 * sum(log(1/y))
cbind(aic, adj.AIC = aic$AIC + c(0, adj))
##           df    AIC adj.AIC
## model      3 364.09   364.1
## model.log  3  38.94   334.8

やはり、補正項を加えても、log 変換後の方が良いモデルとなっています。いいですね。

それでは問題となる2番目の例です。

adj2 <- -2 * sum(log(1/y2))
cbind(aic2, adj.AIC = aic2$AIC + c(0, adj2))
##            df   AIC adj.AIC
## model2      3 96.48   96.48
## model.log2  3 37.99  107.58

こんどは log 変換しないモデルの方が良いと出ました。やったね!

これにより、log 変換前後のモデルを AIC で比較することができるようになりました。

もうすこし一般的に

もう少し一般的に書くと、次のように書けると思います。

library(numDeriv)
trans.func <- sqrt
adj <- -2*sum(log(grad(trans.func, x=y)))

このように補正項を求めれば、log 変換だけでなく、平方根変換とかでも使えるようになります。

まとめ

この記事は R Advent Calendar 2013 の 23 日目の記事です。
こういう情報 Web 上に無いな〜と思って書いてみました。
みなさん変数変換するときはどうしてるんでしょうかね?
出力変数が2変数以上のときはもうちょっと複雑でヤコビアンを使うことになります。
ヤコビアンも上記の numDeriv パッケージにあるので R で簡単に実装できると思います。

以上