#======================================================
#	4 線形回帰と平滑化
#======================================================

#======================================================
#	4.1 線形回帰モデル
#======================================================
#------------<最小二乗法>------------

青森県 <- 市区町村[市区町村$都道府県=="青森県", ]  #青森県のデータを抽出

#単回帰

青森.lm <- lm(log(高齢者人口比率)~log(人口密度), data=青森県)

#図4.1
plot(log(高齢者人口比率) ~ log(人口密度), data=青森県, pch=16, col="steelblue")
abline(青森.lm)
segments(x0=log(青森県$人口密度), x1=log(青森県$人口密度),
	 y0=log(青森県$高齢者人口比率), y1=fitted(青森.lm), lty=2)

#------------<線形回帰からの出力>------------

#lmオブジェクト、summary.lmオブジェクト
names(青森.lm)
names(summary(青森.lm))

#回帰結果の要約
summary(青森.lm)

#------------<QR分解による解法>------------

#QR分解
(X <- model.matrix(青森.lm))
QR <- qr(X)
Q <- qr.Q(QR, complete=T)
(R <- qr.R(QR))

#βの推定
y <- log(青森県$高齢者人口比率)
tQy <- qr.qty(QR, y) 
p <- dim(R)[2]
Qy.f <- tQy[1:p]
(beta.hat <- solve(R) %*% Qy.f)

#Std.Errorの推定
n <- length(y)                  #標本数
r <- tQy[(p + 1):n]             #残差
resvar <- t(r) %*% r / (n - p)  #回帰式の分散
(stddev <- sqrt(resvar))        #回帰式の標準誤差
R.inv <- chol2inv(QR$qr[1:p, 1:p, drop=FALSE]) #逆行列
(se <- sqrt(diag(R.inv)*resvar)) #回帰係数の標準誤差

#t値
(t.value <- beta.hat/se)

#p値
pt(abs(t.value), n-p, lower.tail=F) * 2 

#決定係数
f <- X %*% beta.hat #当てはめ値
rss <- sum(r ^2)    #残差平方和
mss <- sum((f - mean(f))^2)  #回帰平方和
(r.squared <- mss/(mss + rss)) #決定係数
#自由度修正済決定係数
(adj.r.squared <- 1 - (1 - r.squared)*((n-1)/(n-p))) 

#F値
(f.value <- (mss/(p-1))/resvar)
pf(f.value, p-1, n-p, lower.tail=F) #p値

#------------<回帰診断：診断統計量による方法>------------
#------------<様々な外れ値>------------

#図4.2
op <- par(mfrow=c(2, 2), mar=par()$mar-0.4)
#パターン１
set.seed(123)
x <- 1:10 * 0.5 + rnorm(10)
y.o <- y <- 1:10 * 0.5 + rnorm(10)
x[5] <- mean(x)
y.o[5] <- 9.5
plot(x[-5], y.o[-5], xlab="x", ylab="y", ylim=c(0, 10),
	main="パターン１") 
points(x[5], y.o[5], pch=16, col="steelblue")
abline(lm(y.o ~ x))
abline(lm(y[-5] ~ x[-5]), lty=2)

#パターン２
x.o <- x <- 1:10 * 0.5 + rnorm(10)
y.o <- y <- 1:10 * 0.5 + rnorm(10)
x.o[10] <- 9.5
y.o[10] <- 7.3
plot(x[-10], y.o[-10], xlab="x", ylab="y", xlim=c(0, 10), ylim=c(0, 10)
	,main="パターン２")
points(x.o[10], y.o[10], pch=16, col="steelblue")
abline(lm(y.o ~ x.o))
abline(lm(y[-10] ~ x[-10]), lty=2)

#パターン３
x.o <- x <- 1:10 * 0.5 + rnorm(10)
y.o <- y <- 1:10 * 0.5 + rnorm(10)
x.o[10] <- 9.5
y.o[10] <- 12
plot(x[-10], y.o[-10], xlab="x", ylab="y", xlim=c(0, 10), ylim=c(0, 12)
	,main="パターン３")
points(x.o[10], y.o[10], pch=16, col="steelblue")
abline(lm(y.o ~ x.o))
abline(lm(y[-10] ~ x[-10]), lty=2)

par(op)

#------------<各種診断統計量>------------
#計算結果
influence.measures(青森.lm)

#計算の基礎となるオブジェクト
model <- 青森.lm
infl <- influence(model)
p <- model$rank #パラメータの数
e <- weighted.residuals(model)                       #残差
s <- sqrt(sum(e^2, na.rm = TRUE)/df.residual(model)) #残差の標準誤差
mqr <- qr(model.matrix(model))                       #QR分解
xxi <- chol2inv(mqr$qr, mqr$rank)　      #コレスキー分解による逆行列
si <- infl$sigma　#i番目の観測値を除去した場合のモデル残差の標準偏差
h <- infl$hat　   #ハット行列の対角要素

#dfbatas
dfbetas <- infl$coefficients/outer(infl$sigma, sqrt(diag(xxi)))
dfbetas

#dffits
dffits <- e * sqrt(h)/(si * (1 - h))
dffits

#cov.ratio
cov.ratio <- (si/s)^(2 * p)/(1 - h)
cov.ratio

#Cookの距離
cooks.d <- ((e/(s * (1 - h)))^2 * h)/p
cooks.d

#梃子比（ハット値）
h 

#抽出
ans <- influence.measures(青森.lm)
(out <- which(apply(ans$is.inf, 1, any)))

#推定から除外
青森.lm2 <- lm(log(高齢者人口比率)~log(人口密度), data=青森県[-out, ])
summary(青森.lm2)

#図4.3
plot(log(高齢者人口比率)~log(人口密度), data=青森県, 
        main="外れ値除去前後の回帰直線",xlim=c(1.5, 7))
abline(青森.lm2)
abline(青森.lm, lty=2)
legend("topright", legend=c("外れ値除去前" ,"外れ値除去後"), lty=2:1)
text(log(青森県$人口密度[out]), log(青森県$高齢者人口比率[out])+0.03, 
        青森県$市区町村[out])

#------------<回帰診断：図示による方法>------------
#------------<既定のプロット>------------

#図4.4
op <- par(mfrow=c(2, 2))
plot(青森.lm)
par(op)

#------------<パッケージcarにより追加されたプロット>------------
library(car)

#QQプロット
#図4.5
qqPlot(青森.lm, id.n=3, labels=青森県$市区町村)

#dfbetas
#図4.6
dfbetasPlots(青森.lm, id.n=4, labels=青森県$市区町村)

#------------<回帰当てはめ結果による予測>------------

#関数fitted
青森.fit <- fitted(青森.lm) #当てはめ結果を取得
plot(log(高齢者人口比率)~log(人口密度), data=青森県) #データ点
lines(log(青森県$人口密度), 青森.fit) #回帰直線

#同じ
plot(log(高齢者人口比率)~log(人口密度), data=青森県)
abline(青森.lm)

#predictによる信頼区間
#図4.7
rng <- range(青森県$人口密度)
xx <- seq(rng[1], rng[2], length=100) 
青森.conf <- predict(青森.lm, newdata=data.frame(人口密度=xx), 
	interval="confidence")
plot(log(高齢者人口比率)~log(人口密度), data=青森県, main="信頼区間")
matplot(log(xx), 青森.conf, type="l", lty=c(1, 2, 2), add=T) #信頼区間

#predictによる予測区間
#図4.8
青森.pred <- predict(青森.lm, newdata=data.frame(人口密度=xx), 
	interval="prediction")
plot(log(高齢者人口比率)~log(人口密度), data=青森県, main="予測区間")
matplot(log(xx), 青森.pred, type="l", lty=c(1, 2, 2), add=T) #予測区間

#======================================================
#	4.2 重回帰モデルへの拡張
#======================================================
#------------<モデル式>------------
#------------<交互作用>------------

data(car.test.frame, package="rpart")       #データをロードする。 

Japan <- car.test.frame$Country == "Japan"  #日本車かどうか
#交互作用項を含めない場合
fit.0 <- lm(Mileage ~ Disp. + Japan, data=car.test.frame)
summary(fit.0)

#当てはめ結果
japan.mpg <- fitted(fit.0)[Japan==T]
japan.ord <- order(japan.mpg)
non.japan.mpg <- fitted(fit.0)[Japan==F]
non.japan.ord <- order(non.japan.mpg)

#図4.9
#プロット
plot(car.test.frame$Disp.[Japan==T][japan.ord], 
        japan.mpg[japan.ord], type="l", 
        xlim=c(70, 310), ylim=c(15, 35), xlab="排気量", ylab="燃費",
　　    main="燃費：交互作用項なし")
lines(car.test.frame$Disp.[Japan==F][non.japan.ord], 
        non.japan.mpg[non.japan.ord], lty=2)
legend("topright", legend=c("日本車", "日本車以外"), lty=1:2)

#交互作用項を含める場合
fit.1 <- lm(Mileage ~ Disp. + Japan + Disp.:Japan, data=car.test.frame)
summary(fit.1)

#当てはめ結果
japan.mpg <- fitted(fit.1)[Japan==T]
japan.ord <- order(japan.mpg)
non.japan.mpg <- fitted(fit.1)[Japan==F]
non.japan.ord <- order(non.japan.mpg)

#図4.10
#プロット
plot(car.test.frame$Disp.[Japan==T][japan.ord], 
        japan.mpg[japan.ord], type="l", 
        xlim=c(70, 310), ylim=c(15, 35), xlab="排気量", ylab="燃費",
　　    main="燃費：交互作用項あり")
lines(car.test.frame$Disp.[Japan==F][non.japan.ord], 
        non.japan.mpg[non.japan.ord], lty=2)
legend("topright", legend=c("日本車", "日本車以外"), lty=1:2)

summary(car.test.frame$Disp.[Japan==T])

#------------<多項式>------------

#poly
p <- contr.poly(n=5)　#選点直交多項式を計算
round(p, 2)           #係数を小数点第２位まで丸める。
#上記の.L から^4 までは、１次から４次までの多項式を表す。
zapsmall(t(p) %*% p)　#お互いに直交している。

#ethanolの多項式回帰
data(ethanol, package="lattice")　#データをロードする。

poly.lm <- lm(NOx ~ poly(E, degree=4), data=ethanol)
summary(poly.lm)
ord <- order(ethanol$E)
#図4.11
plot(NOx ~ E, data=ethanol, xlab="E", ylab="NOx", pch=16, col="steelblue")
lines(ethanol$E[ord], predict(poly.lm)[ord])

new.x <- seq(from=0.4, to=1.4, by=0.01)
#図4.12
plot(NOx ~ E, data=ethanol, xlab="E", ylab="NOx", 
        pch=16, col="steelblue", xlim=c(0.4, 1.4))
lines(new.x, predict(poly.lm, newdata=data.frame(E=new.x)))

#------------<応答変数の変換>------------
#------------<変換の例>------------

#図4.13
x<-seq(0.01,10,len=100)
op <- par(mfrow=c(2,2))
plot(x, sqrt(x), type="l", lwd=2, 
        main=expression(tilde(y)==y^2))
text(5, 2.5, expression(y==sqrt(x)), srt=40)

plot(x, log(x), type="l", lwd=2, ylim=c(-2,2.5), 
        main=expression(tilde(y)==exp(y)))
text(5, 1.9, expression(y==log(x)), srt=20)

plot(x, x^2, type="l", lwd=2, main=expression(tilde(y)==sqrt(y)))
text(5, 35, expression(y==x^2), srt=45)

plot(x, exp(x), type="l", lwd=2, main=expression(tilde(y)==log(y)))
text(5, 1000, expression(y==exp(x)), srt=5)
par(op)

#------------<ボックス・コックス変換>------------

#数値例
set.seed(123)
x0 <- seq(from=10, to=60, by=0.1)
x1 <- 2 * x0
x2 <- 3 * x0
y <- (x1 + x2 + rnorm(length(x0)))^0.333
fit.0 <- lm(y ~ x1 + x2)
#図4.14
plot(fit.0, which=2)

#boxcox
library(MASS)
#図4.15
lambda <- boxcox(y ~ x1 + x2, lambda = seq(2.5, 3.5, 1/10))
(lambda.max <- lambda$x[which.max(lambda$y)])　#最大値
y.tilde <- y^3

fit <- lm(y.tilde ~ x1 + x2)
plot(fit, which=2)

#------------<変数選択>------------
#------------<例：Boston住宅価格>------------

#データ
library(mlbench)
data(BostonHousing)

#当てはめ
fit.0 <- lm(medv ~ ., data=BostonHousing)
summary(fit.0)

fit.1 <- update(fit.0, . ~ . - age)
summary(fit.1)

fit.2 <- update(fit.1, . ~ . - indus)
summary(fit.2)

#逸脱度
deviance(fit.0)
#検算
res <- resid(fit.0)  #残差
(dev <- sum(res^2))  #残差平方和（RSS） 

#対数尤度
logLik(fit.0)
#検算
p <- fit0$rank #パラメータ数
N <- length(res)　#観測数
(LL <- - N * (log(2 * pi) + 1 - log(N) + log(dev))/2)　#最大対数尤度

#AIC
AIC(fit0)
#検算
-2 *LL + 2 * (1 + p) #AIC

#交互作用有り
fit.3 <- lm(medv ~ . ^ 2, data=BostonHousing)
summary(fit.3)

#stepAIC
fit.best <- stepAIC(fit.3)
summary(fit.best)

#------------<重回帰における回帰診断>------------
#------------<既定のプロット>------------
#図4.16
op <- par(mfrow=c(2, 2))
plot(fit.best)
par(op)

#------------<パッケージcarにより追加されたプロット>------------
library(car)

#residualPlots
#図4.17
residualPlots(fit.best, ~ crim + zn + indus + chas + nox + rm,
	layout=c(2, 3), fitted=F, main="予測変数ごとの残差(a)")

#図4.18
residualPlots(fit.best, ~  age + dis + tax + ptratio + b + lstat, 
       layout=c(2, 3), fitted=F, main="予測変数ごとの残差(b)")

#avPlots
#図4.19
avPlots(fit.best, ~ crim + zn + indus + chas + nox + rm, 
        layout=c(2, 3), main="変数追加プロット(a)", cex.lab=2)
#図4.20
avPlots(fit.best, ~  age + dis + tax + ptratio + b + lstat, 
        layout=c(2, 3), main="変数追加プロット(b)", cex.lab=2)

#influencePlot
#図4.21
influencePlot(fit.best)

#outlierTest
outlierTest(fit.best)

#influence.measures
ans <- influence.measures(fit.best)
out <- which(apply(ans$is.inf, 1, any))
length(out)

#再当てはめ
fit.4 <- update(fit.3, data=BostonHousing[-out, ])
fit.best.2 <- stepAIC(fit.4)
summary(fit.best.2)

#======================================================
#	4.3 平滑化
#======================================================
#------------<スプライン平滑化>------------

#mcycle
require(MASS)
x <- mcycle$times
y <- mcycle$accel

#多項式の当てはめ
test.poly <- lm(y ~ poly(x, 10))
pred.poly <- predict(test.poly)
#図4.22
plot(x, y)
lines(x, pred.poly)

summary(test.poly)

#スプライン回帰
library(mgcv)
kn.10 <- place.knots(x, nk=10) #節点

sm.10 <- smoothCon(s(x, bs="cr", k=10), data=data.frame(x=x),
        knots=data.frame(kn.10))[[1]]
X <- sm.10$X    #基底関数
head(X)         #最初の数行
rowSums(X)      #行合計は、すべて1

lm.cr <- lm(y ~ X - 1) #線形回帰
beta <- coef(lm.cr)    #回帰係数
pred.cr <- X %*% beta  #当てはめ値

#図4.23
plot(x, y)
lines(x, pred.cr)
abline(v = kn.10, lty=2)  #節点の位置

gam.cr <- gam(y ~ s(x, bs="cr", k=10), knots=data.frame(kn.10))
pred.gam <- predict(gam.cr)
lines(x, pred.gam, lty=3, col="blue")

#図4.24
op <- par(mfrow=c(3, 3), mar=par()$mar - 1)
for(i in 1:9) {
    plot(x, X[ ,i], type="l", lwd=3, main=paste("X", i, sep=""))
    abline(v=kn.10, lty=2)	
}
par(op)

#ペナルティ付き回帰スプライン

#ペナルティ最小
#図4.25
sp.min <- gam(y ~ s(x, bs="cr", k=80, sp=0))
sp.min.pred <- predict(sp.min)
plot(x, y)
lines(x, sp.min.pred, lwd=2)

#ペナルティ最大
#図4.26
sp.max <- gam(y ~ s(x, bs="cr", k=3), sp=100)
sp.max.pred <- predict(sp.max)
plot(x, y)
lines(x, sp.max.pred, lwd=2)

#ペナルティ行列
sm.10 <- smoothCon(s(x, bs="cr", k=10), data=data.frame(x=x),
        knots=data.frame(kn.10))[[1]]    #前述と同じ
S <- sm.10$S[[1]]    #ペナルティ行列
round(S, 2)

B <- t(mroot(S, method="svd"))

lambda <- 10    #とりあえず仮置きした値
X.aug <- rbind(X, sqrt(lambda) * B) #拡大デザイン行列
y.aug <- c(y, rep(0, nrow(B)))      #拡大応答変数 
lm.penal <- lm(y.aug ~ X.aug - 1)   #線形回帰
beta <- coef(lm.penal)    #回帰係数
pred.penal <- X %*% beta  #当てはめ値

#図4.27
plot(x, y)
lines(x, pred.penal)

#ペナルティ最適化
sp.opt <- gam(y ~ s(x, bs="cr", k=10))
sp.opt.pred <- predict(sp.opt)
lines(x, sp.opt.pred, lwd=2, lty=2)
legend("topleft", lty=1:2, lwd=1:2, 
        legend=c(expression(lambda==10), "GCV"))

#------------<局所回帰>------------

#図4.28
span <- 0.3
fit1 <- loess(y ~ x, span=span, data=mcycle)
n.num <- length(x) * span
at <- x[110]
lower <- x[110 - floor(n.num/2)]
upper <- x[110 + floor(n.num/2)]
x0 <- x[x > lower & x < upper]
y0 <- y[x > lower & x < upper]
u <- abs(x0 - at) / max(abs(x0 - at))
w <- (1 - u^3)^3
fit0 <- lm(y0 ~ x0, weights=w)
plot(x, y, type="n")
symbols(x = at+2, y = (min(y) + max(y))/2,
        rectangles=matrix(c(max(x0)-min(x0), max(y)-min(y)), ncol=2),
        inches=F, bg="gray85", fg="white", add=T)
points(x, y)
lines(0:60, predict(fit1, newdata=data.frame(x=0:60)), lty=1)
fit1.pred <- predict(fit1, newdata=data.frame(x=at))
fit0.pred <- predict(fit0, newdata=data.frame(x0=at))
diff <- fit1.pred - fit0.pred
lines(x0, predict(fit0, newdata= data.frame(x0=x0)) + diff, lty=2)
lines(x0[order(x0)], w[order(x0)] * 30 - 80, lty=3)
arrows(x0=at, y0=-130, x1=at, y1=-70, length=0.15)

#図4.29と図4.30
library(ggplot2)　　#qplotのため
qplot(x, y, data=mcycle, geom = c("point", "smooth"),
      method="loess", span=0.15, main = "loess:span = 0.15")
qplot(x, y, data=mcycle, geom = c("point", "smooth"),
      method="loess", span=0.6, main = "loess:span = 0.6")

#supsmu
#図4.31
ord <- order(x)
plot(x, y, pch=16, col="steelblue")
lines(supsmu(x[ord], y[ord]), lwd=3, col="red")
grid(col="gray20")


