学部向けのレクチャーノートを見ていると言及されていることが多い一方、使いどころがあるようでないのが加重最小二乗法（WLS）です。連立方程式モデルやパネルデータ分析のランダム効果モデルで拡張的な手法は使うわけですが、WLS自体は使い道が無いと言うか。
不均一分散自体はよく見られる一方、それに配慮したからといって推定結果を大きく左右するかと言うとそうでもなく、逆に補正を誤ると奇妙な推定結果が出てきます。都道府県ごとの希少な病気の発生件数が典型例としてあげられる*1のですが、必要になることは意外に乏しいです。
よく知られたヤマメのデータセットでOLSをかけた後、不均一分散が入っているかBP検定をかけてみましょう。
model <- Sepal.Length ~ Species r_lm <- lm(model, iris) library(lmtest) bptest(r_lm)
studentized Breusch-Pagan test data: r_lm BP = 12.34, df = 2, p-value = 0.002091
帰無仮説の均一分散が1%有意で棄却されました。品種が違えば分散も異なるのは直観的でもあります。不均一分散があると看做しましょう。FGLSになっていますが、加重最小二乗法を用います。
r_res <- lm(log(abs(residuals(r_lm))) ~ Species, iris) # 非負にするために分散の大きさが指数関数に従うと仮定している w <- 1/exp(predict(r_res)) # lm関数のウェイトは平方根をとらなくてよい r_wls <- lm(model, iris, weights = w)
OLSとWLSの推定結果を見比べてみましょう。
coef(summary(r_lm))
Estimate Std. Error t value Pr(>|t|) (Intercept) 5.006 0.07280222019 68.761639227 1.134286186e-113 Speciesversicolor 0.930 0.10295788717 9.032819394 8.770194241e-16 Speciesvirginica 1.582 0.10295788717 15.365505679 2.214821349e-32
coef(summary(r_wls))
Estimate Std. Error t value Pr(>|t|) (Intercept) 5.006 0.05151993771 97.16626654 2.888838414e-135 Speciesversicolor 0.930 0.09250143061 10.05389856 2.011782048e-18 Speciesvirginica 1.582 0.09808115552 16.12950002 2.472870822e-34
推定された係数は一致しています。標準誤差がかなり小さくなっているので、p-hacking的に有用なわけですが、分散の特定化の誤りの可能性も入るのでWLSを信じるべきとも言い難いです。
もちろんシミュレーションでは有用な結果を出せます。
# データセット生成 set.seed(1234) n <- 50 x <- runif(n, 0, 1) z <- c(rep(1, n - 3), rep(10, 3)) y <- 1 + 2*x + rnorm(n, sd = z) # 推定 r_lm <- lm(y ~ x) r_res <- lm(log(abs(residuals(r_lm))) ~ x + z) w <- 1/exp(predict(r_res)) r_wls <- lm(y ~ x, weights = w) # 推定結果の確認 coef(summary(r_lm))
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.528502 0.9675744 0.5462133 0.5874493 x 3.423524 1.8009453 1.9009594 0.0633206
coef(summary(r_wls))
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.5718522 0.3385899 1.688923 0.097721188 x 2.4158260 0.7706556 3.134767 0.002932745
DGPのxの係数は2で、WLSの推定結果の方が近いです。標準誤差もかなり小さいですね。
ただし、DGPの分散を小さくするとOLSと結果が変わらなくなる一方、大きくするとOLSもWLSも大きく外した推定結果になったりするので、実際に使うかは悩みそうです。
教科書に沿った推定
lm関数にweightsをつけるだけでは教科書の説明との対応が悪いので、変数にウェイトをかける方法での推定も示しておきます。
versicolor <- iris$Species == "versicolor" virginica <- iris$Species == "virginica" r_lm <- lm(Sepal.Length ~ versicolor + virginica, iris) r_res <- lm(log(abs(residuals(r_lm))) ~ versicolor + virginica, iris) w <- 1/exp(predict(r_res)) s_w <- sqrt(w) r_wls_ln <- lm(I(s_w*Sepal.Length) ~ 0 + I(s_w) + I(s_w*versicolor) + I(s_w*virginica), iris) coef(summary(r_wls_ln))
Estimate Std. Error t value Pr(>|t|) I(sw) 5.006 0.05151993771 97.16626654 2.888838414e-135 I(sw * versicolor) 0.930 0.09250143061 10.05389856 2.011782048e-18 I(sw * virginica) 1.582 0.09808115552 16.12950002 2.472870822e-34
推定量がlmにweightsをつけた場合と一致しますね。
なお、ウェイトは平方根をとっています。また、ウェイトを変数にかける都合で、ダミー変数を展開しています。
*1:人口が多いところは平均回帰するので分散が小さい一方、人口が少ないところはそうでないので分散が大きくなる。
