からだにいいもの Rで解析：モデルの応答と残差を可視化できる「plotmo」パッケージの紹介
機械学習や統計モデルを構築後、各予測変数が予測値にどのように影響しているかを確かめたい場面があります。しかし、背景となる変数を固定しながら応答曲線を1つずつ描き、さらに残差の分布まで確認するには手間がかかります。
「plotmo」パッケージは、線形モデルや決定木、勾配ブースティングなど幅広いモデルを対象に、予測変数ごとの応答曲線と2変数の交互作用面を描画する機能を提供します。背景変数を中央値で固定する方法と、影響を平均化する部分依存プロットの両方に対応し、残差プロットや累積分布、QQプロットのプロットも可能です。また、勾配ブースティングモデルの学習曲線や、glmnetモデルの係数経路の描画も可能です。本パッケージの利用で、モデルの挙動を確認しながら解釈を進める作業ができるのではないかと考えます。
パッケージバージョンは3.7.1。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。
＜おすすめのRに関する書籍です＞
RとStanではじめる 心理学のための時系列分析入門 (KS心理学専門書)
心理学の縦断データ・時系列データを扱うための一冊。RとStanで実装しながら学ぶ。
www.amazon.co.jp
パッケージのインストール
下記コマンドを実行してください。
# パッケージのインストール
install.packages("plotmo")
# パッケージの読み込み
library("plotmo")
# 必要に応じて
install.packages("gbm")
install.packages("glmnet")後半で扱うplot_gbmコマンドはgbmパッケージ、plot_glmnetコマンドはglmnetパッケージを利用します。いずれもplotmoの推奨パッケージのため、必要に応じてinstall.packages(“gbm”)、install.packages(“glmnet”)で別途導入してください。
スポンサーリンク
コマンド例
詳細はコメント、パッケージのヘルプを確認してください。
plotmoコマンドは、着目する予測変数だけを変化させ、それ以外の変数（背景変数）を代表値に固定して応答を描画します。既定では背景変数を中央値（因子は最頻カテゴリ）に固定しますが、pmethod=”partdep”を指定すると、背景変数の影響を全データで平均化した部分依存プロットになります。1つの変数の主効果を描いた図をdegree1プロット、2つの変数の交互作用を描いた面をdegree2プロットと呼びます。
共通データとして、京都府の宇治茶園を想定した架空の一番茶収量データを作成します。春季の平均気温、月あたりの日照時間、茶園の標高、窒素施肥量、栽培品種から、10アールあたりの収量（kg）を予測する構成です。
# 乱数の種を固定
set.seed(20260903)
# 標本サイズ
n <- 240
# 栽培品種（やぶきた・ごこう・さみどり）
hinshu <- factor(sample(c("yabukita", "gokou", "samidori"), n, replace = TRUE))
# 春季の平均気温（度）
kion <- round(runif(n, 11, 18), 1)
# 月あたりの日照時間（時間）
nissho <- round(runif(n, 110, 220))
# 茶園の標高（メートル）
hyoko <- round(runif(n, 30, 450))
# 窒素施肥量（10アールあたりkg）
hiryo <- round(runif(n, 18, 60))
# 品種ごとの収量への上乗せ量
hinshu_kou <- c(yabukita = 0, gokou = 14, samidori = 6)[as.character(hinshu)]
# 10アールあたりの収量（kg）を生成
shukaku <- 40 +
6.0 * (kion - 14) -
0.18 * (kion - 14)^2 +
0.10 * (nissho - 165) +
0.004 * (kion - 14) * (nissho - 165) -
0.02 * (hyoko - 200) / 10 +
0.9 * hiryo -
0.010 * hiryo^2 +
hinshu_kou +
rnorm(n, 0, 6)
# データフレームにまとめる
chaen <- data.frame(shukaku, kion, nissho, hyoko, hiryo, hinshu)予測変数ごとの応答曲線のプロット：plotmoコマンド
モデルオブジェクトを渡すと、予測変数ごとのdegree1プロットと、選ばれた組み合わせのdegree2プロットを1ページにまとめて描画します。pmethod=”partdep”で部分依存プロットに切り替えられます。
| オプション | 意味 | 初期値 |
|---|---|---|
| object | モデルオブジェクト | stop(“no ‘object’ argument”) |
| type | predictに渡す予測タイプ、既定はモデルに応じて自動選択 | NULL |
| nresponse | predictが複数列を返す場合に使用する列の番号または名前 | NA |
| pmethod | 描画方式、”plotmo”は背景変数を中央値で固定、”partdep”は部分依存、”apartdep”は近似部分依存 | “plotmo” |
| pt.col | 応答点の色、0で非表示 | 0 |
| jitter | pt.col指定時に応答点へ加える散らばりの大きさ、0で無効 | 0.5 |
| smooth.col | 応答点を通す平滑線の色、0で非表示 | 0 |
| level | predictが対応する場合に描く信頼区間・予測区間の水準、0で非表示 | 0 |
| func | 各図に重ねて描く関数func(x) | NULL |
| inverse.func | 描画前に応答へ適用する関数、変換した応答を元の尺度へ戻す用途 | NULL |
| nrug | 図の下端に描くラグの本数、0で非表示、”density”で密度表示 | 0 |
| grid.col | degree1プロットの背景グリッドの色、0で非表示 | 0 |
| type2 | degree2プロットの種類、”persp”・”image”・”contour”のいずれか | “persp” |
| degree1 | 描画するdegree1プロットを選ぶ添字または変数名、0で非表示 | TRUE |
| all1 | TRUEでplotmoが通常選ぶ変数に限らず全予測変数を描画 | FALSE |
| degree2 | 描画するdegree2プロットを選ぶ添字または変数名 | TRUE |
| all2 | TRUEで全ての予測変数の組を描画 | FALSE |
| do.par | 作図パラメータの設定方法、NULL・FALSE・TRUE・2のいずれか | TRUE |
| clip | TRUEで大きく外れた予測値を無視してylimを決める | TRUE |
| ylim | 全図で共通のY軸範囲、NULLで自動、NAで各図個別、c(ymin,ymax)で指定 | NULL |
| caption | 図全体の見出し、既定は自動生成、””で非表示 | NULL |
| trace | 実行時の情報表示の詳しさ、-1でplotmoの通常メッセージを抑制 | 0 |
| grid.func | 軸に無い変数を固定する値を求める関数、既定は中央値と最頻値 | NULL |
| grid.levels | 背景変数ごとに固定値を指定するリスト、grid.funcより優先 | NULL |
| extend | 各図で横軸を広げる割合、0で広げない | 0 |
| ngrid1 | 各degree1プロットで用いる等間隔なx値の個数 | 50 |
| ngrid2 | degree2プロットの格子の細かさ、ngrid2×ngrid2点を描画 | 20 |
| ndiscrete | 異なる値がこの個数以下の変数を階段状に描画 | 5 |
| npoints | 描画する応答点の個数、pt.col指定時のみ有効、TRUEまたは-1で全点 | 3000 |
| center | TRUEで描画する応答を中央化 | FALSE |
| xflip | TRUEでX軸の向きを反転 | FALSE |
| yflip | TRUEでdegree2プロットのY軸の向きを反転 | FALSE |
| swapxy | TRUEでdegree2プロットのX軸とY軸を入れ替え | FALSE |
| int.only.ok | TRUEで切片のみのモデルでも最初の予測変数で1枚だけ描画 | TRUE |
| … | plotやpersp、predictなどに渡すその他の引数 | なし |
# 気温と日照の交互作用項を含む線形モデルを構築
lm_model <- lm(shukaku ~ kion * nissho + hyoko + hiryo + hinshu, data = chaen)
# 予測変数ごとの応答曲線と交互作用面を描画
plotmo(lm_model)
plotmo grid: kion nissho hyoko hiryo hinshu
14.1 166 254 39 gokou
# 背景変数の影響を平均化した部分依存プロットを描画
plotmo(lm_model, pmethod = "partdep", degree2 = FALSE)
calculating partdep for kion
calculating partdep for nissho
calculating partdep for hyoko
calculating partdep for hiryo
calculating partdep for hinshu
＜おすすめのRに関する書籍です＞
サラっとできる!フリー統計ソフトEZR(Easy R)でカンタン統計解析
EZRでの統計解析をコンパクトにまとめた一冊。必要な手順だけを素早く確認したいときに。
www.amazon.co.jp
回帰モデルの残差のプロット：plotresコマンド
モデル図、残差の累積分布、予測値に対する残差、QQプロットの4枚を既定で描画します。whichで描く図を選び、info=TRUEで学習済み決定係数などの補足情報を重ねられます。
| オプション | 意味 | 初期値 |
|---|---|---|
| object | モデルオブジェクト | stop(“no ‘object’ argument”) |
| which | 描画する図の種類、既定の1:4はモデル図・累積分布・予測値に対する残差・QQプロット | 1:4 |
| info | TRUEで残差分布や学習済み決定係数などの補足情報を表示 | FALSE |
| versus | 残差をプロットする対象、1で予測値、2で観測番号、3で応答、4でてこ比、文字列で予測変数 | 1 |
| standardize | TRUEで残差を標準化、一部のモデルのみ対応 | FALSE |
| delever | TRUEで残差をてこ比で補正、一部のモデルのみ対応 | FALSE |
| level | モデルが対応する場合に描く信頼区間・予測区間の水準、0で非表示 | 0 |
| id.n | ラベルを付ける残差の大きい順の個数、TRUEまたは-1で全件 | 3 |
| labels.id | 残差点に付けるラベル、id.nが正のときのみ使用 | NULL |
| smooth.col | 残差点を通す平滑線の色、0で非表示 | 2 |
| grid.col | 背景グリッドの色、0で非表示 | 0 |
| jitter | 残差点へ加える散らばりの大きさ、0で無効 | 0 |
| do.par | 作図パラメータの設定方法、NULL・FALSE・TRUE・2のいずれか | NULL |
| caption | 図全体の見出し、既定は自動生成、””で非表示 | NULL |
| trace | 実行時の情報表示の詳しさ | 0 |
| npoints | 描画する点の個数、大きい残差20点前後は必ず含む、TRUEまたは-1で全点 | 3000 |
| center | TRUEで残差プロットの横軸を中央に寄せ、分布の非対称を見やすくする | TRUE |
| type | residualsおよびpredictに渡すタイプ、既定は自動選択 | NULL |
| nresponse | residualsやpredictが複数列を返す場合に使用する列の番号または名前 | NA |
| object.name | エラーやトレース表示に用いるオブジェクト名 | quote.deparse(substitute(object)) |
| … | plotmoやpredictなどに渡すその他の引数 | なし |
# 残差の4種類の診断図を描画
plotres(lm_model)
# 学習済み決定係数などの補足情報を重ねて表示
plotres(lm_model, info = TRUE)
＜おすすめのRに関する書籍です＞
R言語ではじめるプログラミングとデータ分析 | 馬場 真哉
Rの文法から始めて、データ分析の基礎までひとつながりで学べる入門書。
勾配ブースティングモデルの学習曲線のプロット：plot_gbmコマンド
gbmコマンドで構築した勾配ブースティングモデルについて、木の本数に対する学習誤差・検証誤差・交差検証誤差・OOB誤差の曲線を描画します。既定では最小点付近を拡大して表示し、各曲線の最小点における木の本数を戻り値として返します。
| オプション | 意味 | 初期値 |
|---|---|---|
| object | gbmモデル | stop(“no ‘object’ argument”) |
| smooth | train・test・CV・OOBの各曲線を平滑化するかを表す4要素ベクトル | c(0, 0, 0, 1) |
| col | train・test・CV・OOBの各曲線の色を表す4要素ベクトル、0でその曲線を非表示 | c(1, 2, 3, 4) |
| ylim | Y軸範囲、”auto”で最小点付近を詳しく表示、NULLで曲線の全範囲 | “auto” |
| legend.x | 凡例のX位置、既定は自動配置、NAで凡例なし | NULL |
| legend.y | 凡例のY位置 | NULL |
| legend.cex | 凡例の文字の大きさ | 0.8 |
| grid.col | 背景グリッドの色、既定はグリッドなし | NA |
| n.trees | plotres用、残差計算に使った木の本数を示す縦線のX位置 | NA |
| col.n.trees | plotres用、n.treesを示す縦線の色 | “darkgray” |
| … | plot関数に渡すその他の引数 | なし |
# gbmパッケージの読み込み
library("gbm")
Loaded gbm 2.3.1
This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
# 勾配ブースティングモデルを構築
gbm_model <- gbm(shukaku ~ ., data = chaen, distribution = "gaussian",
n.trees = 600, shrinkage = 0.05, interaction.depth = 3,
train.fraction = 0.8, verbose = FALSE)
# 学習曲線を描画し、最小点の木の本数を確認
plot_gbm(gbm_model)
＜おすすめのRに関する書籍です＞
R言語ではじめるプログラミングとデータ分析 | 馬場 真哉
Rの文法から始めて、データ分析の基礎までひとつながりで学べる入門書。
glmnetモデルの係数経路のプロット：plot_glmnetコマンド
glmnetコマンドで推定した正則化回帰について、正則化パラメータの変化に対する各係数の推移を描画します。plot.glmnetコマンドと異なり、既定の横軸は減少するlog lambdaで、右側に表示する変数名の数をlabelで調整します。
| オプション | 意味 | 初期値 |
|---|---|---|
| x | glmnetモデル | stop(“no ‘x’ argument”) |
| xvar | 横軸に取る量、”rlambda”で減少するlog lambda、”lambda”でlog lambda、”norm”で係数のL1ノルム、”dev”で説明された逸脱度の割合 | c(“rlambda”, “lambda”, “norm”, “dev”) |
| label | 図の右側に表示する変数名の数、FALSEで非表示、TRUEで全て | 10 |
| nresponse | 多変量応答モデルでどの応答を描画するか | NA |
| grid.col | 背景グリッドの色、既定はグリッドなし | NA |
| s | plotres用、残差計算に使ったlambda sを示す縦線のX位置 | NA |
| … | plot.glmnetに渡すその他の引数 | なし |
# glmnetパッケージの読み込み
library("glmnet")
Loading required package: Matrix
Loaded glmnet 5.0
# 因子を含む説明変数を数値行列へ変換
x <- model.matrix(~ kion + nissho + hyoko + hiryo + hinshu, data = chaen)[, -1]
# 目的変数のベクトル
y <- chaen$shukaku
# Lasso回帰の係数経路を推定
glmnet_model <- glmnet(x, y)
# 正則化パラメータに対する係数の推移を描画
plot_glmnet(glmnet_model)
この記事が誰かの役に立ちますように。
スポンサーリンク
