正則化・次元削減回帰 実習資料 (ISLR Chapter 6)
Yusuke Matsui
2025-07-02
目的 本資料では R コードを実際に入力しながら，Ridge／Lasso／Elastic Net に加え，主成分回帰（PCR）および部分最小二乗回帰（PLS）の理念・実装・ハイパーパラメータ調整・結果解釈を系統的に学習する。
目次
- R 文法ウォームアップ
- データ前処理
- Ridge 回帰（L2 正則化）
- Lasso 回帰（L1 正則化）
- Ridge と Lasso の比較
- Elastic Net
- 主成分回帰（PCR）
- 部分最小二乗回帰（PLS）
- まとめおよび課題
0 R 文法ウォームアップ — 5 分のスピード復習
このセクションで学ぶこと
以降の演習で迷わないよう、R の最頻出構文とデータ操作パターン を一気に復習する。
各チャンクは「▶ 入力 → ★ 標準出力」の最小例を示し、動かした直後にポイントを確認できるよう設計されている。
読み進め方の指針
- コードをコピーして RStudio Console か Quarto で実行。
- すぐ下の出力や ポイント を見て「何が起きたか」を言語化。
- 疑問があれば同じチャンク内で値を書き換え→再実行して挙動を確認。
ここで扱うのは 代入・ベクトル演算・データフレーム／tibble・パイプ・ユーザ関数・formula・設計行列 の７つ。 これらは後続の glmnet／pls など高度な関数でも “共通語” として登場するため、必ず手を動かして感覚を取り戻してほしい。
以下、各サブセクションで “コード → 出力 → ポイント” の順に進む。最初の 0-1 から 0-7 までは 1 行ずつタイプ することでショートカットキーや自動保完にも慣れておくと後が楽になる。
0-1 オブジェクトの作成と代入演算子
目的（後続との関連） R の基本データ型と代入記号
<-を確認する。以降の λ グリッド生成 や 行列・ベクトル操作 はすべてベクトルオブジェクトの作成が出発点となる。
# ベクトルを作成
x <- 1:5 # “<-” が R の標準的な代入記号
x # → 1 2 3 4 5## [1] 1 2 3 4 5ポイント
<-と=の両方で代入は可能だが、R 文化圏では<-がデフォルト。1:5は 整数シーケンス を作るショートハンド。seq(1, 5)と同義。
0-2 ベクトル演算とリサイクル規則
目的（後続との関連） ベクトル化演算 と リサイクル規則 を体験し、後で用いる
mean((y-pred)^2)などの 要素ごと演算 + 集約 の挙動を理解する。
y <- x * 2 # 要素ごとに倍
y # → 2 4 6 8 10## [1] 2 4 6 8 10x + c(1, 2) # ベクトル長が合わない場合は “リサイクル” される## Warning in x + c(1, 2): longer object length is not a multiple of shorter
## object length## [1] 2 4 4 6 6# → 2 4 4 6 6 # (1,2,1,2,1 と繰り返し)ポイント
- R では ベクトル化演算 が基本。
*や+はループ不要で全要素に作用。- 長さが合わないと短い方が繰り返される（警告が出る場合あり）— “リサイクル規則”。
0-3 tibble（近代的データフレーム）
目的（後続との関連）
tidyverse系データ構造に慣れる。後続の 欠損処理 (drop_na) や パイプ操作 は tibble が前提となる。
library(tidyverse) # dplyr 等を一括ロード## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorstbl <- tibble(a = 1:3,
b = c("A", "B", "C"))
tbl## # A tibble: 3 × 2
## a b
## <int> <chr>
## 1 1 A
## 2 2 B
## 3 3 C# A tibble: 3 × 2
a b
<int> <chr>
1 1 A
2 2 B
3 3 Cポイント
tibbleは 列型を印字 し、省略表示がスマート。tidyverse関数群（dplyr::filter()など）と高い親和性を持つ。- 行番号付与や因子自動変換がないので、古典的
data.frameより安全。
0-4 パイプ演算子 — データ加工の読みやすさ向上
目的（後続との関連） ネイティブパイプ
|>の読み順を習得し、欠損除去 → 計算 → 可視化 の一連フローを後で簡潔に記述できるようにする。
iris |>
head() |>
dplyr::summarise(across(where(is.numeric), mean))## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 4.95 3.383333 1.45 0.2333333ポイント
|>は base R 4.1+ で公式採用されたネイティブパイプ。- “左で作ったものを右の関数の 第一引数 に渡す” だけの極めてシンプルな規則。
magrittr流の%>%とほぼ同じだが、ネイティブの方が依存が少ない。
0-5 ユーザ定義関数
目的（後続との関連） 関数定義の文法を確認し、後で汎化性能を測る
test_mse()などの小さなユーティリティ関数を自作できるようにする。
square <- function(z) {
z ^ 2
}
square(4) # → 16## [1] 16ポイント
function(引数){ 処理 }で匿名関数を作り、オブジェクトに代入。- R は 第一級関数 を持つ言語：関数を変数に格納・返却・引数渡しできる。
0-6 formula（~）記法と線形回帰
目的（後続との関連） formula で 応答 ~ 説明変数 を宣言する書式を学ぶ。PCR・PLS・
model.matrix()でも同じ構文を用いるため必須知識。
fit <- lm(Sepal.Length ~ Species, data = iris)
summary(fit)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.006 0.07280222 68.761639 1.134286e-113
## Speciesversicolor 0.930 0.10295789 9.032819 8.770194e-16
## Speciesvirginica 1.582 0.10295789 15.365506 2.214821e-32ポイント
y ~ x1 + x2で 目的変数と説明変数を宣言。lm()は重回帰までシンプルに書ける。Speciesは自動でダミー変数化（因子扱い）。
0-7 model.matrix() で設計行列を直接取得
目的（後続との関連） 数値行列入力 を要求する
glmnet()に渡す設計行列を生成する手順を身につける。切片列の除外[ , -1]は正則化回帰で定石。
X <- model.matrix(Sepal.Length ~ ., iris)[, -1] # 切片列を除外
dim(X) # → 150 × 4## [1] 150 5ポイント
- glmnet 等の 行列入力専用 関数にそのまま渡せる。
- ダミー変数（0/1）や spline 基底なども formula で指定すれば自動生成。
0-8 練習 — ここまでの要素を 1 行で
目的（後続との関連） ここまで習得した パイプ・グループ化・無名関数・ユーザ関数 を複合し、後で行う モデル評価ループ の下地となるコード読解力を仕上げる。
iris |>
dplyr::group_by(Species) |>
dplyr::summarise(across(where(is.numeric), ~ square(mean(.x))))## # A tibble: 3 × 5
## Species Sepal.Length Sepal.Width Petal.Length Petal.Width
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 setosa 25.1 11.8 2.14 0.0605
## 2 versicolor 35.2 7.67 18.1 1.76
## 3 virginica 43.4 8.84 30.8 4.10各種テクニック（パイプ・グループ化・ユーザ関数・無名関数
~）を複合。 実行して出力を確認し、構文を身体で思い出そう。
1 準備 & データ確認
1-1 パッケージ読み込み
目的 必要パッケージをロードし，乱数シードを固定して解析の再現性を確保する。
処理内容
- （事前）
install.packages()で依存パッケージを取得。library()で各パッケージの名前空間を登録。set.seed()で疑似乱数生成器 (RNG) の状態を固定。使用関数（抜粋）
library(<pkg>)— パッケージ<pkg>をロードset.seed(integer)— RNG 状態を固定結果の読み取り方
library()がエラーなく完了すればロード成功。 同じシードを用いることで，本資料の数値結果を受講者側でも再現できる。
# install.packages(c("ISLR2", "glmnet", "pls", "tidyverse"))
library(ISLR2) # 演習データ集（ISLR 第 2 版）
library(glmnet) # 正則化回帰（Ridge / Lasso / Elastic Net）## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-8library(pls) # PCR・PLS 回帰##
## Attaching package: 'pls'## The following object is masked from 'package:stats':
##
## loadingslibrary(tidyverse) # dplyr, ggplot2, tidyr など一括ロード
set.seed(1)1-2 データのクリーニング
目的
- 解析前に 欠損状況を可視化 し，どの変数に何件の NA が含まれるかを把握する。
- 応答変数 Salary に NA を持つレコードを除外し，完全データセットを作成する。
処理内容 Step 0 欠損チェック
colSums(is.na(Hitters))で各列の欠損件数を列挙し，特にSalaryに NA が含まれることを確認する。 Step 1 欠損除去drop_na(Salary)で Salary 列に NA を含む行 を削除。 Step 2 除去後の寸法確認dim()で行数・列数を確認し，欠損除去が想定どおり行われたか検証する。使用関数（抜粋）
is.na(x)— 要素が NA か TRUE/FALSE を返すcolSums(mat)— 列方向に合計（TRUE は 1 として加算）drop_na(cols)— 指定列で NA を持つ行を削除（tidyr）dim(object)— 行数・列数を整数ベクトルで返す
結果の読み取り方
colSums(is.na(ISLR2::Hitters))## AtBat Hits HmRun Runs RBI Walks Years CAtBat ## 0 0 0 0 0 0 0 0 ## CHits CHmRun CRuns CRBI CWalks League Division PutOuts ## 0 0 0 0 0 0 0 0 ## Assists Errors Salary NewLeague ## 0 0 59 0の出力で
Salaryが 59 件の NA を含むことを確認。 つづいてdim(hitters)が263 20ならば 322 − 59 = 263 行・20 列（応答 1 + 説明 19）の完全データが得られた。
## Step 0 : 欠損チェック ---------------------------------------------
na_count <- colSums(is.na(ISLR2::Hitters))
na_count[na_count > 0] # 欠損のある列だけ表示## Salary
## 59## Step 1 : Salary 欠損を除去 --------------------------------------
hitters <- ISLR2::Hitters |>
drop_na(Salary)
## Step 2 : 除去結果を確認 -----------------------------------------
dim(hitters) # 263 20## [1] 263 201-3 行列入力を準備
目的
glmnet()が受け付ける 数値行列xと 数値ベクトルyを作成する。処理内容
model.matrix()でダミー変数化済み設計行列を生成。- 1 列目（Intercept）を除外して説明変数行列
xとする。- データフレーム列
Salaryをベクトル抽出して応答yとする。使用関数（抜粋）
model.matrix(formula, data)— 設計行列 (X) を返す[, -1]— 行列/データフレームの列スライス$— データフレーム列アクセス結果の読み取り方
x： 263 行 × 19 列の数値行列y： 長さ 263 の数値ベクトル注意
glmnet()はxに数値行列（matrixあるいはdgCMatrix）を要求する点に留意。
x <- model.matrix(Salary ~ ., hitters)[, -1] # 263 × 19
y <- hitters$Salary # 長さ 263注意
glmnet()は引数xに数値行列（matrixあるいはdgCMatrix）を要求する。
補足：model.matrix() と formula 記法の基礎
1.
model.matrix()とは
- 設計行列 (design matrix, \(X\)) を生成 する関数。
- 引数は
formulaとdata。- 因子（カテゴリカル変数）を ダミー変数（0/1 ベクトル） に自動変換。
- 切片列（全て 1 の列）を既定で先頭に付与（
-1または0 +で削除可能）。- 交互作用項 (
x1:x2) や高次項 (I(x^2)) も formula で指定すれば列を追加。2. formula 記法の最小セット
構文 意味 例における列追加 y ~ .右辺の .は「データ内の全列」を展開Hitters の 19 説明変数 +列を追加 ~ x1 + x2-列を除外 ~ . - CHits（CHits 列を外す）-1または0 +切片を除外 ~ x1 + x2 - 1x1:x2交互作用 Speed:Spin列を生成I(x^2)変数変換（数学演算） Age^2列を生成3. 回帰モデルとの対応 線形回帰の一般形は
\[ y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \varepsilon,\qquad \mathbf{y} = X\boldsymbol{β} + \boldsymbol{ε} \]
ここで \(X\) が
model.matrix()で生成した行列になる。
lm()など formula 入力 を受け付ける関数は内部でmodel.matrix()を呼び出し自動で \(X\) を構築。glmnet()は 行列入力 専用のため，事前に自分でmodel.matrix()を実行し，切片列を外す（[, -1]）必要がある。4. 実行して確かめる
mm <- model.matrix(Salary ~ CHits + CRuns, data = Hitters)
head(mm, 3)## (Intercept) CHits CRuns
## -Alan Ashby 1 835 321
## -Alvin Davis 1 457 224
## -Andre Dawson 1 1575 828
- 先頭列が切片、続く列が
CHits,CRuns。- 因子を含めば自動でダミーダミー列が追加される。
以上を踏まえ，
glmnet()用のxはmodel.matrix(式, data)[, -1]で切片列を除外して作成するのが定石である。
2 Ridge Regression（L2 ペナルティ）
2-1 λ グリッドの設定
目的 L2 正則化パスを探索するため、広い対数レンジで λ 値（正則化強度）を準備する。
処理内容
seq(10, -2, length = 100)で 10 から −2 まで等間隔の 100 点を生成し、これを底 10 の指数に変換することで 10¹⁰ 〜 10⁻² を得る。使用関数
seq(from, to, length)：指定区間を等間隔ベクトルに変換^：ベクトル化された累乗演算子結果の読み取り 得られる
lambda_gridは長さ 100 の数値ベクトル。値が大→小へ連続的に減少していることをhead()/tail()で確認するとよい。
lambda_grid <- 10 ^ seq(10, -2, length = 100)2-2 モデル学習
目的 Ridge 回帰（α = 0）を λ グリッド全体でフィッティングし、係数パスを取得する。
処理内容
glmnet()に説明行列x、応答y、α = 0、およびlambda_gridを与えて学習。使用関数
glmnet(x, y, alpha, lambda)：Elastic-Net 正則化付き線形モデルを一括学習結果の読み取り 返り値
ridge_fitは S3 クラス “glmnet”。coef()で λ ごとの係数行列、ridge_fit$lambdaで実際に学習された λ が確認できる。
ridge_fit <- glmnet(x, y, alpha = 0, lambda = lambda_grid) # α = 0 → Ridge2-3 係数パスの可視化
目的 λ の増減に伴う各係数の shrink 挙動を視覚的に確認する。
処理内容
plot()メソッド（glmnet 用）を呼び出し、x 軸に log(λ)、y 軸に係数を描画。label = TRUEで変数名を末端に付与。使用関数
plot.glmnet(object, xvar, label)：glmnet オブジェクトの係数パス描画結果の読み取り 右端（大きい λ）では全係数 ≈ 0、左へ行くほど係数の絶対値が大きくなる。曲線が互いに滑らかで交差が少ないほど多重共線性が軽微であることを示唆する。
plot(ridge_fit, xvar = "lambda", label = TRUE)
2-4 交差検証による λ の選択
交差検証 (Cross-Validation) をゼロから 交差検証とは，「データを k 個の等分割（fold）に分け，そのうち 1 つをテスト，残り k–1 個を訓練として回帰を学習・評価する」手続きを k 回 繰り返し，平均テスト誤差を汎化誤差の推定量とする方法である。
- 目的：未知データに対する予測性能（汎化性能）を事前に見積もる。
- 仕組み：各観測が 1 回だけテストデータに回るため，サンプルの使い回し効率が高い。
- 10-fold CV：デフォルトで k = 10。経験則としてデータサイズが 100 以上あればバランスが良い。
- λ の選択：CV 平均 MSE が最小の λ（
lambda.min）と，その ±1 標準誤差以内で最も単純な λ（lambda.1se）を採用することが多い。後者は過学習をさらに抑制し，係数がより縮小される。
図 各色が異なるテスト fold を示す。1 回の学習では 1 色部分のみテスト，残り 9 色が訓練となる。
目的 汎化誤差を最小化する λ（
lambda.min）と 1SE ルールの λ（lambda.1se）をデータ駆動で決定する。処理内容
cv.glmnet()で 10-fold CV を実行し、各 λ に対する平均 MSE と標準誤差を算出。プロットで CV 曲線を確認。使用関数
cv.glmnet(x, y, alpha)：k-fold CV に基づく λ 最適化plot.cv.glmnet()：CV 曲線の可視化結果の読み取り
- λmin — CV 誤差が最小
- λ1se — 最小誤差＋1SE 以内で最も単純（係数縮小が大きい） 実務ではモデル簡潔性を優先して λ1se を選択する場面も多い。
cv_ridge <- cv.glmnet(x, y, alpha = 0) # 既定で 10-fold CV
plot(cv_ridge)
ridge_lambda_min <- cv_ridge$lambda.min # 最小 MSE
ridge_lambda_1se <- cv_ridge$lambda.1se # 1SE ルール2-5 テスト MSE の評価
目的 学習済み Ridge モデルの汎化性能をホールドアウト法で定量評価する。
処理内容
- サンプルの 50 % を訓練データにランダム抽出。
- 補集合をテストデータとし、選択した λ で予測。
- 平均二乗誤差 (MSE) を計算。
使用関数
sample(n, size)：乱数サンプリングpredict(glmnet_obj, s, newx)：指定 λ で予測mean()：要素平均結果の読み取り 出力されるスカラーがテスト MSE。Lasso や Elastic Net と比較し、最も小さい手法が “そのデータにおけるベースライン” となる。
set.seed(123)
train_id <- sample(nrow(x), nrow(x) * 0.5)
test_mse <- function(fit, lambda, test_id) {
pred <- predict(fit, s = lambda, newx = x[test_id, ])
mean((y[test_id] - pred) ^ 2)
}
test_id <- setdiff(seq_len(nrow(x)), train_id)
test_mse(ridge_fit, ridge_lambda_min, test_id)## [1] 120512.4演習 1（5 分） λ を 1/10 倍および 10 倍に変更し，テスト MSE の感度を確認せよ。
3 Lasso（L1 ペナルティ）
3-1 モデル学習
目的 L1 正則化（Lasso）を λ グリッド上で学習し，係数がどの λ で 0 へ縮退するかを可視化する。
処理内容
glmnet()にalpha = 1を指定し Lasso モデルを一括学習する。plot()メソッドを呼び出し，横軸：log λ・縦軸：係数 のパス図を描画する。使用関数
glmnet(x, y, alpha, lambda)— Elastic-Net 系モデルを学習plot.glmnet(object, xvar, label)— 係数パスを可視化結果の読み取り 右端（大きい λ）では全係数 ≈ 0，左へ移るにつれ重要変数の係数が順に非 0 となる。非 0 係数が”何本残るか”が Lasso の変数選択効果となる。
lasso_fit <- glmnet(x, y, alpha = 1, lambda = lambda_grid) # α = 1 → Lasso
plot(lasso_fit, xvar = "lambda", label = TRUE)
3-2 交差検証で λ を選択
目的 10-fold 交差検証で汎化誤差最小の λ（
lambda.min）を決定する。処理内容
cv.glmnet()にalpha = 1を渡し，λ ごとの平均 CV MSE と標準誤差を計算する。使用関数
cv.glmnet(x, y, alpha)— k-fold CV に基づく λ 最適化結果の読み取り
lambda.minが最小 MSE を与える λ，plot(cv_lasso)で λ–MSE 曲線を確認できる。
cv_lasso <- cv.glmnet(x, y, alpha = 1)
lasso_lambda_min <- cv_lasso$lambda.min3-3 非ゼロ係数の抽出
目的
lambda.minに対応する係数ベクトルから 0 でない項目 を抽出し，Lasso が残した変数を特定する。処理内容
coef()で係数ベクトルを取得。- 値が 0 でない要素のみをフィルタし，名前付きベクトル
nonzeroを生成する。使用関数
coef(cv_obj, s)— 指定 λ の係数を抽出結果の読み取り
length(nonzero)が採択された変数数。Ridge（全係数 ≠ 0）と比較すれば，Lasso による変数選択効果が定量的に把握できる。
lasso_coef <- coef(cv_lasso, s = "lambda.min")
nonzero <- lasso_coef[lasso_coef != 0]
nonzero## [1] 123.7520756 -1.5473426 5.6608972 4.7296908 -9.5958375
## [6] 0.5108207 0.6594856 0.3927505 -0.5291586 32.0650811
## [11] -119.2990171 0.2724045 0.1732025 -2.0585083演習 2
length(nonzero)を報告し，Ridge の係数と比較して Lasso がどの変数を”切り捨てた”か考察せよ。
3-4 テスト MSE の計算
目的 ホールドアウト集合で Lasso モデルの汎化性能（MSE）を算出し，他手法と比較する。
処理内容
test_mse()ヘルパー関数に Lasso 学習器・lambda.min・テスト ID を渡し，平均二乗誤差を取得する。使用関数
predict(glmnet_obj, s, newx)— 指定 λ で予測mean()— 誤差平方を平均して MSE を計算結果の読み取り 出力スカラーが Lasso のテスト MSE。Ridge や Elastic Net の MSE と並べ，最小値のモデルを選択基準とする。
mse_lasso <- test_mse(lasso_fit, lasso_lambda_min, test_id)
mse_lasso## [1] 117670.34 Ridge と Lasso の比較
目的 データのスパース性・効果サイズに応じた正則化手法の選択指針を整理し，実務での意思決定を支援する。
処理内容 テーブル に「データ特性 ↔︎ 推奨手法 ↔︎ 根拠」をまとめ，簡潔なルール・オブ・サムへ昇華する。
使用要素 Markdown テーブル（
\|区切り）結果の読み取り方 表を参照し，自データの状況に最も近い行を選択する。Elastic Net は Ridge と Lasso の中庸であり，α と λ を同時に交差検証することで最適化が可能。
| データ特性 | 推奨手法 | 根拠 |
|---|---|---|
| 効果のある変数が 少数 | Lasso | L1 ペナルティにより多数の係数が 0 となり，モデル解釈性が高い |
| 多数の変数が 弱く作用 | Ridge | 係数を均一に縮小し，分散を抑制する効果が高い |
| データ構造が不明／両者の利点を折衷したい | Elastic Net（α と λ を同時最適化） | L1・L2 のハイブリッドでスパース性と安定性を両立；CV で α を学習可能 |
5 Elastic Net
5-1 α の設定例
目的 Ridge（α = 0）と Lasso（α = 1）の中間となる α = 0.5 を指定し，Elastic Net モデルを学習する。
処理内容
glmnet()にalpha = 0.5を渡し，既定の λ シーケンスでフィッティング。使用関数
glmnet(x, y, alpha)— α が 0–1 の範囲で Elastic Net を実行。結果の読み取り方 返り値
enet_fitは係数パスを含む “glmnet” オブジェクト。coef(enet_fit, s = enet_fit$lambda[10])などで任意 λ の係数を確認できる。
enet_fit <- glmnet(x, y, alpha = 0.5) # α = 0.5 で Elastic Net5-2 caret による α, λ の同時探索（参考コード）
目的 交差検証を用いて α と λ を同時に最適化 し，汎化誤差を最小化する組み合わせを取得する。
処理内容
expand.grid()で α–λ グリッドを作成。trainControl(method = "cv", number = 10)で 10-fold CV を設定。train()でmethod = "glmnet"を指定し，グリッドサーチを実行。使用関数
caret::train()— 多様なモデリング手法を統一 API で学習・チューニングtrainControl()— 交差検証設定expand.grid()— パラメータグリッド生成結果の読み取り方
caret_fit$bestTuneに最良の α・λ が格納される。caret_fit$resamplesを参照すれば，各パラメータ組における CV 結果を比較可能。
# library(caret)
# grid <- expand.grid(alpha = seq(0, 1, 0.1),
# lambda = 10 ^ seq(2, -2, length = 50))
# ctrl <- trainControl(method = "cv", number = 10)
# caret_fit <- train(x, y, method = "glmnet",
# trControl = ctrl, tuneGrid = grid)
# caret_fit$bestTune # 最適 α・λ演習 3 上記の α–λ グリッドサーチを実行し，α ∈ {0, 0.5, 1} のうちテスト MSE を最小とする組み合わせを報告せよ。
6 主成分回帰（PCR）
はじめに — “似た説明変数が多すぎる” ときの抜本策 野球選手の成績を例にすると，Hits・Runs・RBIs のように “成績が良い選手ほど全部高くなる” 変数はたくさん存在する。 こうした 相関がきつい変数を同時に回帰に入れる と，係数が不安定になり「どの変数が効いているのか」も判別しづらくなる――これが 多重共線性 と呼ばれる問題である。
主成分回帰 (Principal Components Regression, PCR) は，
相関の強い変数群をまとめて “新しい軸” (主成分) に回転し，
上位 k 本の主成分だけを説明変数として使って回帰する ――という二段構えで，多重共線性を丸ごと緩和しつつ 次元も圧縮 する方法である。
- 絵でいうと：データの点群を “雲” と見立て，その雲に最も沿った向きを新しい X 軸（PC1），その直交方向を Y 軸（PC2）… に取り直すイメージ。
- メリット：上位主成分がデータ分散を 80–90 % 説明するなら，元の 19 変数を 3〜4 軸に圧縮しても情報はほぼ保たれる。
- デメリット：主成分は元変数の線形結合なので “変数名” がなくなり，解釈がやや抽象的になる。
ゴール は「必要最小限の主成分数 k をクロスバリデーションで選び，汎化誤差を最小にする」ことである。
背景 主成分回帰は説明変数を主成分空間に射影し，多重共線性を緩和しつつ次元削減を図る。
6-1 モデル学習と交差検証
目的 主成分数ごとの交差検証 MSEP を算出し，最適次元を探索する。
処理内容
pcr()にvalidation = "CV"を指定して 10-fold CV を実行。summary()で各主成分の累積寄与率を確認。使用関数
pcr(formula, data, scale, validation)— PCR モデルを学習summary(pcr_obj)— 寄与率・CV 結果を要約結果の読み取り方 累積寄与率が 80–90 % に達する主成分数が一つの目安。 ただし最終的な選択は CV MSEP に基づく。
set.seed(2)
pcr_fit <- pcr(Salary ~ ., data = hitters, scale = TRUE,
validation = "CV")
summary(pcr_fit) # 主成分ごとの寄与率## Data: X dimension: 263 19
## Y dimension: 263 1
## Fit method: svdpc
## Number of components considered: 19
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 452 351.9 353.2 355.0 352.8 348.4 343.6
## adjCV 452 351.6 352.7 354.4 352.1 347.6 342.7
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 345.5 347.7 349.6 351.4 352.1 353.5 358.2
## adjCV 344.7 346.7 348.5 350.1 350.7 352.0 356.5
## 14 comps 15 comps 16 comps 17 comps 18 comps 19 comps
## CV 349.7 349.4 339.9 341.6 339.2 339.6
## adjCV 348.0 347.7 338.2 339.7 337.2 337.6
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 38.31 60.16 70.84 79.03 84.29 88.63 92.26 94.96
## Salary 40.63 41.58 42.17 43.22 44.90 46.48 46.69 46.75
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 96.28 97.26 97.98 98.65 99.15 99.47 99.75
## Salary 46.86 47.76 47.82 47.85 48.10 50.40 50.55
## 16 comps 17 comps 18 comps 19 comps
## X 99.89 99.97 99.99 100.00
## Salary 53.01 53.85 54.61 54.616-2 最適主成分数の選択
目的 交差検証曲線 (MSEP vs. 主成分数) を可視化し，最小 MSEP 点を特定する。
処理内容
validationplot()で CV 誤差をグラフ表示。使用関数
validationplot(pcr_obj, val.type = "MSEP")— CV 指標を描画結果の読み取り方 縦軸：MSEP 横軸：主成分数。最小値を与える成分数がモデルの候補。 過学習を避けるため “肘” の位置（誤差が下げ止まる点）を採用しても良い。
validationplot(pcr_fit, val.type = "MSEP")
演習 4 最小 MSEP を与える主成分数を決定し，そのモデルのテスト MSE を算出せよ。
6-3 テスト MSE の評価例
目的 ホールドアウト法で PCR モデルの汎化性能を評価する。
処理内容
- データの 50 % を訓練用にサンプリング。
- 訓練集合で CV 付き再学習 (
subset = train_id)。- 任意主成分数
ncompでテスト集合を予測し MSE を算出。使用関数
sample()— ランダム抽出pcr(..., subset = )— 部分集合で再学習predict(pcr_obj, newdata, ncomp)— 指定成分数で予測mean()— MSE 計算結果の読み取り方 出力 MSE を Lasso・Ridge・PLS と比較し，最も小さい手法を選択基準とする。
set.seed(123)
train_id <- sample(nrow(hitters), nrow(hitters) * 0.5)
test_id <- setdiff(seq_len(nrow(hitters)), train_id)
pcr_fit_cv <- pcr(Salary ~ ., data = hitters, subset = train_id,
scale = TRUE, validation = "CV")
# 例: CV 図より 5 主成分を採用
y_pred <- predict(pcr_fit_cv, hitters[test_id, ], ncomp = 5)
mean((y_pred - y[test_id]) ^ 2)## [1] 137597.87 部分最小二乗回帰（PLS）
はじめに — “応答変数と最も仲良しな軸を探す”
多数の説明変数が 互いにも応答にも強い相関 を持つ場合，
- PCR は「説明変数どうしの分散」を最優先に主成分を構築するため，応答に無関係な方向も保持しがちである。
- PLS は「説明変数と応答の共分散」を最大化する軸 (latent component) を逐次抽出することで，応答に直結した低次元表現 を学習する。
イメージとしては，点群 \((x_1, x_2, y)\) を 3D 空間に配置し，そのうち 「最も \(y\) と一緒に動く方向」 を第 1 成分（PLS1）として選び，残差空間に対して同様の基準で第 2 成分（PLS2）… を定義する操作に相当する。 このため PLS は 少数の成分 でも応答変動を高効率で説明できる傾向があり，とりわけ説明変数が多く，かつサンプルサイズが限られる状況で威力を発揮する。
図 第 1 PLS 成分 \(t_1\) と応答 \(y\) が高い直線相関を示す。 主成分回帰では「データ分散」を最大化する軸を選ぶのに対し，PLS は「応答との共分散」を直接最大化する点が本質的な違いである。
背景 PLS は説明変数と応答変数の 共分散 を最大化する潜在成分を逐次抽出するため，応答に直接関連する低次元表現を学習できる。
7-1 モデル学習と交差検証
目的 PLS 成分数ごとの CV MSEP を算出し，最適成分数を探索する。
処理内容
plsr()にvalidation = "CV"を指定し 10-fold CV を実施。使用関数
plsr(formula, data, scale, validation)— PLS モデルを学習summary(pls_obj)— 成分寄与・CV 結果を要約結果の読み取り方 初期成分で CV 誤差が急減し，その後飽和する傾向があれば，肘点を採用するのが実務的。
set.seed(1)
pls_fit <- plsr(Salary ~ ., data = hitters, scale = TRUE,
validation = "CV")
summary(pls_fit)## Data: X dimension: 263 19
## Y dimension: 263 1
## Fit method: kernelpls
## Number of components considered: 19
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 452 348.7 347.8 349.0 350.7 358.3 357.6
## adjCV 452 348.3 347.0 347.9 349.4 356.0 354.9
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 351.8 346.7 347.7 347.2 344.7 344.9 347.3
## adjCV 349.5 344.6 345.5 345.1 342.7 342.8 345.0
## 14 comps 15 comps 16 comps 17 comps 18 comps 19 comps
## CV 347.1 348.2 348.3 347.8 347.7 349.4
## adjCV 344.8 345.8 345.9 345.4 345.3 346.9
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 38.08 51.03 65.98 73.93 78.63 84.26 88.17 90.12
## Salary 43.05 46.40 47.72 48.71 50.53 51.66 52.34 53.26
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 92.92 95.00 96.68 97.68 98.22 98.55 98.98
## Salary 53.52 53.77 54.04 54.20 54.32 54.47 54.54
## 16 comps 17 comps 18 comps 19 comps
## X 99.24 99.71 99.99 100.00
## Salary 54.59 54.61 54.61 54.617-2 最適成分数の選択
目的 PLS の CV MSEP 曲線を可視化し，最小誤差成分数を決定する。
処理内容
validationplot()をval.type = "MSEP"で描画。使用関数
validationplot(pls_obj, val.type = "MSEP")結果の読み取り方 PCR と同様に最小 MSEP 点または”誤差低下が頭打ちになる点”が候補。
validationplot(pls_fit, val.type = "MSEP")
7-3 テスト MSE の評価例
目的 ホールドアウト集合で PLS モデルの汎化性能を算出し，PCR と比較する。
処理内容
- 6-3 と同じ
train_id/test_idを使用。- 訓練集合で PLS を再学習 (
subset = train_id)。- 最適成分数
ncompを指定して予測し MSE を算出。使用関数
plsr()・predict()・mean()— 手順は PCR と同様結果の読み取り方 得られたテスト MSE を PCR の値と並べ，どちらが小さいか判断する。
pls_fit_cv <- plsr(Salary ~ ., data = hitters, subset = train_id,
scale = TRUE, validation = "CV")
# 例: 1 成分が最適と仮定
pred_pls <- predict(pls_fit_cv, hitters[test_id, ], ncomp = 1)
mean((pred_pls - y[test_id]) ^ 2)## [1] NA演習 5 最小 MSEP を与える成分数を用いて PCR・PLS のテスト MSE を比較し，優劣を論じよ。
8 まとめおよび課題
要点整理
- 正則化回帰（Ridge / Lasso / Elastic Net）は汎化誤差の低減と変数選択に有効。
- 次元削減回帰（PCR / PLS）は多重共線性の緩和と高次元データへの適用で有用。
- 手法選択は データ構造・解釈性要件・交差検証性能 の三点で総合判断する。
総合課題
ISLR2::Creditデータで Lasso・PCR・PLS を実装し，テスト MSE を比較せよ。- p ≫ n のシミュレーションを構築し，Ridge・Lasso・PCR・PLS を横断比較せよ。