Quantile regression の理解のためには、quantile の最小値問題による
定義の理解 が必要である。
データ \( a = a_1,...,a_n \) が与えられたとき、\( t \)の関数を２つ考える.
\[ f(t) = \sum_{i=1}^n (a_i - t)^2 , \]
\[ g(t) = \sum_{i=1}^n |a_i - t| \]
これらの関数を最小にする \( t \) の値は何であろう。\( f \)については、\( t \)で
微分すれば直ちに解るように、\( t = mean(a) \) において最小となる。
最も簡単な最小二乗法である。
では、\( g \) ではどうなるだろうか。実は\( t = median(a) \)において最小となる。
これを実験してみよう。Poisson分布の乱数を発生させ、\( y=f(t), \, y=g(t) \)
のグラフを描く(左図）。これから理解されると期待する。
但し、データ数が偶数の場合、medianは中央の２つの値の平均とするが、
\( g \) はその２点の間で一定値となるので、どこを選択するかは自由度がある。
\( g \)の最小値を与える\( t \)の特徴付けを考えてみるとよい。
また、1st quartile(\( a_{1st} \)), 3rd quartile(\( a_{3rd} \)) も実は関数を設定すればその最小化として実現できる（右図）。
\[ h_1(t) = \frac{1}{4}\sum_{a_i >= t} |a_i - t| + \frac{3}{4}\sum_{a_i < t} |a_i - t| \] \[ h_3(t) = \frac{3}{4}\sum_{a_i >= t} |a_i - t| + \frac{1}{4}\sum_{a_i < t} |a_i - t| \]
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = TRUE, fig.width = 10)
options(width = 116, scipen = 3) # 1000000 (5+(3-1))桁
set.seed(1234)
opar <- par(mfrow=c(1,2), mar=c(4,4,3,1), family="serif")
a <- c(rnorm(20,3),rnorm(30,5, 2))
xdat <- seq(0,max(a),length=200)
fdat <- sapply(xdat, function(x) sum((a-x)^2) )
gdat <- sapply(xdat, function(x) sum(abs(a-x)) )
#dat <- cbind(b=c(a,a),val=c(fdat,gdat),code=rep(c("f","g"),each=20))
dat <- cbind(a=a,fdat=fdat,gdat=gdat)
fgdat <- c(fdat,gdat)
plot(0,0,type="n", xlim=c(0,max(a)),ylim=c(0,500),
xlab='seq(0,max(a),length=200)', ylab='f & g',
main='Mean and Median by L2 and L1 theories')
lines(xdat, fdat, lwd=2)
lines(xdat, gdat, lwd=2, lty=2, col=2)
meda <- median(a)
abline(v=mean(a), lwd=2)
abline(v=meda, col=2, lty=3, lwd=2)
abline(h=0)
rug(a)
boxplot(a, horizontal=TRUE, add=TRUE, notch=TRUE, at=1000/3, boxwex=100)
legend("topleft", c("mean","median"), text.col=1:2, col=1:2, lty=1:2, lwd=2)
a1 <- quantile(a, 0.25)
a3 <- quantile(a, 0.75)
auh <- median(a[a>=meda])
alh <- median(a[a<=meda])
h1dat <- sapply(xdat, function(x) 3*sum(abs(a[a<x]-x))/4+
sum(abs(a[a>=x]-x))/4)
h3dat <- sapply(xdat, function(x) sum(abs(a[a<x]-x))/4+
3*sum(abs(a[a>=x]-x))/4)
h13dat <- c(h1dat, h3dat)
plot(0,0,type="n", xlim=c(0,max(a)),ylim=c(0,150),
xlab='seq(0,max(a),length=200)', ylab='h1 & h3',
main='Q1 and Q3 by L1 theory')
lines(xdat, h1dat, lwd=2, col=3)
lines(xdat, h3dat, lwd=2, lty=2, col=4)
abline(v=a1, col=3, lwd=2)
#abline(v=alh, col=5, lwd=2)
abline(v=a3, col=4, lty=3, lwd=2)
#abline(v=auh, col=6, lty=3, lwd=2)
abline(h=0)
rug(a)
boxplot(a, horizontal=TRUE, add=TRUE, notch=TRUE, at=100, boxwex=30)
legend("topleft", c("1st Q","3rd Q"), text.col=3:4, col=3:4, lty=1:2, lwd=2)
一般に、\( q \)-quantile \( (0 < q < 1) \) は、次式の最小値を与える \( t \)である。 特に上記の \( g \) は定数倍を除き \( h_{1/2} \) である。
\[ h_q(t) = q \sum_{a_i >= t}(a_i-t) + (1-q) \sum_{a_i < t}(t-a_i)= q \sum_{i=1}^n (a_i-t) + \sum_{a_i < t}(t-a_i) \]
すると、通常の２乗和\( f \)の1乗版である\( g \)が、さらに一般化されて\( h_q \)になる
ということであり、通常の最小二乗法で行っているパラメータ推定を、
\( h_a \)の類似によって行うことができる可能性が示唆される。実際に通常の
線形回帰分析に対応して、分位回帰(Quantile Regression)分析を行うことが
できる。上記のことから推測できるように、通常の回帰分析が平均的な位置の
直線（超平面）を求めるのに対して、q-分位点的な位置を求めるものと
想定される。
たとえば、時間的に測定されたデータ列があるとする。平均が
各時点のデータにおいてほぼ0であるとすれば、回帰式はほぼy=0である。
しかし、分散が時間経過に従い大きくなっているとすれば、3rd Quartile点は
時間経過とともに大きくなっているであろう。そしてそれが、分位回帰に
よってわかるということになる。このとき、回帰分析の残差グラフは
時間の経過に従い残差が大きくなるという傾向を示しているであろう。
set.seed(123)
tseq <- list()
for (cnt in 1:10){
tseq[[cnt]] <- rnorm(20, 0, 10*cnt)
# tseq[[cnt]] <- tseq[[cnt]] - mean(tseq[[cnt]])
}
tdat <- cbind.data.frame(time=rep(1:10,each=20), dat=unlist(tseq))
par(mfrow=c(2,3), family="serif")
plot(tdat)
plot(density(tdat[,2]))
rug(tdat[,2])
tlm <- lm(dat ~ time, data=tdat)
for (cnt in c(1:3,5)) plot(tlm, which=cnt, cex.id=1.5
,cook.levels=c(4/200, 2.5^2/200, 3^2/200, 0.5,1)) #
Quantile Regression の生態学での適用例をみれば、実際の必要性が解る。
Quantile Regression をパッケージ quantreg の rq() で行い、
Quantile値 \( \tau=0.1,0.25,0.5,0.75,0.9 \) の設定、即ち decile の0.1, 0.9 と
1st quartile, median, 3rd quartile の回帰を行うと次の通りである。
点線が通常の線形回帰であり、それに比較して傾向を捉えていることが
わかる。なお、plot(
パッケージ VGAM には、LMS quantile regression の関数 lmsqreg
がある。
Random Forests で回帰を行うことに対応して、Quantile Regression Forests
という手法も開発されており、パッケージ quantregForests で実行できる。
ただし、quantiles は 0.1, 0.5. 0.9 の３点に固定されている。
また、Quantile regression は線形計画法による変数選択、Dantzig selector
の実現にも使うことができる。そもそも、Quantile Regression そのものが
実験計画法として定式化できる。Dantzig selector には色々な意見があり
特に Lasso, Lars より使うべきものともされていない。Outlier の影響を
受けやすいということもある。
Lasso系 は quantreg においても、lassoとSCAD(smoothly clipped absolute deviation penalty)を実行できる。パラメータ設定において method=“lasso”
“scad"とすればよい。但し、パッケージ lasso, lars の方が、種々の出力を
得るのに便利である。また、flare(Family of Lasso Regression)は
Lasso, Dantzig Selector, LAD Lasso, SQRT Lasso, Lq Lasso と、sparse な
行列に向けた TIGER, CLIME using L1 をサポートしている。
Q-Q plot の命令 qqplot(u, v)は、データ列uとvを quantile に変換して
比較した散布図を描く。整合的に対応をとるので\( length(u) \ne length(v) \)であっても
問題ない。同じ種類の分布であればほぼ直線上に並ぶ。その形状は
qu <- quantile(u, seq(0,1,0.01)) ; qv <- quantile(v, seq(0,1,0.01))
plot(qu, qv)
とすれば全容がわかる（ここでは100点とした）。特に、uを標準正規分布として、x軸を理論的なquantile としたものを、Normal Q-Q plot とよび、正規分布
との一致度を見るのに用いる。qqnorm(v); qqline(y) と用いる。
この場合、qqnorm(v)は、z=ppoints(length(n)) によって定まる
確率点に対し plot(qnorm(z),sort(v)) としている。またqqplot(v)は
vの1st Quartile と 3rd Quartileの点を結ぶ直線を描く。
パッケージ car の qqPlot は、更に情報を添加する。理論分布は
dist="norm”, “unif” などと指定すればよい。またはずれが大きいデータを
確認するためにラベルを付与する個数を、id.n=5 のように指定する。
コンソールには、付加したラベルと、そのデータの昇順の順位が
出力される。
赤の点線（色はcol.linesで指定可能）は引き渡したオブジェクトによって
異なる。単なるデータ数列の場合はline=“quartile” 即ち四分位数のペア
(1st, 3rd)となり、クラスlmの場合は line=“robust” であり、robust
regrassion(MASS::rlm を用いる)による残差回帰線である。なおこの
クラスlmの場合、dist=“t"、つまりdefaultで理論分布はt分布が採用される。
require(quantreg)
require(car)
qrg <- rq(dat ~ time, tau=c(0.1, 0.25,0.5,0.75, 0.9), data=tdat)
summary(qrg)
# plot.rqs(qrg, mfrow=c(1,2)) ## HTML 変換時に異常が起こるので省略
par(mfrow=c(3,2), family="serif")
plot(tdat, main="Quantile regression. tau=0.1, 0.25, 0.5, 0.75, 0.9")
for (cnt in 1:5) abline(coef(qrg)[,cnt], col=cnt, lwd=2)
abline(tlm, lty=2, col=6)
legend("topleft",c("0.1","1stQR","medianR","3rdQR","0.9", "LM"),
lty=c(rep(1,5),2), lwd=c(rep(2,5),1), col=1:6)
set.seed(123)
u <- rnorm(20)
v <- runif(30)
qu <- quantile(u, seq(0,1,0.01))
qv <- quantile(v, seq(0,1,0.01))
qqnorm(v, sub="vs runif(30)") ; qqline(v)
qqplot(u, v, main="qqplot(rnorm(20),runif(30))")
plot(qu, qv, main="plot(qu, qv) 100pts", pch=20)
par(new=TRUE)
qqplot(u, v, col=2, cex=1.5, xlab="", ylab="")
## Anderson-Darling検定は略
qqPlot(residuals(qrg), id.n=5, main='qqPlot(residuals(qrg), id.n=5)')
qqPlot(v, dist="unif", id.n=5, main='qqPlot(v, dist="unif", id.n=5, pch=16, cex=2, col=8)',
pch=16, cex=2, col=8)
par(mfrow=c(1,1), mar=c(5,4,3,2))
次の図は通常の環境では出せないので、コマンドは省略する。
上三角部は相関係数（上）と偏相関係数（下）、
対角部は種別の標本密度分布（全体で１となるように、件数による重み処理済）。
前節で car::qqPlot を用いた。carにはlm, glm オブジェクトに対する
診断プロットがいくつかある。ここでそのいくつかを紹介する。
まずは、通常の plot.lm による６種類の診断プロットを掲げる。
その次の６図のうち４番目の influencePlotは、はずれ値、影響の強いデータを判別するために
ハット行列 \( H = \tilde X (t(\tilde X) \tilde X)^{-1} t(\tilde X) \)
の対角線成分 \( h_{ii} \) を\( x \)軸に、スチューデント化残差を\( y \)軸に、
そしてクックの距離を
円の半径で表示している。極端な元を検出することが目的である。
\[ tr(H)=tr((t(\tilde X) \tilde X)^{-1} t(\tilde X)\tilde X )=tr(I_{p+1})=p+1 \]
より、Hat Value \( h_{ii} \)の平均は\( (p+1)/n \)となる。よって、その２倍を
こえるものが大きいとし、nとpが小さい時には３倍をこえる大きいとする。
即ち閾値を\( (2p+2)/n, \, (3p+3)/n \)とする。
下記の実例では
\( n=150, p=2 \) であるから、\( 2 \times (2+1)/150 = 0.04 \)となる。
その位置と３倍の\( 0.06 \)に破線が立てられている。これらを超えると回帰に
大きな影響を与えるとみなす。但し、後述の\( h_{ii}-1/n \geq 0 \)を考慮すると
実質的な２倍、３倍の過剰閾値は\( (2p+1)/n, \, (3p+1)/n \) である。
ここで、Hの対角線成分について \( 0 < h_{ii} < 1 \) を示す。まず容易に
解るように、\( H \) は \( H^2 = H, t(H)=H \) を満たす（即ち、射影子である）。
今\( (p+1)\times (p+1) \)直交行列\( U \)を取り、\( Z=HU \) を考える。このとき、
\[ Zt(Z) = HUt(U)t(H)=H I_n H = H \]
となる、対角第\( i \)成分を比較すれば
\[ \sum_{j=1}^n z_{ij}^2 = h_{ii} \]
従って、\( h_{ii} > 0 \). また \( I-H \)も \( (I-H)^2=I-H, t(I-H)=I-H \)を
満たすので、同様に\( 1-h_{ii} > 0 \) を得る。Q.E.D.
実はより強く \( h_{ii} \geq 0 \) が成立する。ここで等号成立は
データ\( x_i \)が、\( x \)の平均値\( \bar x \)に一致するとき、その時に限る。
さらに、\( h_{ii} - 1/n \) 、あるいは \( (h_{ii} - 1/n)/(1-h_{ii}) \) は、
平均\( \bar x \) 、あるいは \( x_i \)を除いた\( n-1 \)個のデータにおける平均\( \bar x_{(-i)} \)から\( x_i \)へのマハラノビス汎距離の２乗に比例する。
\[ h_{ii}-\frac{1}{n} = \frac{D_n(x_i; \bar x)^2}{n-1}, \,\,
\frac{h_{ii} - 1/n}{1 - h_{ii}} =
\frac{(n-1)D_{n-1}(x_i; \bar x_{(-i)})^2}{n(n-2)}. \]
回帰の残差推定量 \( \hat e_i \) の分散の期待値はモデルでの\( e_i \)の分散
\( \sigma^2 \)ではなく、\( (1-h_{ii}) \sigma^2 \)になることが解る。
\( \sigma^2 \)の推定値\( s_e^2 \)は、残差平均平方和\( S_e/(n-p-1) \)である。
\[ r_i = \frac{\hat e_i}{s_e \sqrt{1-h_{ii}} } \]
をスチューデント化残差(Studentized Residual)とよぶ。\( r_i \)は\( t \)分布、
\( t(n-p-1) \) に従う。標準化残差と同様に、\( r_i \)は
絶対値で２, あるいは３を超えると、Outlierであるとみなす。正確には、
\( qt(0.975, n-p-1), qt(0.995, n-p-1) \) を境界とすればよい。
なお、計算の詳細は次の通り。
クラス lm (glm) に対するメソッドとして、標準の stats に rstudent,
cooks.distance(=\( CookD^2 \))があり、内部で関数 lm.influenceを呼び出している。メソッドの内容は getS3method("rstudent”, “lm”) により閲覧できる。
infl = lm.influence(model) とすれば、
\( s_e \) = infl$sigma, Hat = diag(H) = infl$hat, \( \hat e \) = infl$wt.res
となっている。従って、
\[ StudRes = rstudent(model) = \frac{\hat e}{s_e \sqrt{1 - Hat}} \]
\[ CookD^2 = cooks.distance(model) = \frac{1}{p+1} StudRes^2 \frac{Hat}{1-Hat} \]
である。これより、Cook's distance は x=Hat, y=StudRes 座標でいえば
次の等高線を描けばよい。
\[ y^2 = CookD^2 (p+1)(\frac{1}{x}-1) \]
plot.lm(obj, which=5)での赤破線がそれである。
一方、which=6 で描かれる 'Cook's distance vs leverage/(1-leverage) plot'
は、\( x=Hat/(1-Hat), y=CookD^2 \, \) のグラフであり、その傾斜\( StudRes^2/(p+1) \)
が黒破線で入っている。
なお、次元から言って CookD を'クックの距離'と
呼ぶべきであるが、その２乗の cooks.distance を 'クックの距離' と呼ぶことも多いので注意を要する。
“Kenneth A. Bollen & Robert W. Jackman Regression Diagnostics: An Expository Treatment of Outliers and Influencial cases”
\( Modern Methods of Data Analysis \) (1990, p.266, p.268 Table 6.1)
では、\( CookD > 2/\sqrt{n} \) という Belsley et al.(1980, p.28)に基づく
基準を提案している。n=150 ならば閾値は0.1633(2乗では0.0267) となる。
僅かな例ではあるが経験上では\( CookD > 2.5/\sqrt{n} \)の方が適切に思える。今回も、この２つと\( 3/\sqrt{n} \) の３本の線を試みに入れている。
一方、nにかかわらず、\( CookD^2 > 1 (or 0.5) \) という意見もあるようである。
ie.lm <- lm(Sepal.Length ~ Sepal.Width + Petal.Width, data=iris)
#ie.lm <- lm(Kakaku ~ Tochi + Toho, data=IeKakaku)
par(mfrow=c(3,2), mar=c(4,4,2,1.5),family="serif")
for (cnt in 1:6) plot(ie.lm, which=cnt, col=iris[,5], cex.id=1.5,
id.n=5, cook.levels=c(0.0267, 0.0417, 0.06,0.5,1.0)) #
# 4/150=0.026667, 2.5^2/150=0.0417, 3^2/150=0.06
abline(h=0.0267, lty=2, col=2)
abline(h=0.0417, lty=2, col=2)
abline(h=0.06, lty=2, col=2)
residualPlots(ie.lm, col=iris$Species, layout=NA, main="residualPlots",
id.n=3, id.cex=1.5, id.col="magenta")
Test stat Pr(>|t|)
Sepal.Width 1.004 0.317
Petal.Width -3.340 0.001
Tukey test -3.116 0.002
source("influencPlot2.R")
influencePlot2(ie.lm, col=iris$Species, main="influencePlot",
id.n=3, id.cex=1.5, id.col="blue")
StudRes Hat CookD
16 0.2244 0.07040 0.03577
61 -0.5181 0.05640 0.07332
101 -2.0662 0.03703 0.23136
107 -2.7617 0.01792 0.21066
115 -2.5304 0.02332 0.22171
118 1.3436 0.05522 0.18703
123 2.6927 0.01440 0.18404
132 2.2619 0.04788 0.28882
#ceresPlots(ie.lm, col=iris$Species, layout=NA)
crPlots(ie.lm, col=iris$Species, layout=NA, # main="crPlots",
id.n=3, id.cex=1.5, id.col="magenta")
#leveragePlots(ie.lm, col=iris$Species, layout=NA)
avPlots(ie.lm, col=iris$Species, layout=c(1,2), main="avPlots",
id.n=3, id.cex=1.5, id.col="magenta")
Bootstrap により係数の変動を検討することができる。やはり car に 含まれる car::Boot を用いる。
par(mfrow=c(1,2))
ie.boot <- car::Boot(ie.lm, R=2000)
summary(ie.boot)
R original bootBias bootSE bootMed
(Intercept) 2000 3.457 -0.003959 0.3046 3.451
Sepal.Width 2000 0.399 0.000709 0.0879 0.401
Petal.Width 2000 0.972 0.001870 0.0565 0.973
confint(ie.boot)
Bootstrap quantiles, type = bca
2.5 % 97.5 %
(Intercept) 2.8538 4.0719
Sepal.Width 0.2226 0.5715
Petal.Width 0.8614 1.0865
plot(ie.boot)
par(mfrow=c(2,2))
hist(ie.boot, layout=NA)
dataEllipse(ie.boot$t[,2], ie.boot$t[,3],
xlab="Sepal.Width", ylab="Petal.Width",
cex=.5,levels=c(.5,.95,.99),robust=TRUE)
par(mfrow=c(1,1),mar=c(5,4,3,2))
# MASS::stepAIC(ie.lm)