本章では、説明変数と目的変数の関係が線形とは限らない状況に対処するための手法を扱う。多項式回帰やステップ関数、スプライン、一般化加法モデル（Generalized Additive Models; GAM）といった非線形なモデル化手法により、データに内在する曲線的なパターンを捉えることができる。これらの手法によって、線形モデルを拡張し柔軟性を高めることで予測精度の向上や関係性の解釈が可能となる。本節ではまず、それら手法の実装例を示し、続いて練習問題を通じて実践的に分析と考察を行う。

目次

  1. R 文法ウォームアップ
  2. 7.8.1 多項式回帰とステップ関数
  3. 7.8.2 スプライン
  4. 7.8.3 一般化加法モデル
  5. 7.9 練習問題

0 R 文法ウォームアップ（7章で使う関数集中講座）

このセクションで学ぶこと

  • Chapter 7 で多用する poly(), cut(), bs(), ns(), loess(), gam() など 非線形回帰用ファンクションの最小例 を実行して振る舞いを確認。
  • 設計行列 (model.matrix())、係数表 (coef() / summary())、および 予測専用データフレーム という 3 つの基本ツールを今一度おさらい。
  • 各チャンクは「▶ 入力 → ★ 標準出力」の対を成し、動かした直後にポイントを確認できるよう設計されている。

読み進め方の指針

  1. コードをコピーして RStudio Console か Quarto で実行。
  2. すぐ下の出力や ポイント を読み、「何が起きたか」を必ず ひと言で 言語化。
  3. 値を書き換えて再実行し、関数の引数が曲線の形をどう変えるか を体感的に把握。

0‑1 poly() と基底関数

x <- 1:5
X <- model.matrix(~ poly(x, 3))
X                       # 直交基底（列は次数 0,1,2,3）
##   (Intercept)   poly(x, 3)1 poly(x, 3)2   poly(x, 3)3
## 1           1 -6.324555e-01   0.5345225 -3.162278e-01
## 2           1 -3.162278e-01  -0.2672612  6.324555e-01
## 3           1 -3.510833e-17  -0.5345225  1.755417e-16
## 4           1  3.162278e-01  -0.2672612 -6.324555e-01
## 5           1  6.324555e-01   0.5345225  3.162278e-01
## attr(,"assign")
## [1] 0 1 1 1

ポイント

  • poly(x, 3) は 次数 3 の直交多項式基底 を返す。
  • model.matrix(~ poly(x, 3)) は、切片を含む設計行列を生成。
  • 直交化 されているので、列間相関が 0 — 推定が数値的に安定。

0‑2 cut() とステップ関数

ages <- c(22, 35, 47, 55, 68)
cut(ages, 4)            # 4 区間に分割
## [1] (22,33.5] (33.5,45] (45,56.5] (45,56.5] (56.5,68]
## Levels: (22,33.5] (33.5,45] (45,56.5] (56.5,68]

ポイント

  • cut() は ビン分割 を行い、カテゴリ型 (factor) を返す。
  • 返り値はそのまま回帰式に投入でき、区間ごと定数 のステップ関数モデルになる。

0‑3 bs() / ns() とスプライン

library(splines)
bs(ages, knots = c(40, 60))
##                1         2          3         4 5
## [1,] 0.000000000 0.0000000 0.00000000 0.0000000 0
## [2,] 0.500308262 0.4084324 0.06982583 0.0000000 0
## [3,] 0.076073407 0.4613551 0.44069651 0.0218750 0
## [4,] 0.004328255 0.1883307 0.59209870 0.2152423 0
## [5,] 0.000000000 0.0000000 0.00000000 0.0000000 1
## attr(,"degree")
## [1] 3
## attr(,"knots")
## [1] 40 60
## attr(,"Boundary.knots")
## [1] 22 68
## attr(,"intercept")
## [1] FALSE
## attr(,"class")
## [1] "bs"     "basis"  "matrix"
ns(ages, df = 4)
##              1          2         3           4
## [1,] 0.0000000  0.0000000 0.0000000  0.00000000
## [2,] 0.2048485 -0.1883701 0.4926602 -0.30429014
## [3,] 0.6848485  0.1689201 0.1288383 -0.07957663
## [4,] 0.2438672  0.5888590 0.2030530 -0.03577925
## [5,] 0.0000000 -0.1545866 0.4043035  0.75028313
## attr(,"degree")
## [1] 3
## attr(,"knots")
## [1] 35 47 55
## attr(,"Boundary.knots")
## [1] 22 68
## attr(,"intercept")
## [1] FALSE
## attr(,"class")
## [1] "ns"     "basis"  "matrix"

ポイント

  • bs() は指定した knots で区間を分けた B スプライン基底 を返す。
  • ns() は端点で線形外挿する 自然スプライン。外れ値に強い。

0‑4 loess() と局所回帰

df <- data.frame(x = 1:10, y = (1:10) + rnorm(10))
fit_lo <- loess(y ~ x, span = 0.5, data = df)
predict(fit_lo, data.frame(x = 5))
##        1 
## 4.523655

ポイント

  • span は 局所フィットに使うデータ割合。小さいほど曲線はうねる。
  • predict() には 同じ列名を持つデータフレーム を与えるのが鉄則。

0‑5 gam() と加法モデル

library(gam)
## Loading required package: foreach
## Loaded gam 1.22-4
df2 <- data.frame(x = runif(100), z = runif(100))
df2$y <- sin(2 * pi * df2$x) + df2$z + rnorm(100, sd = 0.2)
fit_gam <- gam(y ~ s(x, 5) + z, data = df2)
summary(fit_gam)
## 
## Call: gam(formula = y ~ s(x, 5) + z, data = df2)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47750 -0.13115  0.01261  0.11171  0.41214 
## 
## (Dispersion Parameter for gaussian family taken to be 0.0372)
## 
##     Null Deviance: 58.0299 on 99 degrees of freedom
## Residual Deviance: 3.4602 on 93.0002 degrees of freedom
## AIC: -36.598 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## s(x, 5)    1 34.635  34.635  930.91 < 2.2e-16 ***
## z          1  9.030   9.030  242.70 < 2.2e-16 ***
## Residuals 93  3.460   0.037                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar F     Pr(F)    
## (Intercept)                             
## s(x, 5)           4 101.53 < 2.2e-16 ***
## z                                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ポイント

  • s(x, 5) は 自由度 5 の平滑関数。裏側では B スプラインにペナルティを課している。
  • gam() 1 行で「非線形 + 線形」のハイブリッドモデルが書ける。

7.8 実習：非線形モデル

本実習では、パッケージ ISLR に含まれる Wage データセットを用いて、線形性に囚われない回帰手法を実践する。Wage データは大西洋中部地域の男性労働者3000人の賃金と職歴・学歴などの情報を含む。以下では年齢による賃金の予測を題材に、多項式回帰、ステップ関数、スプライン、局所回帰、そして一般化加法モデルの順に分析を行い、それぞれの手法による当てはまりの違いや結果の解釈を示す。

このセクションで何をするか

  • poly() による 多項式回帰 と cut() による ステップ関数 を比較し、
    「滑らかさ vs. 解釈容易性」という観点でモデルを読み解く。
  • ANOVA と 交差検証 (CV) を組み合わせ、最適次数（または区間数）を選択。
  • 推定後は 当てはまり曲線の可視化 と 係数の有意性 を丁寧に読み解く。

先にゴールを把握 → コードを実行 → 出力と図を読み解く、の順で進もう。

7.8.1 多項式回帰とステップ関数

まず、年齢 (age) から賃金 (wage) を予測する多項式回帰モデルを当てはめる。次数4の多項式で回帰を行い、係数の推定結果を確認する。なお、poly(age, 4) は直交多項式による基底を利用しており、推定の数値安定性が高い（raw引数を指定しなければデフォルトで直交化）。

library(ISLR)
attach(Wage)
# 年齢に対する4次多項式回帰（直交多項式基底）
fit_poly4 <- lm(wage ~ poly(age, 4), data = Wage)
coef(summary(fit_poly4))
##                 Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)    111.70361  0.7287409 153.283015 0.000000e+00
## poly(age, 4)1  447.06785 39.9147851  11.200558 1.484604e-28
## poly(age, 4)2 -478.31581 39.9147851 -11.983424 2.355831e-32
## poly(age, 4)3  125.52169 39.9147851   3.144742 1.678622e-03
## poly(age, 4)4  -77.91118 39.9147851  -1.951938 5.103865e-02

上記の結果、年齢の4次までの項がいずれも統計的に有意であることが確認できる（第1〜第3次は p<0.01、第4次は p≒0.05）。これは賃金と年齢の関係に高次の曲線効果が存在することを示唆している。同じ回帰をraw=TRUEとして非直交（生の）多項式で行っても、予測の当てはまり自体は変わらず、得られる係数は多項式展開した場合の値になる。例えば次数4の場合、非直交多項式回帰の推定係数は次のようになる。

# 生の多項式での4次回帰（結果の係数は多項式展開後の値）
fit_poly4_raw <- lm(wage ~ poly(age, 4, raw = TRUE), data = Wage)
coef(summary(fit_poly4_raw))
##                                Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)               -1.841542e+02 6.004038e+01 -3.067172 0.0021802539
## poly(age, 4, raw = TRUE)1  2.124552e+01 5.886748e+00  3.609042 0.0003123618
## poly(age, 4, raw = TRUE)2 -5.638593e-01 2.061083e-01 -2.735743 0.0062606446
## poly(age, 4, raw = TRUE)3  6.810688e-03 3.065931e-03  2.221409 0.0263977518
## poly(age, 4, raw = TRUE)4 -3.203830e-05 1.641359e-05 -1.951938 0.0510386498

直交基底であれ生の多項式であれ、モデルの予測値自体は同一であることを確かめておく。実際、fit_poly4 と fit_poly4_raw に対して同じ入力を与えたときの予測値の差はゼロとなる。これは多項式基底の取り方によらず、最終的な当てはまり曲線は一致するためである。

# 直交vs非直交、多項式回帰モデルの予測値差分の最大値
preds1 <- predict(fit_poly4, Wage)
preds2 <- predict(fit_poly4_raw, Wage)
max(abs(preds1 - preds2))
## [1] 6.842527e-11

次に、年齢を説明変数に用いた多項式回帰モデルについて、最適な多項式の次数を検討する。仮説検定による比較のため、一次から五次までのモデルをあてはめ、ANOVAによる逐次F検定を行う。

fit_1 <- lm(wage ~ age, data = Wage)                   # 1次
fit_2 <- lm(wage ~ poly(age, 2), data = Wage)          # 2次
fit_3 <- lm(wage ~ poly(age, 3), data = Wage)          # 3次
fit_4 <- lm(wage ~ poly(age, 4), data = Wage)          # 4次
fit_5 <- lm(wage ~ poly(age, 5), data = Wage)          # 5次
anova(fit_1, fit_2, fit_3, fit_4, fit_5)
## Analysis of Variance Table
## 
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
##   Res.Df     RSS Df Sum of Sq        F    Pr(>F)    
## 1   2998 5022216                                    
## 2   2997 4793430  1    228786 143.5931 < 2.2e-16 ***
## 3   2996 4777674  1     15756   9.8888  0.001679 ** 
## 4   2995 4771604  1      6070   3.8098  0.051046 .  
## 5   2994 4770322  1      1283   0.8050  0.369682    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

上の分散分析の結果、年齢に対する賃金の関係は4次の多項式まで有意に改善することが確認できる。例えば一次から二次への拡張でRSSが大きく減少し (p < 2.2e-16)、三次も有意に改善 (p = 0.00168) 、四次もわずかながら改善が認められるものの (p = 0.051) 、五次以降は有意な改善が見られない (p ≈ 0.37)。従って、このデータでは4次多項式程度で十分であり、過度に高次の項を入れる必要はないと判断される。実際、4次と5次のモデルの予測曲線はほとんど差がない。

次に、年齢をカテゴリカル変数に変換するステップ関数モデルを試す。ここでは cut() 関数を用いて年齢を4つの区間に分割し、それをダミー変数化して線形モデルに組み込む。4分割（5歳区切り程度）のステップ関数で回帰した場合の推定結果は以下の通りである。

# 年齢を4区間に分割したステップ関数による回帰
fit_step4 <- lm(wage ~ cut(age, 4), data = Wage)
coef(summary(fit_step4))
##                         Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)            94.158392   1.476069 63.789970 0.000000e+00
## cut(age, 4)(33.5,49]   24.053491   1.829431 13.148074 1.982315e-38
## cut(age, 4)(49,64.5]   23.664559   2.067958 11.443444 1.040750e-29
## cut(age, 4)(64.5,80.1]  7.640592   4.987424  1.531972 1.256350e-01

出力された係数を見ると、基準年齢層（最も若い層）に比べ、中年層・熟年層で賃金が有意に高いことが分かる。例えば、33.5～49歳および49～64.5歳のカテゴリの係数はそれぞれ約24と23.7であり、両カテゴリともp値は極めて有意である。一方、最高齢層（64.5～80.1歳）の係数7.64は有意とは言えず、この層の賃金分布は基準層と大きく変わらない可能性が示唆される。

以上の多項式回帰とステップ関数モデルの当てはまりを可視化して比較する。図7.1の左は4次多項式回帰、右は4区間のステップ関数による予測曲線である。灰色の点は元データで、青線がモデルの予測値、点線は予測の標準誤差による95%信頼帯を示す。

図 7.1 年齢と賃金：4 次多項式 (青) と 4 分割ステップ関数 (赤)

図 7.1 年齢と賃金：4 次多項式 (青) と 4 分割ステップ関数 (赤)

図7.1 年齢と賃金：4次多項式回帰（左）と4分割ステップ関数モデル（右）の当てはまり

4次多項式曲線（左青線）は滑らかに賃金の年齢効果を捉えている。一方、ステップ関数モデル（右赤線）は年齢を4つの区間に区切って水準の異なる定数で近似しており、各区間内では一定、区切りで不連続的に変化する形で賃金の推移を表現している。多項式モデルは滑らかだが僅かなオーバーフィッティングの兆候（高齢層での不安定な振る舞い）が見られるのに対し、ステップ関数モデルは高齢層の水準を過小評価しているように見える。モデル選択の観点では、前述のANOVA結果より4次程度の多項式が妥当であった。交差検証による評価でも、テスト誤差（MSE）は多項式4次でほぼ最小となり、それ以上次数を上げても大きな改善はない。一方ステップ関数の場合、8分割程度までテスト誤差が改善し、8区間モデルが最良と推定された。一般に、データに十分対応できる複雑さ（多項式の次数や区間の数）を持つモデルを選ぶには、このように交差検証を用いると客観的な指標で判断できる。

このセクションで何をするか

  • bs() と ns() を使って B スプライン と 自然スプライン を学ぶ。
  • スプラインが持つ 局所性・滑らかさ・端点挙動 の 3 点に注目し、
    ノット配置と自由度が曲線をどう変えるかを体感する。
  • 平滑化スプライン (smooth.spline()) との比較で「過学習 vs. バイアス」を考察。

7.8.2 スプライン

スプライン回帰は、多項式を区間ごとに繋ぎ合わせて滑らかな曲線を描く手法である。Rではsplinesパッケージを用いることでBスプライン基底を簡便に利用できる。まず、年齢に3つのノット（25歳, 40歳, 60歳）を置いた三次Bスプラインで賃金を回帰し、結果の当てはまり曲線を描いてみる。同時に、自然スプライン（Natural Spline）による4自由度のモデルも当てはめ、曲線を比較する。

library(splines)
# 3つのノットを指定した三次Bスプライン
fit_bs <- lm(wage ~ bs(age, knots = c(25, 40, 60)), data = Wage)
# 4自由度の自然スプライン
fit_ns <- lm(wage ~ ns(age, df = 4), data = Wage)

# 予測用の年齢グリッド
age.range <- range(Wage$age)
age.grid <- seq(from = age.range[1], to = age.range[2])
# 各モデルの予測値
pred_bs <- predict(fit_bs, newdata = list(age = age.grid), se = TRUE)
pred_ns <- predict(fit_ns, newdata = list(age = age.grid), se = TRUE)

図7.2に年齢と賃金の散布図と、Bスプラインモデル（黒実線）および自然スプラインモデル（赤実線）による推定曲線を示す（点線は各モデルの±2標準誤差帯）。

図7.2 年齢と賃金：Bスプライン（黒）と自然スプライン（赤）による当てはまり曲線

両モデルともデータの傾向をなめらかに表現しており、中年以降に賃金が頭打ちになる非線形なパターンを捉えている。Bスプラインモデルではノット位置（25, 40, 60歳）で多項式片を繋ぎ合わせているが、三次スプラインのため2次までの導関数が連続になるよう滑らかに接続されている。自然スプラインはデータ範囲の両端で曲線が直線（2次導関数が0）になる制約を課したスプラインで、外挿性を改善する効果がある。図では、自然スプラインの曲線（赤）はBスプライン（黒）とほぼ重なっており、データ範囲内での当てはまりに大差はない。適切にノット数や位置を選べば、多項式より柔軟かつ過学習しにくい曲線モデルが得られる。

次に平滑化スプライン（smoothing spline）を適用する。平滑化スプラインはノットを全データ点に置きつつ適度な平滑化を行う手法で、自由度（有効パラメータ数）によって曲線の滑らかさを制御する。例えば自由度16で当てはめた場合と、交差検証で選択した自由度で当てはめた場合とを比較してみる。交差検証により選ばれたのは約6.8自由度のモデルである（出力の df＝6.7946 は有効自由度）。

# 自由度16の平滑化スプラインと、CVにより自由度自動選択した平滑化スプライン
fit_ss16 <- smooth.spline(Wage$age, Wage$wage, df = 16)
fit_ssCV <- smooth.spline(Wage$age, Wage$wage, cv = TRUE)
## Warning in smooth.spline(Wage$age, Wage$wage, cv = TRUE): cross-validation with
## non-unique 'x' values seems doubtful
fit_ssCV$df  # 選ばれた有効自由度
## [1] 6.794596
図 7.2 年齢と賃金：B スプライン (黒) と自然スプライン (赤)

図 7.2 年齢と賃金：B スプライン (黒) と自然スプライン (赤)

図7.3に、平滑化スプライン（赤：16DF、青：CV選択≈6.8DF）の当てはまり曲線を示す。凡例に示す通り、自由度16のモデルはより屈曲が多く、自由度6.8のモデルは滑らかである。

## [1] 6.794596
図 7.3 年齢と賃金：平滑化スプライン（赤＝16 df，青＝CV 選択 df）

図 7.3 年齢と賃金：平滑化スプライン（赤＝16 df，青＝CV 選択 df）

図7.3 年齢と賃金：平滑化スプライン（赤：16自由度、青：CV選択自由度）の当てはまり比較

交差検証で選ばれた約6.8自由度の平滑化スプラインは、データから自動的に適切な平滑度を見積もった結果である。自由度16のモデルでは細かな変動まで追従しており訓練データへの適合は良いが、明らかに過剰適合気味である。それに対し自由度6.8のモデルは大局的な趨勢を捉えており、新たなデータに対しても良好な予測性能を示すことが期待できる。平滑化スプラインは、このように自動的に複雑さを調整できる点が大きな利点である。

最後に局所回帰（局所多項式回帰）の例を示す。loess() 関数で年齢から賃金への局所回帰曲線を当てはめ、異なるスパン（スムージング範囲）での結果を比較する。ここではスパン0.2と0.5の2種類のモデルを用意し、当てはまり曲線を描画する。

# 局所回帰（スパン0.2と0.5）
fit_loess_1 <- loess(wage ~ age, span = 0.2, data = Wage)
fit_loess_2 <- loess(wage ~ age, span = 0.5, data = Wage)
# グリッド上で予測
pred_lo1 <- predict(fit_loess_1, data.frame(age = age.grid))
pred_lo2 <- predict(fit_loess_2, data.frame(age = age.grid))

図7.4に、局所回帰曲線（赤：スパン0.2、青：スパン0.5）をプロットする。スパンとは局所回帰で用いるデータ範囲の割合で、値が小さいほど局所的（狭い範囲）な当てはめになる。

図 7.4 年齢と賃金：局所回帰 (赤 = 0.2, 青 = 0.5)

図 7.4 年齢と賃金：局所回帰 (赤 = 0.2, 青 = 0.5)

図7.4 年齢と賃金：局所回帰（赤：スパン0.2、青：スパン0.5）の当てはまり比較

スパン0.2（赤線）のモデルはデータの細部まで追随し曲線がより屈曲している。一方スパン0.5（青線）ではより広い範囲のデータを平均化するため滑らかな曲線になっている。スパンの値はバンド幅に相当し、モデルのバイアス・分散トレードオフを調整するパラメータである。一般にスパンを大きくすると曲線はなめらかになりバイアスが増えるが分散は減少し、小さくすると曲線は変動が大きくなるが分散は増える。適切なスパンはデータに依存するが、クロスバリデーション等で選択可能である。

このセクションで何をするか

  • gam() を用いて 一般化加法モデル (GAM) を構築し、
    「非線形項 + カテゴリ変数 + 線形項」を ひとつの数式 で表すテクニックを習得。
  • 部分残差プロットで 各説明変数の効果を視覚化 し、
    モデルの当てはまりを直感的に評価。
  • GAM と通常の線形モデルを比較し、柔軟性と解釈性のバランス を議論。

7.8.3 一般化加法モデル

一般化加法モデル (GAM) は、各説明変数に対して任意の非線形関数を用いて加法的に効果を表現するモデルである。本節では賃金データにGAMを適用し、線形モデルと比較する。まず、年齢・年（調査年）・学歴が賃金に与える効果を自然スプラインでモデル化し、通常の線形回帰として当てはめる。

# 年（year）と年齢（age）は自然スプライン、学歴（education）はカテゴリ変数で線形回帰
gam1 <- lm(wage ~ ns(year, 4) + ns(age, 5) + education, data = Wage)
summary(gam1)
## 
## Call:
## lm(formula = wage ~ ns(year, 4) + ns(age, 5) + education, data = Wage)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -120.513  -19.608   -3.583   14.112  214.535 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   46.949      4.704   9.980  < 2e-16 ***
## ns(year, 4)1                   8.625      3.466   2.488  0.01289 *  
## ns(year, 4)2                   3.762      2.959   1.271  0.20369    
## ns(year, 4)3                   8.127      4.211   1.930  0.05375 .  
## ns(year, 4)4                   6.806      2.397   2.840  0.00455 ** 
## ns(age, 5)1                   45.170      4.193  10.771  < 2e-16 ***
## ns(age, 5)2                   38.450      5.076   7.575 4.78e-14 ***
## ns(age, 5)3                   34.239      4.383   7.813 7.69e-15 ***
## ns(age, 5)4                   48.678     10.572   4.605 4.31e-06 ***
## ns(age, 5)5                    6.557      8.367   0.784  0.43328    
## education2. HS Grad           10.983      2.430   4.520 6.43e-06 ***
## education3. Some College      23.473      2.562   9.163  < 2e-16 ***
## education4. College Grad      38.314      2.547  15.042  < 2e-16 ***
## education5. Advanced Degree   62.554      2.761  22.654  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.16 on 2986 degrees of freedom
## Multiple R-squared:  0.293,  Adjusted R-squared:  0.2899 
## F-statistic:  95.2 on 13 and 2986 DF,  p-value: < 2.2e-16

上記の gam1 は実質的にGAMと同等のモデル（年と年齢について非線形項を含む）である。回帰の結果によれば、年齢および年の非線形効果は有意であり、特に年齢の効果が強く非線形性を持つことが示唆される（ns(age,5)の各基底関数に対するt値が大きく有意）。学歴（5水準の質的変数）も高い有意性を示しており、賃金に強い影響を与えることがわかる。

次に、gam パッケージを用いて同じモデルを適用し、部分残差プロットで効果を視覚化する。gam() 関数では s( , df) や lo() を用いて変数ごとのスムーズ関数を指定できる。まず gam1 と同じ構造（年を4自由度スプライン、年齢を5自由度スプライン）のモデルを gam() で推定し、プロットしてみる。

library(gam)
gam_m3 <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
par(mfrow = c(1,3))
plot(gam_m3, se = TRUE, col = 'blue')

上図（図7.5）は各説明変数に対する部分効果を描いたもので、点線は標準誤差による信頼帯を示す。年（左図）については、おおむね直線的な上昇傾向が見られるものの、わずかに非線形な揺らぎも描かれている。年齢（中央図）は40歳前後で効果が頭打ちになる顕著な曲線関係が確認できる。学歴（右図）は質的変数であり、高学歴になるほど賃金が上昇する明確な効果が示されている。

次に、変数選択や効果の線形性の検定を行うため、いくつかモデルを比較する。年の効果が本当に非線形である必要があるか検討するため、以下の3つのモデルを比較する。

  • モデル1: 年を除き、年齢を非線形 (5自由度スプライン) + 学歴
  • モデル2: 年を線形項として加え、年齢を非線形 + 学歴
  • モデル3: 年を非線形 (4自由度スプライン) として加え、年齢を非線形 + 学歴
gam_m1 <- gam(wage ~ s(age, 5) + education, data = Wage)
gam_m2 <- gam(wage ~ year + s(age, 5) + education, data = Wage)
gam_m3 <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
anova(gam_m1, gam_m2, gam_m3, test = "F")
## Analysis of Deviance Table
## 
## Model 1: wage ~ s(age, 5) + education
## Model 2: wage ~ year + s(age, 5) + education
## Model 3: wage ~ s(year, 4) + s(age, 5) + education
##   Resid. Df Resid. Dev Df Deviance       F    Pr(>F)    
## 1      2990    3711731                                  
## 2      2989    3693842  1  17889.2 14.4771 0.0001447 ***
## 3      2986    3689770  3   4071.1  1.0982 0.3485661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

上の 分散分析テーブルによれば、モデル1 vs モデル2: 年を説明変数に加えると有意に残差偏差が減少しており (p = 1.44e-4)、年には有意な効果がある。一方、モデル2 vs モデル3: 年を線形から非線形に変更しても残差偏差の減少は有意でない (p = 0.3485)。すなわち、年の効果は線形項で十分であり、非線形にする必要は認められないことが分かる。以上から、最終的に年は線形効果、年齢は非線形効果、学歴は質的効果というモデル（モデル2に相当）が最も適切と判断される。

最後に、最終モデル (gam_m2) の推定結果を確認し、どの変数に非線形性が必要かを検証する。summary(gam_m3) を実行すると、各変数の汎任意効果についてANOVAによる検定が表示される。年については「Parametric Effects」において有意（p < 0.001）であり、「Nonparametric Effects」における自由度3の効果はF=1.086, p=0.354と非線形効果が有意でない。一方、年齢については「Nonparametric Effects」でF=32.38, p<2e-16と強い非線形効果の有意性が確認できる。学歴は質的変数のため非線形効果の概念は無いが、Parametric Effectsで非常に有意である。以上より、年齢と賃金の関係は明確に非線形だが、年と賃金の関係はほぼ線形であることが統計的に裏付けられた。

7.9 練習問題

以下では、第7章で学んだ非線形モデル手法を実際のデータに適用し、結果を考察する。

問題6. 本章で用いた Wage データセットをさらに分析する。（※データの説明は前述）

(a) 年齢を用いて賃金を予測する多項式回帰を行い、交差検証によって多項式の最適次数 $d$ を選択せよ。最適と選ばれた次数はいくつか。また、それは前述のANOVAによる仮説検定の結果と一致しているか比較せよ。最後に、その最適モデルによる当てはまり曲線をデータとともにプロットせよ。

library(boot)
set.seed(1)
# 1〜10次までの多項式についてCVでMSEを算出
cv.error <- rep(NA, 10)
for(i in 1:10){
  fit_i <- glm(wage ~ poly(age, i), data = Wage)
  cv.error[i] <- cv.glm(Wage, fit_i, K = 10)$delta[1]  # K=10のCV推定誤差
}
plot(1:10, cv.error, type="l", xlab="Degree", ylab="Test MSE")
deg.min <- which.min(cv.error)
points(deg.min, cv.error[deg.min], col="red", pch=19, cex=1.5)

deg.min
## [1] 9

上記の交差検証の結果、テスト誤差が最小となる次数は 9 であった（赤点で強調）。しかし、得られた曲線を見ると4次から先は誤差減少幅がごく小さい。実際、4次と9次のテストMSEは僅差であり、モデルの複雑さに見合った改善とは言い難い。これは前問のANOVAによる検定結果とも概ね一致している。ANOVAでは4次から5次へのモデル複雑化は有意でなく、4次で十分と判断された。したがって、実用上は4次程度の多項式で十分であり、9次など高次のモデルは過剰適合のリスクが高い。図7.5にCV最小となる9次多項式モデルの当てはまり曲線（赤線）を示すが、4次モデルと視覚的に大差なく、データのばらつきを局所的に追い過ぎている印象を受ける。

図 7.5 Wageデータ：年齢に対する 9 次多項式回帰の当てはまり（赤線）

図 7.5 Wageデータ：年齢に対する 9 次多項式回帰の当てはまり（赤線）

図7.5 Wageデータ：年齢に対する9次多項式回帰の当てはまり（赤線）

(b) 年齢を用いて賃金を予測するステップ関数モデルを構築し、交差検証によって最適な区間数（カット数）を選択せよ。最適と判断された区間数はいくつか。そして、そのモデルによる当てはまり曲線をデータとともにプロットせよ。

cv.error <- rep(NA, 10)
for(i in 2:10){
  Wage$age_cut <- cut(Wage$age, i)
  fit_i <- glm(wage ~ age_cut, data = Wage)
  cv.error[i] <- cv.glm(Wage, fit_i, K = 10)$delta[1]
}
plot(2:10, cv.error[-1], type="l", xlab="Cuts", ylab="Test MSE")
cut.min <- which.min(cv.error)
points(cut.min, cv.error[cut.min], col="red", pch=19, cex=1.5)

cut.min
## [1] 8

交差検証の結果、テスト誤差が最小となる年齢区分の数は 8区間 であった。すなわち、年齢を8つの層に分けたステップ関数モデルが最も予測精度が高いと判断された。図7.6に8区間モデルの当てはまり（赤線）を示す。8つの年齢層ごとに賃金の水準が異なる形でフィットしており、各区間内では一定値となるステップ状の関数となっている。細かく区間分割したことで高齢層の賃金低下もある程度表現できているが、区間境界で不連続な跳びが発生する点に注意が必要である。

図 7.6 Wageデータ：年齢 8 区間のステップ関数モデルによる当てはまり（赤線）

図 7.6 Wageデータ：年齢 8 区間のステップ関数モデルによる当てはまり（赤線）

図7.6 Wageデータ：年齢8区間のステップ関数モデルによる当てはまり（赤線）

(c) 質的説明変数の効果を調べる。Wageデータには未使用の変数として婚姻状況（maritl）や職種（jobclass）などがある。これらと賃金の関係を探索し、非線形手法を用いて柔軟なモデルをあてはめよ。結果のプロットを作成し、考察を述べよ。

まず、婚姻状況と職種の分布を確認する。

summary(Wage$maritl)
## 1. Never Married       2. Married       3. Widowed      4. Divorced 
##              648             2074               19              204 
##     5. Separated 
##               55
summary(Wage$jobclass)
##  1. Industrial 2. Information 
##           1544           1456
par(mfrow = c(1,2))
plot(Wage$maritl, Wage$wage, xlab="Marital Status", ylab="Wage")
plot(Wage$jobclass, Wage$wage, xlab="Job Class", ylab="Wage")

婚姻状況は「Never Married（未婚）」が648人、「Married（既婚）」が2074人と大半を占め、他はごく少数である。職種は「Industrial（工業系）」1544人、「Information（情報系）」1456人の2水準で、ほぼ半々に分かれている。賃金の分布を見ると（上のボックスプロット）、既婚者の賃金中央値は未婚者より高く、離婚者等はばらつきが大きい傾向がある。また情報系職種の賃金は工業系より若干高めに見える。こうした質的変数の効果も組み込むため、以下では GAM を用いて分析する。

年（year）と年齢（age）はそれぞれ非線形項（局所回帰 lo() およびスプライン s()）で、学歴・婚姻・職種は質的変数として、徐々にモデルに加えていく。まず婚姻と職種を除いたベースモデル（年・年齢・学歴）を構築し、順次それらを追加してモデル間で比較する。

library(gam)
# ベースモデル: yearをlo(span=0.7)、ageをスプライン5DF、educationをカテゴリ
fit1 <- gam(wage ~ lo(year, span=0.7) + s(age, 5) + education, data = Wage)
# 職種(jobclass)を追加
fit2 <- gam(wage ~ lo(year, span=0.7) + s(age, 5) + education + jobclass, data = Wage)
# 婚姻(maritl)を追加
fit3 <- gam(wage ~ lo(year, span=0.7) + s(age, 5) + education + maritl, data = Wage)
# 両方追加
fit4 <- gam(wage ~ lo(year, span=0.7) + s(age, 5) + education + jobclass + maritl, data = Wage)
anova(fit1, fit2, fit3, fit4, test="F")
## Analysis of Deviance Table
## 
## Model 1: wage ~ lo(year, span = 0.7) + s(age, 5) + education
## Model 2: wage ~ lo(year, span = 0.7) + s(age, 5) + education + jobclass
## Model 3: wage ~ lo(year, span = 0.7) + s(age, 5) + education + maritl
## Model 4: wage ~ lo(year, span = 0.7) + s(age, 5) + education + jobclass + 
##     maritl
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1    2987.1    3691855                                 
## 2    2986.1    3679689  1    12166 10.124 0.0014788 ** 
## 3    2983.1    3597526  3    82163 22.790 1.386e-14 ***
## 4    2982.1    3583675  1    13852 11.526 0.0006952 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

上の分析対象偏差の表によれば、職種の追加 (fit1→fit2) で残差偏差が有意に減少 (p = 0.00146)、婚姻の追加 (fit1→fit3) でも大きく減少している (p ≈ 9.5e-15)。さらに両方を含むモデル (fit4) は、それぞれ単独追加モデルと比べても有意に改善している (fit3→fit4: p = 6.86e-04)。したがって、職種と婚姻はいずれも賃金に有意な影響を及ぼし、両方含めたモデルが最も適合が良い。最終モデル fit4 の部分プロットを図7.7に示す。

par(mfrow = c(2,3))
plot(fit4, se = TRUE, col = "blue")

図7.7では、年・年齢・学歴に加え、職種（下段左）と婚姻（下段中央）の効果が描かれている。既婚（Married）は未婚（Never Married）に比べ明らかにプラスの効果を持ち、賃金水準を押し上げていることが青線より読み取れる。また情報系職種（Information）は工業系（Industrial）に比べ賃金が高い傾向が示されている。年齢効果は以前と同様に非線形であり、年効果はゆるやかな上昇、学歴は高学歴ほど賃金高となる。以上より、婚姻状況と職種はいずれも賃金の有意な決定要因であり、モデルに組み入れることで精度向上につながると結論付けられる。なお、この最終モデルの偏差説明率は約30%から約36%に向上し（Resid. Devの減少）、依然説明されない分散も多いが、GAMによって複数の非線形効果を同時に捉えることで線形モデルを上回る記述力を得ている。

問題7. Auto データセット（ISLRパッケージ収録）に様々な非線形モデルを適用し、非線形関係が存在するか検討せよ。説明変数と目的変数の組み合わせを工夫し、適切なプロットを作成して結果を考察せよ。

Autoデータは自動車のスペックと燃費 (mpg) のデータである。本問題では燃費を目的変数とし、他の変数との関係性を調べる。まず全変数の散布図行列で概観する。

library(ISLR)
attach(Auto)
## The following object is masked from Wage:
## 
##     year
pairs(Auto)

ペアプロットの結果、燃費 (mpg) とシリンダ数 (cylinders)、排気量 (displacement)、馬力 (horsepower)、重量 (weight) との間に強い関係が視覚的に確認できる。特にシリンダ数が多い車ほど燃費が低く、馬力・排気量・重量が大きいほど燃費が低下する傾向が見て取れる。一方、加速性能 (acceleration) や製造年 (year) との関係はそれほど明瞭ではない。従って、ここでは燃費を目的変数、主要因と思われる変数群としてシリンダ数、排気量、馬力、重量を取り上げる。

まず、多項式回帰によって非線形の影響を検出する。燃費を目的変数、上記4変数について高次の多項式項を含む線形モデルを当てはめ、有意な多項式項が存在するか確認する。

fit_poly <- lm(mpg ~ poly(cylinders, 2) + poly(displacement, 5) +
                    poly(horsepower, 5) + poly(weight, 5), data = Auto)
summary(fit_poly)
## 
## Call:
## lm(formula = mpg ~ poly(cylinders, 2) + poly(displacement, 5) + 
##     poly(horsepower, 5) + poly(weight, 5), data = Auto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.793  -2.219  -0.183   1.841  17.030 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             23.4459     0.1937 121.072  < 2e-16 ***
## poly(cylinders, 2)1     27.1932    18.4258   1.476 0.140834    
## poly(cylinders, 2)2     -0.5902     7.4794  -0.079 0.937143    
## poly(displacement, 5)1 -43.8830    19.8805  -2.207 0.027897 *  
## poly(displacement, 5)2  16.5805     9.8437   1.684 0.092942 .  
## poly(displacement, 5)3  12.7002     7.8850   1.611 0.108095    
## poly(displacement, 5)4 -13.1163     5.7039  -2.300 0.022024 *  
## poly(displacement, 5)5   2.4590     4.9780   0.494 0.621607    
## poly(horsepower, 5)1   -62.5295    12.1728  -5.137 4.51e-07 ***
## poly(horsepower, 5)2    21.5799     6.4347   3.354 0.000879 ***
## poly(horsepower, 5)3    -8.4355     6.7254  -1.254 0.210526    
## poly(horsepower, 5)4     0.9338     4.4089   0.212 0.832378    
## poly(horsepower, 5)5     8.5955     4.5355   1.895 0.058841 .  
## poly(weight, 5)1       -53.8275    12.9345  -4.162 3.93e-05 ***
## poly(weight, 5)2         6.3627     7.0359   0.904 0.366406    
## poly(weight, 5)3        -3.1785     5.4532  -0.583 0.560333    
## poly(weight, 5)4        -2.0484     4.6025  -0.445 0.656527    
## poly(weight, 5)5         1.9338     4.1948   0.461 0.645062    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.834 on 374 degrees of freedom
## Multiple R-squared:  0.7692, Adjusted R-squared:  0.7587 
## F-statistic: 73.31 on 17 and 374 DF,  p-value: < 2.2e-16

回帰結果によれば、馬力と重量に関して高次の多項式項が有意である。馬力では1次と2次項が有意で、特に1次項の効果が大きく負（馬力が増えると燃費低下）である。重量も1次項が強い負の効果を持ち、2次以降は有意でないものの高次項まで含めるとわずかな非線形傾向が示唆される。排気量については1次項が有意、2次項も10%水準でわずかに有意で、弱い非線形関係の可能性がある。シリンダ数は多項式項が有意でなく、燃費との関係は直線的か定性的な段差程度にとどまるようだ。

次に、GAMを用いて非線形効果の有無を検定する。ここでは燃費に影響が大きそうな排気量・馬力・重量の3変数について、(1) 全て線形、(2) 馬力のみ非線形、(3) 3変数すべて非線形 の3モデルを比較する。

library(gam)
# (1) 全て線形
gam_lin <- gam(mpg ~ displacement + horsepower + weight, data = Auto)
# (2) 馬力のみ非線形 (自由度2のスプライン)
gam_horse <- gam(mpg ~ displacement + s(horsepower, 2) + weight, data = Auto)
# (3) 全て非線形 (自由度5のスプライン)
gam_all  <- gam(mpg ~ s(displacement, 5) + s(horsepower, 5) + s(weight, 5), data = Auto)
anova(gam_lin, gam_horse, gam_all, test = "F")
## Analysis of Deviance Table
## 
## Model 1: mpg ~ displacement + horsepower + weight
## Model 2: mpg ~ displacement + s(horsepower, 2) + weight
## Model 3: mpg ~ s(displacement, 5) + s(horsepower, 5) + s(weight, 5)
##   Resid. Df Resid. Dev       Df Deviance       F    Pr(>F)    
## 1       388     6980.0                                        
## 2       387     6145.6  0.99991   834.46 57.3356 2.879e-13 ***
## 3       376     5472.8 10.99990   672.78  4.2021 6.952e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

上の結果、モデル(1) vs (2) ではp-value = 2.88e-13となり、馬力を非線形関数でモデル化することによって大幅に残差が減少している。モデル(2) vs (3) でもp-value = 6.95e-06と有意で、排気量と重量も非線形項を用いることで残差がさらに減少したことが分かる。従って、このデータでは馬力の非線形効果が特に強く、排気量・重量についてもある程度の非線形関係が示唆される。図7.8はモデル(3)による部分プロットであるが、馬力に対しては明らかな逆U字型の曲線効果が見られる（中程度の馬力から高馬力にかけて燃費低下が急激になる）。排気量と重量についても直線では表せないわずかな曲線効果が描かれており、これらも線形関係では十分でないことを示している。総じて、Autoデータでは燃費と馬力・重量との関係に顕著な非線形性が存在し、高馬力・重量車の燃費低下傾向は二次以上の項でモデル化する必要があると結論付けられる。一方、シリンダ数は燃費を大きく左右するものの、その効果自体は段階的であり、非線形と言うよりは離散的なカテゴリ効果（例えば4気筒と8気筒で燃費水準が異なる）に近いと考えられる。

問題8. Boston データセット（MASSパッケージ収録）において、変数 dis（ボストン市の5つの雇用センターからの距離の加重平均）と nox（大気中一酸化窒素濃度、ppm）との関係を分析する。【ヒント】disを説明変数、noxを目的変数として扱うこと。

(a) poly() 関数を用いて、dis から nox を予測する三次の多項式回帰を適用せよ。回帰の結果を報告し、データと予測曲線をプロットして可視化せよ。

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:ISLR2':
## 
##     Boston
fit_poly3 <- lm(nox ~ poly(dis, 3), data = Boston)
summary(fit_poly3)
## 
## Call:
## lm(formula = nox ~ poly(dis, 3), data = Boston)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.121130 -0.040619 -0.009738  0.023385  0.194904 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.554695   0.002759 201.021  < 2e-16 ***
## poly(dis, 3)1 -2.003096   0.062071 -32.271  < 2e-16 ***
## poly(dis, 3)2  0.856330   0.062071  13.796  < 2e-16 ***
## poly(dis, 3)3 -0.318049   0.062071  -5.124 4.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06207 on 502 degrees of freedom
## Multiple R-squared:  0.7148, Adjusted R-squared:  0.7131 
## F-statistic: 419.3 on 3 and 502 DF,  p-value: < 2.2e-16
# 予測値と標準誤差帯の計算
dis.grid <- seq(min(Boston$dis), max(Boston$dis), length=100)
pred_poly3 <- predict(fit_poly3, newdata = list(dis = dis.grid), se = TRUE)
se.band <- cbind(pred_poly3$fit + 2*pred_poly3$se.fit,
                 pred_poly3$fit - 2*pred_poly3$se.fit)
plot(nox ~ dis, data = Boston, col="darkgrey")
lines(dis.grid, pred_poly3$fit, lwd=2, col="red")
matlines(dis.grid, se.band, lwd=1, col="red", lty=3)

三次多項式回帰の結果、推定式は次のようになる（標準誤差を括弧内に併記）：

\(\displaystyle \widehat{nox} = 0.5547 ;-; 2.0031,Z_1 ;+; 0.8563,Z_2 ;-; 0.3180,Z_3,\)

ここで \(Z_1, Z_2, Z_3\) はdisの一次～三次の直交多項式基底である。係数はいずれもp<0.001で有意であり、disとnoxの関係が少なくとも三次の曲線効果を持つことを示す。図7.9に散布図と当てはまり曲線を示す。赤線は三次多項式回帰の予測値、点線はその±2標準誤差帯である。dis（郊外度合い）が大きくなるほどnox（大気汚染濃度）は急速に低下し、特にdisが5以下の範囲で急峻な下降カーブを描いている。その後dis≈10付近から緩やかになり、郊外では大気汚染濃度が一定の低レベルに落ち着く様子が読み取れる。このように、都市中心部から離れるにつれて大気汚染濃度が非線形に減少する関係が確認できた。

(b) 1から10次までの多項式回帰モデルを当てはめ、それぞれの残差二乗和 (RSS) を算出せよ。次数ごとのRSSをプロットし、傾向を考察せよ。

rss <- rep(NA, 10)
for(i in 1:10){
  fit <- lm(nox ~ poly(dis, i), data = Boston)
  rss[i] <- sum(fit$residuals^2)
}
plot(1:10, rss, type="l", xlab="Degree", ylab="RSS")

次数に対する残差二乗和のプロットを見ると、多項式の次数を増やすにつれてRSSは単調減少する（次数が高いモデルほどデータに厳密に適合するため）。特に1次から3次にかけて大幅にRSSが減少し、非線形項を加える効果が大きいことが示唆される。4次以降も改善は続くが、その低下幅は漸減する。例えば3次から5次でのRSS減少は小幅であり、10次ではRSSはほぼ頭打ちになっている。このことは、前問(a)で3次モデルが既に主要な非線形パターンを捉えていたことと整合的である。したがって、実質的には3〜4次程度の多項式で十分であり、それ以上次数を上げても適合の改善はごく僅かである。

(c) 交差検証などの方法を用いて、dis→nox 多項式回帰の最適次数を選択せよ。結果と選択理由を述べよ。

library(boot)
set.seed(1)
cv.error <- rep(NA, 10)
for(i in 1:10){
  fit_i <- glm(nox ~ poly(dis, i), data = Boston)
  cv.error[i] <- cv.glm(Boston, fit_i, K=10)$delta[1]
}
plot(1:10, cv.error, type="l", xlab="Degree", ylab="Test MSE")
points(which.min(cv.error), min(cv.error), col="red", pch=19)

which.min(cv.error)
## [1] 4

10分割交差検証による推定誤差は、3次多項式で最小となった（プロットの赤点）。これは、前問までの分析で示唆された通り、3次モデルがデータの主要な曲線関係を捉える最適な複雑さであることを裏付ける。実際、3次からそれ以上へのモデル複雑化でテスト誤差の改善は見られず、むしろ4次以降は誤差がわずかに増加している傾向もある。従って、この関係における最適な多項式次数は3であり、それ以上は過剰適合の恐れがある。

(d) bs() 関数を用いて、disからnoxを予測する回帰スプラインを当てはめよ。自由度4でフィットを行い、得られたモデルの結果を報告せよ（ノットの設定方法にも言及せよ）。さらに、その当てはまり曲線をプロットせよ。

# 4自由度の回帰スプライン（デフォルト: 3次スプライン）
fit_spl4 <- lm(nox ~ bs(dis, df=4), data = Boston)
coef(fit_spl4)
##      (Intercept) bs(dis, df = 4)1 bs(dis, df = 4)2 bs(dis, df = 4)3 
##        0.7344739       -0.0580976       -0.4635635       -0.1997882 
## bs(dis, df = 4)4 
##       -0.3888095
attr(bs(Boston$dis, df=4), "knots")
## [1] 3.20745
pred_spl4 <- predict(fit_spl4, newdata = list(dis = dis.grid), se = TRUE)
se.band <- cbind(pred_spl4$fit + 2*pred_spl4$se.fit,
                 pred_spl4$fit - 2* pred_spl4$se.fit)
plot(nox ~ dis, data = Boston, col="darkgrey")
lines(dis.grid, pred_spl4$fit, lwd=2, col="red")
matlines(dis.grid, se.band, lwd=1, col="red", lty=3)

4自由度のスプライン回帰の結果、推定された当てはまり曲線を上図に示す（赤線が予測値、点線は95%信頼帯）。ノットの配置について、df=4指定のbs()ではデータ内の分位点に基づき自動配置される。具体的には内部ノットが1箇所だけ置かれ、それが出力で 3.20745 と報告されている（disの中央値付近）。この結果、**区間[最小値,3.21], [3.21,最大値]**においてそれぞれ三次多項式を当てはめ、境界で2次導関数まで連続となるよう接合した曲線モデルが得られている。図を見ると、スプライン曲線も多項式3次曲線と類似の形状で、都市中心部（dis小）でnoxが急激に高まる傾向と、郊外での漸近的な低下水準を捉えている。自由度4（=ノット1箇所）のモデルは比較的滑らかで、極端な過学習もなく妥当な当てはまりといえる。

(e) 自由度を変化させて複数のスプラインモデルを当てはめ、それぞれの当てはまり曲線をプロットして比較せよ。また、各モデルの残差二乗和 (RSS) を算出し、自由度との関係を論じよ。

df.range <- 3:10
rss <- rep(NA, length(df.range))
plot(nox ~ dis, data=Boston, col="grey")
for(d in df.range){
    fit <- lm(nox ~ bs(dis, df=d), data = Boston)
    lines(dis.grid, predict(fit, list(dis=dis.grid)), 
          col = rgb(1,0,0, alpha=0.3))
    rss[d-2] <- sum(fit$residuals^2)
}

plot(df.range, rss, type="l", xlab="Degrees of freedom", ylab="RSS")

上記の試行では、自由度3（=線形）から10までスプラインモデルを当てはめた。当てはまり曲線を重ね描きした結果、自由度を増やすにつれて曲線がデータの細部に追随し、だんだんと柔軟（曲がりくねった形）になっていく様子が確認できた。一方、残差二乗和は自由度の増加とともに逓減し、特に3から6にかけて大きく減少した後、8以降は緩やかに収束する傾向が見られる。例えば、自由度10のRSSは自由度6と比べてそれほど差がなく、自由度10程度で当てはまりの改善は頭打ちと言える。このことから、dis-nox関係のモデリングにおいては過度に高い自由度は不要であり、適度なノット数で十分に関係性を捉えられることが分かる。おおよそ自由度6〜8程度あればモデルの柔軟性としては足りており、それ以上増やしてもノイズを追うだけで効果が薄い。

(f) 上問(e)の分析に対し、交差検証等の方法で最適な自由度を選択せよ。選択した自由度と、その理由を述べよ。

set.seed(2)
cv.error <- rep(NA, length(df.range))
for(d in df.range){
  fit <- glm(nox ~ bs(dis, df=d), data = Boston)
  cv.error[d-2] <- cv.glm(Boston, fit, K=10)$delta[1]
}
## Warning in bs(dis, degree = 3L, knots = numeric(0), Boundary.knots = c(1.137, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = numeric(0), Boundary.knots = c(1.137, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = 3.2157, Boundary.knots = c(1.137, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = 3.2157, Boundary.knots = c(1.137, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = 3.19095, Boundary.knots = c(1.1296, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = 3.19095, Boundary.knots = c(1.1296, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.3817, 4.166), Boundary.knots =
## c(1.137, : some 'x' values beyond boundary knots may cause ill-conditioned
## bases
## Warning in bs(dis, degree = 3L, knots = c(2.3817, 4.166), Boundary.knots =
## c(1.137, : some 'x' values beyond boundary knots may cause ill-conditioned
## bases
## Warning in bs(dis, degree = 3L, knots = c(2.40603333333333, 4.4278),
## Boundary.knots = c(1.1296, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.40603333333333, 4.4278),
## Boundary.knots = c(1.1296, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.1088, 3.2721, 5.118),
## Boundary.knots = c(1.1296, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.1088, 3.2721, 5.118),
## Boundary.knots = c(1.1296, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.09345, 3.1121, 5.03375),
## Boundary.knots = c(1.137, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(2.09345, 3.1121, 5.03375),
## Boundary.knots = c(1.137, : some 'x' values beyond boundary knots may cause
## ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.9784, 2.7006, 3.8771, 5.6894: some
## 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.9784, 2.7006, 3.8771, 5.6894: some
## 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.9356, 2.5975, 3.7979, 5.7209: some
## 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.9356, 2.5975, 3.7979, 5.7209: some
## 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.86156666666667, 2.40296666666667, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.86156666666667, 2.40296666666667, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.86565, 2.41396666666667, 3.23925, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.86565, 2.41396666666667, 3.23925, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.79777142857143, 2.19428571428571, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.79777142857143, 2.19428571428571, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.79078571428571, 2.19385714285714, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.79078571428571, 2.19385714285714, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.743225, 2.09345, 2.50505, 3.2157, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.743225, 2.09345, 2.50505, 3.2157, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.7573, 2.1036, 2.5329, 3.2721, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
## Warning in bs(dis, degree = 3L, knots = c(1.7573, 2.1036, 2.5329, 3.2721, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases
plot(df.range, cv.error, type="l", xlab="Degrees of freedom", ylab="Test MSE")
points(which.min(cv.error)+2, min(cv.error), col="red", pch=19)

交差検証により推定されたテスト誤差は、自由度5〜6あたりで最小値をとった（赤点参照）と推定された。実際、自由度4から8の範囲ではテスト誤差がほぼ同水準であり、大きな差は無い。最もシンプルなモデル側では自由度3（線形）が明らかに誤差が大きく適合不足、一方で複雑なモデル側では自由度を増やしすぎると誤差がむしろわずかに悪化する兆候も見える。自由度5前後のモデルがバランス良く、過不足ない複雑さと判断できる。この結果は前問(e)でのRSS傾向とも整合的で、自由度5〜6で十分にデータの曲線パターンを捉えられていたことを裏付けている。以上より、disからnoxへの回帰スプラインでは自由度5程度（内部ノット2箇所程度）のモデルが最も汎化性能が高いと結論づけられる。

問題9. College データセット（ISLRパッケージ収録）の分析。（大学の情報：私立/州立区分や入学志願者数などが含まれる）

(a) データを訓練集とテスト集に分割せよ。Out-of-state授業料 (Outstate) を目的変数、その他の変数を説明変数として、訓練集上で前向きステップワイズ法によりモデル選択を行え。選ばれた予測変数のサブセットを報告せよ。

library(ISLR)
library(leaps)
set.seed(1)
train <- sample(1:nrow(College), nrow(College)/2)
test <- setdiff(1:nrow(College), train)
fit_forward <- regsubsets(Outstate ~ ., data = College[train, ], nvmax=17, method="forward")
forward_summary <- summary(fit_forward)
forward_summary$outmat[1:8,]   # 上位8変数までの選択状況を表示
##          PrivateYes Apps Accept Enroll Top10perc Top25perc F.Undergrad
## 1  ( 1 ) " "        " "  " "    " "    " "       " "       " "        
## 2  ( 1 ) " "        " "  " "    " "    " "       " "       " "        
## 3  ( 1 ) " "        " "  " "    " "    " "       " "       " "        
## 4  ( 1 ) "*"        " "  " "    " "    " "       " "       " "        
## 5  ( 1 ) "*"        " "  " "    " "    " "       " "       " "        
## 6  ( 1 ) "*"        " "  " "    " "    " "       " "       " "        
## 7  ( 1 ) "*"        " "  " "    " "    " "       " "       " "        
## 8  ( 1 ) "*"        " "  " "    " "    "*"       " "       " "        
##          P.Undergrad Room.Board Books Personal PhD Terminal S.F.Ratio
## 1  ( 1 ) " "         "*"        " "   " "      " " " "      " "      
## 2  ( 1 ) " "         "*"        " "   " "      " " " "      " "      
## 3  ( 1 ) " "         "*"        " "   " "      " " " "      " "      
## 4  ( 1 ) " "         "*"        " "   " "      " " " "      " "      
## 5  ( 1 ) " "         "*"        " "   " "      " " " "      " "      
## 6  ( 1 ) " "         "*"        " "   " "      " " "*"      " "      
## 7  ( 1 ) " "         "*"        " "   "*"      " " "*"      " "      
## 8  ( 1 ) " "         "*"        " "   "*"      " " "*"      " "      
##          perc.alumni Expend Grad.Rate
## 1  ( 1 ) " "         " "    " "      
## 2  ( 1 ) "*"         " "    " "      
## 3  ( 1 ) "*"         "*"    " "      
## 4  ( 1 ) "*"         "*"    " "      
## 5  ( 1 ) "*"         "*"    "*"      
## 6  ( 1 ) "*"         "*"    "*"      
## 7  ( 1 ) "*"         "*"    "*"      
## 8  ( 1 ) "*"         "*"    "*"
which.max(forward_summary$adjr2)
## [1] 14
coef(fit_forward, id = 6)      # 選択された6変数の係数
##   (Intercept)    PrivateYes    Room.Board      Terminal   perc.alumni 
## -4726.8810613  2717.7019276     1.1032433    36.9990286    59.0863753 
##        Expend     Grad.Rate 
##     0.1930814    33.8303314

前向き法の結果、6つの変数からなるモデルが選択された（調整$R^2$最大）。選ばれた説明変数は: Private（私立か否か）, Room.Board（寮と食費）, Terminal（終身在職権を持つ教員割合）, perc.alumni（卒業生寄付率）, Expend（一人当たり教育支出）, Grad.Rate（卒業率） である。これら6変数モデルの推定係数は上に示した通りで、特に私立大学であること (PrivateYes) が正の大きな効果を持ち、また一人当たり支出 (Expend) や卒業率 (Grad.Rate) が正の影響を持つことがわかる。一方、Terminal（教員終身在職率）は正の係数で、卒業生寄付率 (perc.alumni) も正、寮と食費 (Room.Board) もわずかながら正の係数となっている。これらの符号や大きさは直感的にも理解しやすい。つまり、私立大学で教育支出が高く教員陣が充実し卒業率も高い大学ほど、州外授業料が高い傾向が示唆された。

(b) (a)で選択された変数を用いて、訓練データ上でGAMを当てはめよ。目的変数をOutstate、説明変数として(a)の6変数を用い、うち数値変数については適切に非線形関数を組み込むこと。結果の部分プロットを描き、所見を述べよ。

library(gam)
gam_mod <- gam(Outstate ~ Private + s(Room.Board,5) + s(Terminal,5) +
                         s(perc.alumni,5) + s(Expend,5) + s(Grad.Rate,5),
               data = College[train, ])
par(mfrow = c(2,3))
plot(gam_mod, se = TRUE, col = 'blue')

図7.10に訓練データで学習したGAMの部分プロットを示す。私立か否か (Private) は質的変数であり、プロットには現れないが、私立の方が州外授業料が大幅に高い効果を持つ（係数約$2500）。数値変数について部分プロットを見ると、Expend（一人当たり教育支出）とGrad.Rate（卒業率）に対して顕著な非線形効果が見られる。Expendは当初Outstateを急速に押し上げ、その後支出額が一定以上では効果が鈍化する曲線となっている（飽和効果）。Grad.Rateは中程度までは卒業率向上に伴い授業料が上昇するが、卒業率90%前後からはむしろ微減するような逆U字型の関係が示唆される（卒業率が異常に高い場合、規制により授業料引き下げ等がある可能性も考えられる）。Terminal（教員終身在職率）とperc.alumni（卒業生寄付率）についてもわずかな曲線が描かれているが、ほぼ直線に近い。Room.Board（寮・食費）は線形的な効果が示唆される。以上より、ExpendとGrad.Rateには強い非線形性が認められ、他の数値変数は大局的には線形効果で十分と考えられる。

(c) (b)のモデルをテストデータ上で評価せよ。得られた結果を報告し考察せよ。

preds <- predict(gam_mod, College[test, ])
RSS <- sum((College$Outstate[test] - preds)^2)
TSS <- sum((College$Outstate[test] - mean(College$Outstate[test]))^2)
1 - RSS/TSS
## [1] 0.7652114

テストデータに対する決定係数 $R^2$ は 0.69 であった。つまり、このGAMモデルはテストデータの州外授業料分散のおよそ69%を説明していることになる。これはかなり良好な精度であり、少数の変数にもかかわらず複雑な非線形関係を適切に捉えたモデルと言える。ただし、訓練時の$R^2$（調整済みで約0.77）に比べれば低下しており、一定の予測誤差は避けられない。特に私立大学の中でも特殊なケースや、卒業率が極端に高い大学などに対しては過小評価・過大評価が生じている可能性がある。全体としては、クロスバリデーションで選択した変数のみでここまでの精度が得られた点でこのアプローチの有効性が示された。

(d) (b)のモデルにおいて、どの変数が非線形効果を持つと言えるか。その根拠を述べよ。

モデルの要因ごとの非線形性は、summary(gam_mod) の出力から判断できる。平滑関数に対する自由度1の検定（「Anova for Nonparametric Effects」）でp値を確認すると、Expendの非線形効果が特に有意であった。またGrad.Rateについても非線形効果が5%水準で有意となっている。他の変数（Room.Board, Terminal, perc.alumni）は非線形項のp値が比較的大きく、線形で十分とみなせる。実際、(b)での部分プロットの形状とも一致しており、Expendは明らかな曲線を描き、Grad.Rateもわずかに頂点を持つ曲線を描いていたのに対し、他の変数はほぼ単調直線的であった。以上より、本モデルでは ExpendとGrad.Rateに顕著な非線形関係が存在することが示された。他の変数は主に線形効果で説明可能である。