#======================================================
#	5 ロバスト回帰
#======================================================
library(robust) #パッケージのロード

#======================================================
#	5.1 ロバスト回帰の必要性
#======================================================
#------------<外れ値と線形回帰>------------

#図5.1
set.seed(123)
#xy <- rmvnorm(180,sigma=matrix(c(1, 0.8, 0.8, 1), nrow=2))　#通常値
xy <- mvrnorm(180,mu=c(0, 0), Sigma=matrix(c(1, 0.8, 0.8, 1), nrow=2))　#通常値
x <- c(rnorm(20)-2, xy[ ,1]+2)　#外れ値は最初の20
y <- c(rnorm(20)+2, xy[ ,2]-2)　#外れ値は最初の20
fit.lm <- lm(y~x)　　　　　　　 #一般的な線形モデル
fit.rob <- lmrob(y~x)           #ロバストな線形モデル
plot(x[1:20], y[1:20], pch=17, col="skyblue", 
	xlim=range(x), ylim=range(y), xlab="x", ylab="y") #外れ値のプロット
points(x[21:200], y[21:200], pch=16, col="springgreen")  #通常値
abline(fit.lm, lty=2, col="steelblue", lwd=2)　 #一般的な回帰直線
abline(fit.rob, col="darkgreen", lwd=2)　　　　 #ロバストな回帰直線
text(-2, 0, "LM", srt=-20)
text(-1, -4, "Robust", srt=45)

#------------<外れ値診断の限界>------------
#診断統計量
influence.measures(fit.lm)

#診断プロット
#図5.2
wagner.lm <- lm(y ~ ., data=wagnerGrowth)　  #一般的な線形モデル
wagner.rob <- lmRob(y ~ ., data=wagnerGrowth)#ロバストな線形モデル
models <- fit.models(wagner.lm, wagner.rob)  #モデル比較オブジェクト
windows(10, 7)                               #ウィンドウサイズを変える。
plot(models, which=8)　　　　　　　　　　　　#修正残差比較
plot(models, which=2)                        #QQプロット比較

#======================================================
#	5.2 ロバスト回帰に関連する各関数
#======================================================
#------------<bisquare>------------

#図5.3
x <- seq(-5, 5, length=200)

op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))

y <-Mchi(x, cc=4.685, psi="bisquare")
plot(x, y, type="l", ylab="rho", main=expression(rho(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mpsi(x, cc=4.685, psi="bisquare")
plot(x, y, type="l", ylab="psi", main=expression(psi(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mpsi(x, cc=4.685, psi="bisquare", deriv=1)
plot(x, y, type="l", ylab="psi", main=expression({psi^minute}(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mwgt(x, cc=4.685, psi="bisquare")
plot(x, y, type="l", ylab="weight", main=expression(W(x)==psi(x)/x), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

par(op)

#------------<極小値>------------

telef
cbind(telef[1:12, ], telef[13:24, ])

#psiの解を求める関数
psi.mean <- function(mu, data) {
	madn <- mad(data)
	n <- length(mu)
	n2 <- length(data)
	scaled.data <- numeric(n2)
	psi.data <- numeric(n)
	ans <- numeric(n)
	for(i in 1:n) {
		scaled.data <- (data - mu[i])/madn
		psi.data <- Mpsi(scaled.data, cc=4.685, psi="bisquare")
		ans[i] <- mean(psi.data)
	}
	ans
}

#rhoの極小値と極大値
rho.mean <- function(mu, data) {
	madn <- mad(data)
	n <- length(mu)
	n2 <- length(data)
	scaled.data <- numeric(n2)
	rho.data <- numeric(n)
	ans <- numeric(n)
	for(i in 1:n) {
		scaled.data <- (data - mu[i])/madn
		rho.data <- Mchi(scaled.data, cc=4.685, psi="bisquare")
		ans[i] <- mean(rho.data)
	}
	ans
}

#図5.4　複数の極小値
mu <- seq(-8, 28, length=200)
y <- psi.mean(mu=mu, data=telef$Calls)
z <- rho.mean(mu=mu, data=telef$Calls)

op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 2.1, 1.1))
plot(mu, y, type="l", xlab=expression(mu), 
	ylab=expression(paste("mean(", psi, ")", sep="")))
abline(h=0, lty=2)
text(c(2, 7, 14), rep(-0.1, 3), c("a", "b", "c"))
plot(mu, z, type="l", xlab=expression(mu), 
	ylab=expression(paste("mean(", rho, ")", sep="")))
text(c(1.5, 8, 14), c(0.35, 1, 0.8), c("a", "b", "c"))
par(op)

#======================================================
#	5.3 MM推定
#======================================================
#------------<S推定>------------
set.seed(123)
S.est <- lmrob.S(cbind(1, telef$Year), telef$Calls, 
	control=lmrob.control(seed=1, nResamp=50, trace.lev=2))
S.est$scale
S.est$coef

#------------<MM推定>------------
set.seed(123)
MM.est <- lmrob(Calls ~ Year, data=telef, control=lmrob.control(seed=1))
MM.est$scale
MM.est$coef
summary(MM.est)
#追加
plot(telef)
abline(MM.est)

#------------<適用例>------------
#データ
library(mlbench)
data(BostonHousing)

#当初モデル
boston.lm0 <- lm(medv ~ . ^ 2, data=BostonHousing)

#stepAIC
boston.lm1 <- stepAIC(boston.lm0, trace=0)

#MM推定
set.seed(123)
boston.rob <- lmrob(formula(boston.lm1), data=BostonHousing, 
	control=lmrob.control(seed=1))
#警告メッセージが出るので、max.itを増やす。
set.seed(123)
boston.rob <- lmrob(formula(boston.lm1), data=BostonHousing, 
	control=lmrob.control(seed=1, max.it=100))
summary(boston.rob)

set.seed(123)
boston.rob2 <- update(boston.rob, .~. -dis -crim:nox -crim:rm -crim:dis
    -crim:ptratio -crim:b -crim:lstat -zn:tax -chas:ptratio -chas:lstat
    -nox:dis -rm:dis -rm:b -age:dis -age:tax -dis:rad -dis:tax -dis:lstat
    -rad:lstat -b:lstat)
summary(boston.rob2)

set.seed(123)
boston.rob3 <- update(boston.rob2, .~. -crim:zn -zn:dis -indus:rm
     -nox:lstat)
summary(boston.rob3)   

set.seed(123)
boston.rob4 <- update(boston.rob3, .~. -indus:nox)
summary(boston.rob4)

#比較のためのプロット
models <- fit.models(boston.lm1, boston.rob4) #モデル比較オブジェクト
windows(10, 7)                                #ウィンドウサイズを変える。
#図5.5
plot(models, which=2)　　　　　　　　　　　   #QQプロット
#図5.6
plot(models, which=10, lty=1:2)               #残差分布比較

#======================================================
#	5.4 KS2011とKS2014
#======================================================
#------------<LQQ>------------
#図5.7
x <- seq(-9, 9, length=200)

op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))

y <-Mchi(x, cc=c(-0.5,1.5,0.95,NA), psi="lqq")
plot(x, y, type="l", ylab="rho", main=expression(rho(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mpsi(x, cc=c(-0.5,1.5,0.95,NA), psi="lqq")
plot(x, y, type="l", ylab="psi", main=expression(psi(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mpsi(x, cc=c(-0.5,1.5,0.95,NA), psi="lqq", deriv=1)
plot(x, y, type="l", ylab="psi", main=expression({psi^minute}(x)), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

y <-Mwgt(x, cc=c(-0.5,1.5,0.95,NA), psi="lqq")
plot(x, y, type="l", ylab="weight", main=expression(W(x)==psi(x)/x), col="steelblue")
abline(h=0, v=0, col="gray60", lty=2)

par(op)

#------------<適用例>------------

set.seed(123)
KS.est <- lmrob(Calls ~ Year, data=telef, 
	control=lmrob.control("KS2011", seed=1))
KS.est$scale    #KS2011によるスケール推定値
KS.est$coef     #KS2011による回帰係数推定値
#再掲
MM.est$scale    #MM推定によるスケール推定値
MM.est$coef     #MM推定による回帰係数推定値
#追加
summary(KS.est)
windows()
plot(telef)
abline(MM.est)
abline(KS.est, lty=2)

#KS推定
set.seed(123)
boston.ks <- lmrob(formula(boston.lm1), data=BostonHousing, 
	control=lmrob.control("KS2014", seed=1))
summary(boston.ks)

set.seed(123)
boston.ks2 <- update(boston.ks, .~. -crim:dis -crim:ptratio 
    -chas:lstat -nox:dis -nox:lstat -rm:dis -rm:b -age:dis 
    -age:tax -dis:rad -dis:tax -rad:lstat -b:lstat)
summary(boston.ks2)

#比較のためのプロット
models <- fit.models(boston.rob, boston.ks2) #モデル比較オブジェクト
windows(10, 7)                             #ウィンドウサイズを変える。
plot(models, which=2)　　　　　　　　　　　#QQプロット
plot(models, which=10, lty=1:2)            #残差分布比較

