optimってあんまり信用できないなぁ、って話

タイトルの通りです。

仕事で使うことがあったのでRのoptimを使って回帰を解いてみたのですが、これが意外に安定しません。変数の数なのか、ダミー変数が含まれるとダメなのか、原因についてはよくわかりませんが想定以上に解がバラついてしまいました。

実行したコードは以下の通りです。まずは以下のようにデータとoptimに渡す目的関数を定義します：

## サンプルデータとしてirisを使用
dat <- iris
datMdl <- as.data.frame(model.matrix(~.-1, data = dat))

## Speciesを含まないデータ
dat01 <- datMdl[, 1:4]

## Speceisを含むデータ（virginicaを除外）
dat02 <- datMdl[, -7]

## 最小二乗法を以下のように定義
myLM <- function(b, dat) {
  Y   <- as.matrix(dat[,  1])
  X   <- as.matrix(cbind(1, dat[, -1]))
  b   <- as.matrix(b)
  prd <- X %*% b
  e   <- Y - prd
  return(t(e) %*% e)
}

以上の条件で、まずはSpeciesを含まないデータを用いてlmとoptimをそれぞれ実行します。

## Speciesを含まないデータでSepal.Lengthをyとして回帰をかける
resLM_01 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
               data=dat01)

## optimによる回帰係数の取得
init_coef <- rep(1, 4)
resOpt_01 <- optim(par=init_coef, fn=myLM, dat=dat01[, 1:4])

## 係数の確認。だいたい同じ
cbind(
  summary(resLM_01)$coef[, 1],
  resOpt_01$par)

                   [,1]       [,2]
(Intercept)   1.8559975  1.8589403
Sepal.Width   0.6508372  0.6501880
Petal.Length  0.7091320  0.7087105
Petal.Width  -0.5564827 -0.5559860

コメントでも書いてありますが、係数は概ね一致します。続けてSpeciesを含めたデータによりlmおよびoptimをそれぞれ再度実行します。

## 次にSpeciesを含むデータで同様に回帰を実行
resLM_02 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + 
                 Speciessetosa + Speciesversicolor,
               data=dat02)

## 同じく、optimによる回帰係数の取得
init_coef <- rep(1, 6)
resOpt_02 <- optim(par=init_coef, fn=myLM, dat=dat02)

## 係数の確認。結構違う
cbind(
  summary(resLM_02)$coef[, 1],
  resOpt_02$par)

                        [,1]       [,2]
(Intercept)        1.1477685  0.6123026
Sepal.Width        0.4958889  0.4621860
Petal.Length       0.8292439  0.8993939
Petal.Width       -0.3151552 -0.1967130
Speciessetosa      1.0234978  1.5512063
Speciesversicolor  0.2999359  0.4695160

Rの最適化ソルバーはoptim以外にもありますし(optimizeとか)、optimを使うにしてもパラメータや用いるアルゴリズム(Nelder-Mead、BFGSとか)など設定できる箇所は色々とあります。なのでこの結果だけを元にoptimを信用できないとは言えないのですが、このくらいのデータ・モデルでこれほどlmとの差が出てしまうと複雑な問題への適用はちょっとためらってしまいますね。

なお当初の目的としてはMarketing Mix Modelによる広告効果の推定を行う際に、残存効果の減衰率をoptimで推定してみようというものでした。optimの結果があまり芳しくないためこういった検証を行っていたのですが、この辺はいずれ記事にまとめようと思います。それからoptimはオプションでHessian(ヘッセ行列)を出力することにより標準誤差を計算できますので、その辺もいずれ書いてみようと思います。