# 平滑化スプラインと加法モデル

# 馬場真哉
# 2015年8月30日

# ライブラリの読み込み
# パッケージをまだインストールしていない方は、
# Rを管理者として起動したうえで、
# コメントアウトを外して、下記のコードを実行してください。
# install.packages("mgcv")

# すでにパッケージをインストールされている方は、そのまま下記のコードを実行してください
library(mgcv)


# サンプルデータ
head(airquality)

# モデルの作成
lm.model <- gam(Ozone ~ Temp, data=airquality)        #　線形回帰
gam.model <- gam(Ozone ~ s(Temp), data=airquality)    #　平滑化スプライン

# モデルの要約
summary(lm.model)
summary(gam.model)


# プロット
plot(
  gam.model,
  residuals=T, se=T, pch="。",
  main="平滑化スプラインの結果", cex.main=2
)

# 線形モデルとの比較
# 予測用のデータの作成
new <- data.frame(
  Temp=seq(min(airquality$Temp), max(airquality$Temp), 0.1)
)

# 予測
lm.pred <- predict(lm.model, new)
gam.pred <- predict(gam.model, new)

# 予測結果の図示
plot(
  airquality$Ozone ~ airquality$Temp,
  xlab="Temp", ylab="Ozone", 
  main="推定結果", cex.main=2, cex.lab=1.5
)
lines(lm.pred ~ as.matrix(new), col=2, lwd=2)
lines(gam.pred ~ as.matrix(new), col=4, lwd=2)
legend("topleft", lwd=2, col=c(2, 4), legend=c("線形 回帰", "平滑化スプライン"))


# 検定
anova(lm.model, gam.model, test="F")

# 関数を使わないで検定
names(summary(gam.model))

DF <- summary(gam.model)$n
gam.dev <- sum(resid(gam.model)^2)
lm.dev <- sum(resid(lm.model)^2)
gam.df <- DF-summary(gam.model)$residual.df
lm.df <- DF-summary(lm.model)$residual.df
F <- ((lm.dev-gam.dev)/(gam.df-lm.df)) / (gam.dev/(DF-gam.df))

F
1-pf(F, gam.df-lm.df, DF-gam.df)


# グネグネ度（平滑化パラメータ）の確認
# spを0.001から0.1まで、0.001刻みで変化させる
sp <- seq(from=0.001, to=0.1, by=0.001)

# spごとに、GCVを計算
GCV <- numeric()
for(i in 1:length(sp)){
  g.m <- gam(Ozone ~ s(Temp), sp=sp[i], data=airquality)
  GCV[i] <- g.m$gcv.ubre
}

# spとGCVの関係をプロット
plot(GCV ~ sp, main="GCVとグネグネ度", cex.main=2)

# mgcv関数で選ばれたspを重ねてプロットする
points(gam.model$sp, gam.model$gcv.ubre, col=2, pch=18, cex=2)


# 平滑化パラメータを変更したときの、推定結果の変化の図示
par(mfrow=c(2, 2))
plot(gam(Ozone ~ s(Temp), sp=1000, data=airquality), residuals=T)
plot(gam(Ozone ~ s(Temp), sp=1, data=airquality), residuals=T)
plot(gam(Ozone ~ s(Temp), sp=0.00001, data=airquality), residuals=T)
plot(gam(Ozone ~ s(Temp, k=20), sp=0.0000001, data=airquality), residuals=T)
par(mfrow=c(1, 1))


# むりやり直線にした平滑化スプラインと線形回帰の比較
lm.modoki <- gam(Ozone ~ s(Temp), sp=100000, data=airquality)
plot(
  lm.modoki$fitted.values ~ lm.model$fitted.values,
  xlab="線形回帰の予測結", ylab="グネグネ度最小スプラインの予測結果", cex.lab=1.2
)
abline(a=0, b=1, col=2)

# モデルの評価
gam.check(gam.model)

# 薄板平滑化スプライン
gam2 <- gam(Ozone ~ s(Wind, Temp), data=airquality)
# 3Dで図示
vis.gam(gam2, color="cm", theta=45)

# 等高線
vis.gam(gam2, color="cm", plot.type="contour")

# たくさん変数がある場合の図示
gg <- gam(Ozone ~ s(Solar.R) + s(Wind, Temp), data=airquality)
vis.gam(gg, view=c("Wind", "Temp"), color="cm", theta=45)


# 予測精度の評価
# テストデータ作成
# 欠損値の除去
Air <- na.omit(airquality)
length(Air[, 1])

# データを訓練データとテストデータに分ける
set.seed(1)
S1 <- sample(1:111, size=30)
teach <- Air[-S1, ]
test <- Air[S1, ]

# 予測モデルの作成
l.model <- lm(Ozone ~ Solar.R + Wind + Temp, data=teach)
g.model <- gam(Ozone ~ s(Solar.R) + s(Wind) + s(Temp), data=teach)

# AICによるモデル選択
library(MuMIn)
options(na.action = "na.fail")

# 線形モデルでの変数選択
a1 <- dredge(l.model, rank="AIC")
a1

# 加法モデルでの変数選択
a2 <- dredge(g.model, rank="AIC")
a2

# 予測
lm.pred <- predict(l.model, test)
gam.pred <- predict(g.model, test)

# 予測誤差（残差平方和）の計算
lm.resid <- sum((lm.pred-test$Ozone)^2)
gam.resid <- sum((gam.pred-test$Ozone)^2)

# 予測結果の評価
lm.resid
gam.resid
lm.resid - gam.resid