Technically Impossible

Lets look at the weak link in your statement. Anything "Technically Impossible" basically means we haven't figured out how yet.

GVLMAの紹介-Rでの残差を用いた線形モデルの包括的検査

回帰分析から線形モデルを生成したならば、満たすべき仮定を検査する必要がある。その仮定とは、

  • 独立性
  • 分散の均一性
  • 誤差の正規性
  • モデルの線形性

Rを用いるとき、回帰式より生成したlmオブジェクトをプロットすることで、残差と外れ値を視覚的に検証し、試行錯誤を繰り返すことになる。あるいは、次のページで説明されているような手続きを実施することもある。
ademos.people.uic.edu

GVLMA (Global Validation of Linear Model Assumptions)は、線形モデルの残差に注目し、標準化した残差の正規分布を検証することで、先の仮定を包括的に検証しようとする。具体的には、次の事柄を検証している。

  1. 残差の分布の歪度
  2. 残差の分布の尖度
  3. 目的変数の分布(不適切なリンク関数の使用)
  4. 残差の等分散性

GVLMAの出力結果を参照することで、線形モデルの検査が容易になる。

例によって、参照先となるリンクは投稿末尾にまとめている。

GVLMAを用いない場合

線形モデルを生成し、検証のための情報を出力するまでのコードは次のようになる。RStudioに含まれているデータセット、stacklossを用いて線形モデルを出力し、そのプロットと一個抜き交差検証の結果を出力している。
リサーチャーは各出力結果を考察し、線形モデルの試行錯誤を繰り返すことになる。

R code

mylm = lm(stack.loss ~ stack.x)
summary(mylm)
plot(mylm, which = c(1:6))
influence.measures(mylm)


f:id:espio999:20191022214248p:plain
f:id:espio999:20191022214302p:plain

GVLMAを用いる場合

検証結果

同じ線形モデルを用いて、GVLMAを実行する。線形モデルからGVLMAオブジェクトを生成しプロットする。さらに同オブジェクトから一個抜き交差検証を実施し、その結果をプロットしている。

R code

library(gvlma)

mygvlma = gvlma(mylm)
summary(mygvlma)
par(mfrow = c(3, 3))
plot(mygvlma, onepage = FALSE)
mydel = deletion.gvlma(mygvlma)
summary(mydel)
plot(mydel, onepage = FALSE)


f:id:espio999:20191022214356p:plain
まずGVLMAの検証結果が表示される。それぞれ次の事柄を指している。

Global Stat        2.75775  0.5991 Assumptions acceptable. 包括的検証の結果
Skewness           0.13015  0.7183 Assumptions acceptable. 残差分布の歪度
Kurtosis           0.01009  0.9200 Assumptions acceptable. 残差分布の尖度
Link Function      2.59888  0.1069 Assumptions acceptable. 目的変数の分布
Heteroscedasticity 0.01863  0.8914 Assumptions acceptable. 残差分布の等分散性

GVLMAの論文を引用すると、SkewnessからHeteroscedasticityまで、次の事柄を検定している。S^21からS^24が、それぞれSkewnessからHeteroscedasticityに該当している。これら検定の統計値を合計し、包括的な検定を実施したものがGlobal Statだ。

(i) Skewed error distributions will usually be indicated by large values of the statistic Sˆ21
(ii) Deviations from the normal distribution kurtosis of the true error distribution will be generally revealed by large values of statistic Sˆ22
(iii) The use of a misspecified link function, possibly due to the absence of other predictor variables in the model, will mostly be detected by large values of the statistic Sˆ23
(iv) The presence of heteroscedastic errors and/or dependent errors will typically manifest in large values of the statistic Sˆ24

もし、どこかの項目に「Assumptions NOT satisfied!」と表示されたならば、残差が正規分布していないことを示唆している。
Link Functionは、説明変数の取捨選択を誤り、結果として不適切なリンク関数を用いているのではないか、と言うことを示唆している。この項目が不適切な場合、説明変数を選び直せ、と言うことだ。

検証結果のプロット

GVLMAは次の項目のプロットを出力する。今回のデータは時系列ではないので、時系列データのプロットは無視する。これは後述する外れ値に関するプロットでも同様だ。

(a) the response versus each of the predictors in the model
(b) the response versus the time sequence in the gvlma object (gvlmaobj\$GlobalTest\$timeseq), which is the time sequence used for computing the directional test statistic S^24,
(c) the standardized residuals vs the fitted values,
(d) a histogram of the standardized residuals,
(e) a normal probability plot of the standardized residuals, and
(f) a plot of the standardized residuals versus the time sequence.

GVLMAを用いない場合と異なり、残差の分布を視覚的に把握することができる。加えて、検証結果と合わせて考察することで、視覚的な感覚に頼らず、残差を統計量として把握しやすくなる。
f:id:espio999:20191022214841p:plain

外れ値の検証、プロット

GVLMAオブジェクトをdeletion.gvlmaに引き渡し、一個抜き交差検証を実施することができる。この時の出力結果により、先の5つの検証結果毎に外れ値を把握することができる。StatがGlobal Statでの外れ値、Stat1~4がSkewnessからHeteroscedasticityの検証に該当する。
さらに出力オブジェクトをプロットすることで、外れ値の影響力を視覚的に把握することができる。
検証結果に「Assumptions NOT satisfied!」が表示された項目があったならば、該当項目の外れ値をデータから除外してみる、という検討が容易になる。
f:id:espio999:20191022214958p:plain
f:id:espio999:20191022215020p:plain