東京大学・株式会社Nospareの菅澤です．
前回の外れ値とロバスト推定(3)では，ダイバージェンスによる一般的なロバスト推定と，それを線形回帰モデルに適用した場合の例について紹介しました．
今回は，簡単な実験を通して以下の2点について考察していきます．
- 外れ値が混入しているときに通常の線形回帰を当てはめた場合の問題点
- これまでに紹介した線形回帰モデルのロバスト推定法の比較
擬似データの生成
以下の回帰モデルからデータを生成します．
y_i=\beta_0+\sum_{k=1}^p\beta_kx_{ik} + \varepsilon_i, \ \ \ i=1,\ldots,n.
ここで，$n=100$（サンプルサイズ）, $p=10$（説明変数の数）とし，真の回帰係数を$\beta_0=0.5$, $\beta_1=\beta_4=0.3$, $\beta_7=\beta_{10}=1$（残りの係数は0）とします．
また，誤差項$\varepsilon_i$の分布として，以下のような汚染モデルを考えます．
\varepsilon_i \sim (1-\omega) N(0, \sigma^2) + \omega N(10, \sigma^2).
第1成分が通常の誤差分布を表し，第2成分が外れ値を生成する要素に相当します．また，$\omega$が汚染率と呼ばれる量で，外れ値の割合に相当します．今回は$\omega=0.05$ (データのうち$5\%$が外れ値)と設定します．
説明変数ベクトル$(x_{i1},\ldots,x_{ip})$は各ペアごとの相関が$0.2$（すなわち Cor$(x_{ij},x_{ik})=0.2, j\neq k$）となる多変量正規分布から生成します．
上記を実行するコードは以下の通りです．
set.seed(1)
n <- 100 # sample size
p <- 10 # number of covariates
rho <- 0.2
mat <- matrix(NA, p, p)
for(k in 1:p){
for(l in 1:p){ mat[k,l] <- rho^(abs(k-l)) }
}
X <- mvrnorm(n, rep(0, p), mat) # design matrix
sigma <- 0.5 # error standard deviation
Beta <- rep(0, p)
Beta[c(1, 4, 7, 10)] <- c(0.3, 0.3, 1, 1) # non-zero coefficients
int <- 0.5 # intercept
om <- 0.05 # contamination ratio
aa <- 10 # outlier location
ch <- rbinom(n, 1, om)
ep <- (1-ch)*rnorm(n, 0, sigma) + ch*rnorm(n, aa, sigma)
y <- int + as.vector(X%*%Beta) + ep
(非ロバストな)線形回帰モデルの当てはめ
まず通常の線形回帰モデルを当てはめてみます．
fit.LM <- lm(y~X)
summary(fit.LM)
結果は以下のようになります．
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.98539 0.22948 4.294 4.46e-05 ***
X1 -0.25364 0.21177 -1.198 0.234211
X2 0.22963 0.22167 1.036 0.303047
X3 -0.56160 0.26664 -2.106 0.038005 *
X4 0.79686 0.26001 3.065 0.002884 **
X5 0.08327 0.23260 0.358 0.721183
X6 -0.08134 0.24550 -0.331 0.741165
X7 0.78998 0.21659 3.647 0.000446 ***
X8 0.13426 0.21559 0.623 0.535052
X9 -0.07741 0.22223 -0.348 0.728392
X10 0.89313 0.24069 3.711 0.000359 ***
X1, X4, X7, X10に対する係数の真値がそれぞれ0.3, 0.3, 1, 1であったことを思い出すと
- X1は有意な影響を与えているはずなのに，その効果を拾えることができていない
- 逆に，X3は影響を与えないはずの変数であるが，(若干)有意な負の効果を与える結果になっている
これは限定的な実験結果のため統一的なことは言えませんが，外れ値を考慮しない解析によって，本来推定したい構造とは異なった推定結果が得られてしまう可能性があることがわかります．
ロバストな推定手法
density power divergenceによる推定
前回の記事で紹介したように，density power divergenceを用いた($\beta$, $\sigma^2$)の目的関数は以下のような形でした．
D_{\alpha}(\beta,\sigma^2)=\frac{1}{\alpha}\sum_{i=1}^n \phi(y_i;x_i^\top\beta,\sigma^2)^{\alpha} - n(2\pi\sigma^2)^{-\alpha/2}(1+\alpha)^{-3/2}.
(ここで$x_i$は切片項も含んだ表現になっていることに注意してください．)
上記の関数の最小化は以下の手順を収束するまで繰り返すアルゴリズムで推定することができます．
- $i=1,\ldots,n$に対して重み$w_i^{(t)}=\phi(y_i;x_i^\top\beta^{(t)}, \sigma^{(t)2})$を計算する．
- $\beta$の更新
\beta^{(t+1)} \leftarrow \left(\sum_{i=1}^n w_i^{(t)} x_ix_i^\top\right)^{-1}\sum_{i=1}^n w_i^{(t)} x_iy_i
- $\sigma^2$の更新
\sigma^{(t+1)2} \leftarrow \left(\sum_{i=1}^n w_i^{(t)} -\frac{n\alpha}{(1+\alpha)^{3/2}}\right)^{-1}\sum_{i=1}^n w_i^{(t)} (y_i-x_i^\top \beta^{(t+1)})^2
また，$\beta$の推定量$\hat{\beta}_{\alpha}$の漸近分散共分散行列は以下の形で与えられます(Ghosh and Basu (2013)を参照)．
\text{Var}(\hat{\beta}_{\alpha})\approx \frac{C(2\alpha)}{C(\alpha)^2}\hat{\sigma}^2\left(\sum_{i=1}^n x_ix_i^\top\right)^{-1}, \ \ \
C(\alpha)=(2\pi)^{-\alpha/2}(1+\alpha)^{-3/2}.
この行列の対角成分の平方根を計算することで$\hat{\beta}_{\alpha}$の標準誤差を計算することができます．
上記の手順を実装したものが以下のDPD.LM関数です．
# function for fitting linear models based on density power divergence
DPD.LM <- function(y, X, alpha=0.2){
XX <- cbind(1, X)
init.fit <- rlm(y~X, method="MM")
hBeta <- coef(init.fit)
hsig <- summary(init.fit)$sigma
error <- 10^(-8)
maxitr <- 100
for(j in 1:maxitr){
hBeta.old <- hBeta
ww <- dnorm(y, as.vector(XX%*%hBeta), hsig)^alpha
hBeta <- as.vector( solve(t(XX)%*%(XX*ww))%*%t(XX)%*%(y*ww) )
a1 <- sum( ww*(y-as.vector(XX%*%hBeta))^2 )
a2 <- sum(ww - alpha/(1+alpha)^(3/2))
hsig <- sqrt( a1/a2 )
dd <- sum(abs(hBeta-hBeta.old))/sum(abs(hBeta.old))
if(dd<error){ break }
}
c1 <- (2*pi)^(-alpha/2)*hsig^(-alpha-2)*(1+alpha)^(-3/2)
c2 <- (2*pi)^(-alpha)*hsig^(-2*alpha-2)*(1+2*alpha)^(-3/2)
SE <- sqrt( c2/c1^2 * diag( solve(t(XX)%*%XX) ) )
res <- list(beta=rbind(hBeta, SE), sigma=hsig)
return(res)
}
数値的な比較
まずはHuberとTukeyの方法を適用します．
library(MASS)
fit.Huber <- rlm(y~X, psi=psi.huber)
summary(fit.Huber)
fit.Tukey <- rlm(y~X, method="MM")
summary(fit.Tukey)
Huber法の結果
Coefficients:
Value Std. Error t value
(Intercept) 0.5171 0.0561 9.2175
X1 0.1796 0.0518 3.4697
X2 0.0437 0.0542 0.8057
X3 -0.0121 0.0652 -0.1849
X4 0.2828 0.0636 4.4490
X5 -0.0859 0.0569 -1.5111
X6 0.0075 0.0600 0.1247
X7 0.9474 0.0529 17.8946
X8 0.0330 0.0527 0.6255
X9 0.0744 0.0543 1.3692
X10 0.8663 0.0588 14.7242
Residual standard error: 0.4615
Tukey法の結果
Coefficients:
Value Std. Error t value
(Intercept) 0.4577 0.0537 8.5252
X1 0.2438 0.0495 4.9205
X2 0.0132 0.0519 0.2541
X3 0.0675 0.0624 1.0822
X4 0.2358 0.0608 3.8772
X5 -0.1023 0.0544 -1.8801
X6 0.0054 0.0574 0.0936
X7 0.9759 0.0507 19.2618
X8 0.0168 0.0504 0.3323
X9 0.0839 0.0520 1.6134
X10 0.8761 0.0563 15.5607
Residual standard error: 0.4515
どちらの結果も真のデータ生成過程で影響を与えていた変数(X1, X4, X7, X10)は有意な結果になっており，それ以外の変数にはあまり効果がないという結果が得られていることがわかります．
また，$\sigma$の推定値(Residual standard errorの値)も真の誤差分散$\sigma=0.5$に近い値になっています．
次にdensity power divergenceによる方法を適用します．
fit.DPD <- DPD.LM(y, X, alpha=0.2)
DPD.est <- t(fit.DPD$beta)
dimnames(DPD.est)[[1]] <- c("(Intercept)", paste0("X", 1:p))
dimnames(DPD.est)[[2]] <- c("Value", "Std. Error")
DPD.est
fit.DPD$sigma
結果は以下のようになります．
> DPD.est
Value Std. Error
(Intercept) 0.453925340 0.04903743
X1 0.248241901 0.04525350
X2 0.012178367 0.04736875
X3 0.071466057 0.05697883
X4 0.237668410 0.05556123
X5 -0.101432364 0.04970384
X6 0.005358339 0.05245997
X7 0.980595238 0.04628256
X8 0.015517614 0.04606964
X9 0.079880689 0.04748713
X10 0.880408565 0.05143191
> fit.DPD$sigma
[1] 0.4659322
Huber法やTukey法と似たような結果が得られます．
ただし，Huber法やTukey法と比べてDPDの場合は$\alpha$の値を設定する必要があります．
詳細な結果はここでは出しませんが，今回のケースでは$\alpha$を$0.1$から$0.5$まで動かしてもほとんど結果は変わりませんでした．
もう少し定量的に比較をするため，それぞれの方法について二乗誤差を計算してみます．
> sqrt( mean( (coef(fit.LM) - c(int, Beta))^2 ) )
[1] 0.3371287
> sqrt( mean( (coef(fit.Huber) - c(int, Beta))^2 ) )
[1] 0.0686443
> sqrt( mean( (coef(fit.Tukey) - c(int, Beta))^2 ) )
[1] 0.06573471
> sqrt( mean( (fit.DPD$beta[1,] - c(int, Beta))^2 ) )
[1] 0.06435719
この結果からHuber<Tukey<DPDの順で精度の高い推定値を与えていることがわかります．（TukeyとDPDの差はかなり小さいですが．）
また，通常の最小二乗法による推定結果の精度はかなり低いこともわかります．
汚染率を上昇させてみる
最後に，汚染率$\omega$を上昇させたときの結果について簡単に紹介します．
汚染率以外の設定はそのままで，$\omega=0.1, 0.15, 0.2$の場合にデータを生成し，3つのロバストな手法を適用します．
回帰係数$\beta$の二乗誤差と$\sigma$の推定値は以下のようになりました．
| 誤差 ($\omega=0.1$) | $\hat{\sigma}$($\omega=0.1$) | 誤差($\omega=0.15$) | $\hat{\sigma}$($\omega=0.5$) | 誤差($\omega=0.2$) | $\hat{\sigma}$($\omega=0.2$) | |
|---|---|---|---|---|---|---|
| Huber | 0.066 | 0.502 | 0.079 | 0.621 | 0.102 | 0.701 |
| Tukey | 0.071 | 0.472 | 0.072 | 0.540 | 0.068 | 0.645 |
| DPD | 0.068 | 0.468 | 0.070 | 0.477 | 0.067 | 0.488 |
この結果から，汚染率が高くなるとHuber法はあまり機能しない可能性が示唆されますが，回帰係数の推定精度の点ではTukey法とDPDはほぼ同等のパフォーマンスになっています．
一方で，$\sigma$の推定に関しては，他の方法に比べてDPDによる推定精度が高いことが確認できます．これは前回の記事で紹介したように，DPDが「$\sigma$を含めた統計モデルとしての線形回帰モデルのロバスト推定」を考えていることに起因していると考えられます．
おわりに
今回は簡単な実験を通して，線形回帰モデルにおける外れ値の影響といくつかのロバスト推定法の違いについて考察しました．
今回の実験結果では，回帰係数の推定精度に関してはHuber法<Tukey法<DPD法（Tukey法はDPD法ほぼ同等）といった関係性が確認できました．
回帰係数のみに興味がある場合はTukey法でも十分そうですが，この方法は線形回帰モデル以外にそのまま適用することはできません．一方で，DPDは統計モデルをロバストに推定するための一般的な概念のため，一般化線形モデルにも適用することができます．(Ghosh and Basu (2016)などを参照．)
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
