English

はじめに

このパッケージはクロスバリデーションを用いて統計・機械学習モデルの予測能力の推定や、ハイパーパラメーターの調整を行うパッケージです。 現在のバージョンはまだAPIや仕様の設計段階なので、将来仕様が変わるかもしれません。 もし、バグにお気づきだったりご要望があったりしたらメール(アドレスは以下のコードをRにペースト:「rawToChar(as.raw(c(109, 111, 103, 64, 102, 102, 112, 114, 105, 46, 97, 102, 102, 114, 99, 46, 103, 111, 46, 106, 112)))」)もしくはGitHubにご連絡下さい。

インストール

パッケージをインストールするには、以下のコマンドをRのコンソールにコピー&ペーストしてください。

install.packages(
    c("R6", "model.adapter", "cv.models"), type = "source",
    repos = c(
        "http://florivory.net/R/repos", options()$repos
    )
)

基本的な使い方

クロスバリデーションによるモデルの予測能力評価

モデルの予測能力を評価するにはcv.models関数を使います。 一番基本的な使い方は以下の例のようにcv.models関数の中に統計・機械学習モデルの呼び出しを直接入れます。

library(cv.models)
data(iris)
cv <- cv.models(glm(Petal.Length ~ ., data = iris))
print(cv)
## Result of cross validation
## Cross validation metrics:
##          mse      rmse r.squared  spearman   kendall q.squared    sd.mse
## 1 0.07160118 0.2637418 0.9773273 0.9338794 0.8322722 0.9740298 0.0257532
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.04762617  0.008654745  0.05044482 0.07480678   0.01150527

モデルを作る際にはglm(iris$Petal.Length ~ iris$Petal.Width)ではなくglm(Petal.Length ~ Petal.Width, data = iris)のように、dataにデータを入れてモデル作成関数を呼び出してください。前者をどうしても使わないといけない場合があるなら対応するので教えてください。

また、モデル作成関数の呼び出しを含むcallオブジェクトを作成して、cv.models関数に渡すこともできます。

# substitute()関数を使う場合。
call.glm <- substitute(glm(Petal.Length ~ ., data = iris))
cv <- cv.models(call.glm)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared spearman  kendall q.squared     sd.mse
## 1 0.0720893 0.2642072 0.9786271  0.95254 0.855111 0.9758841 0.02458692
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.05037491  0.006491544  0.02646257 0.06361177  0.007255751
# call()関数を使う場合。
call.glm <- call("glm", Petal.Length ~ ., data = iris)
cv <- cv.models(call.glm)
print(cv)
## Result of cross validation
## Cross validation metrics:
##          mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.07225901 0.2609248 0.9764348 0.9429454 0.8348614 0.9741313 0.03779596
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.06812785   0.01431566  0.01652288 0.04044752   0.01500444

一部のモデルはモデルオブジェクトを使ってクロスバリデーションを行うこともできます。 ただし、オブジェクトがcallを保存している関数(e.g., glmgbmなど)にしか対応していません。

model <- glm(Petal.Length ~ ., data = iris)
cv <- cv.models(model)
print(cv)
## Result of cross validation
## Cross validation metrics:
##          mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.07144265 0.2634059 0.9786586 0.9402786 0.8399003 0.9752645 0.02485682
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.04784199  0.008179856  0.02563387 0.04283829   0.01095734

計算される指標の一覧

cv.modelsはモデルの種類(回帰モデル・識別モデル)を自動的に判定し、指標を計算します。 現在、cv.modelsが計算する指標は以下の通りです。

回帰モデル

指標名 列名 定義/説明
Mean squared error (MSE) “mse” \(MSE = mean((prediction - response) ^ 2)\)
Root mean squared error (RMSE) “rmse” \(RMSE = sqrt(mean((prediction - response) ^ 2))\)
R二乗 “r.squared” \(R ^ 2 = cor(prediction, response) ^ 2\)
参考のために計算していますが、相関係数の2乗は切片と傾きが再調整されるので使用をおすすめしません。予測値と実測値をプロットした点が、y = xの線上に乗るのが理想的なモデルですが、例えば予測値と実測値をプロットした点がy = -2x + 1の線に完全に乗っている場合でも(つまり予測がぜんぜんダメでも)相関係数の2乗は1になります。
Spearmanのρ “spearman” \(\rho = cor(prediction, response, method = "spearman")\)
Kendallのτ “kendall” \(\tau = cor(prediction, response, method = "kendall")\)
Q二乗 “q.squared” \(Q ^ 2 = 1 - \sum((prediction - response) ^ 2) / \sum((response - mean(response)) ^ 2)\),
Consonni et al. (2009)を参照。
現在のバージョンではデータの平均(\(mean(response)\))にデフォルトでは各foldの平均値を、aggregate.method = "join"(後述)を使ったときにはデータ全体の平均を使っています。一貫性がないので次のバージョンあたりで整理する予定です。

識別モデル

指標名 列名 定義/説明
最適な閾値 “threshold” 在/不在を分ける最適な閾値。デフォルトではYouden’Jで決定される。
Specificity “specificity” https://en.wikipedia.org/wiki/Sensitivity_and_specificity
Sensitivity “sensitivity”
Accuracy “accuracy”
True negative count (TN) “tn”
True positive count (TP) “tp”
False negative count (FN) “fn”
False positive count (FP) “fp”
Negative predictive value (NPV) “npv”
Positive predictive value (PPV) “ppv”
Diagnostic likelihood ratio (DRL)+ “dlr.positive” \(DRL.positive = Sensitivity / (1 - Specificity)\),
Lopez-Raton et al. (2014)を参照。
Diagnostic likelihood ratio (DRL)- “dlr.negative” \(DRL.negative = (1 - Sensitivity) / Specificity\),
Lopez-Raton et al. (2014)を参照。
Matthews correlation coefficient (MCC) “mcc”
Informedness “informedness” \(Informedness = Sensitivity + Specificity - 1\)
Markedness “markedness” \(Markedness = PPV + NPV - 1\)
Log-likelihood “loglik” \(Log_likelihood = \sum(log(p * y + (1 - p) * (1 - y)))\),
Lawson et al. (2014)を参照。
尤度に基づいたR二乗 “rsq.loglik” Lawson et al. (2014)を参照。
Cohen’s Kappa “kappa” Cohen (1960)を参照。

predictに引数が必要なモデルを使う

gbmのような一部のモデルはpredict関数による予測値の計算時にもパラメーターを指定する必要があります。 このようなモデルをcv.modelsで扱うには、以下の例のようにpredictに渡される引数をcv.modelsの引数に追加します。

library(gbm)
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.1110352 0.3225373 0.9664271 0.9175476 0.8206688 0.9592213 0.05756581
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.08822272   0.02499711   0.1067213  0.1406065   0.02188559

陽性として扱うクラスを制御する

cv.modelsはデフォルトでは以下のルールで陽性として扱うレベルを決定しています。

  1. 応答変数が数値型の時には1を使う(整数型の二項データを仮定)。
  2. 応答変数が論理型の時にはTRUEを使う。
  3. 応答変数が因子型の時には一番目のレベル(levels(response)[1])を使う。
  4. 応答変数が文字列型の時には一意な値の1番目を使う。

因子型や文字列型の変数を応答変数として扱うモデルの場合、positive.class引数を使うことで、どのレベルを陽性として扱うかを指定することができます。 妖精として扱われなかったクラスは全て陰性として扱われます。 以下の例ではirisデータを使い、versicolorを陽性、setosaとvirginicaを陰性として扱い、計算を行います。

cv <- cv.models(
    gbm(Species ~ ., data = iris, n.trees = 50, distribution = "multinomial"),
    n.cores = 1, n.trees = 50, positive.class = "versicolor"
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##   method threshold sensitivity specificity       ppv npv dlr.positive
## 1 Youden 0.6129907           1   0.9725758 0.9383333   1          Inf
##   dlr.negative  fp fn       auc tp  tn accuracy informedness markedness
## 1            0 0.3  0 0.9901717  5 9.7     0.98    0.9725758  0.9383333
##         mcc    loglik rsq.loglik     kappa sd.threshold sd.sensitivity
## 1 0.9547984 -2.467852  0.6411764 0.9514063    0.4612871              0
##   sd.specificity    sd.ppv sd.npv sd.dlr.positive sd.dlr.negative
## 1      0.0443321 0.1012423      0             NaN               0
##       sd.fp sd.fn     sd.auc    sd.tp    sd.tn sd.accuracy sd.informedness
## 1 0.4830459     0 0.01836965 1.563472 1.418136  0.03220306       0.0443321
##   sd.markedness     sd.mcc sd.loglik sd.rsq.loglik   sd.kappa
## 1     0.1012423 0.07331211  2.515344     0.2936008 0.07889443

性能評価指標の計算方法を制御する

cv.modelsは性能評価指標を計算するとき、デフォルトでは各foldごとに指標を計算し、それを平均します。 しかし、この方法はLeave-one-out Cross Validation (LOOCV)を行う時や在・不在データのバランスが偏っている時などにはうまく行かないことがあります。 このような問題をどう扱うのが正しいのか、現在情報収集中なのですが、とりあえずaggregate.method = "join"を指定することで、全体の予測値・応答変数のデータを結合した上で、性能評価指標を計算することができます。

cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    aggregate.method = "join", n.trees = 50
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared sd.mse
## 1 0.1169926 0.3420418 0.9627806 0.9412695 0.7954671 0.9622056     NA
##   sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1      NA           NA          NA         NA           NA

この場合、各指標のSDは計算されません(NAになります)。

Class Stratificationを使う

応答変数のクラスごとのサンプル数に偏りがある場合、foldの分割をランダムに行うと、あるクラスのサンプルが全てのfoldに含まれないことがあります。 このような場合、Class Stratification(応答変数の各クラスが各foldに均等に含まれるようにデータを分割する分割方法)を使うと、問題が解決する場合があります。

Class Stratificationを使うにはcv.modelsの引数でstratify = TRUEを指定します。

library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
# setosaが少ない偏ったデータを作成。
test.data <- subset(iris, Species != "setosa")
test.data <- rbind(test.data, head(subset(iris, Species == "setosa"), 10))

# クロスバリデーション。
cv <- cv.models(
    randomForest(Species ~ ., data = iris, n.trees = 50),
    aggregate.method = "join", stratify = TRUE, n.cores = 1
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##   method threshold sensitivity specificity ppv npv dlr.positive
## 1 Youden      0.95           1           1   1   1          Inf
##   dlr.negative fp fn auc tp  tn accuracy informedness markedness mcc
## 1            0  0  0   1 50 100        1            1          1   1
##       loglik rsq.loglik kappa sd.threshold sd.sensitivity sd.specificity
## 1 -0.1965458  0.9875976     1           NA             NA             NA
##   sd.ppv sd.npv sd.dlr.positive sd.dlr.negative sd.fp sd.fn sd.auc sd.tp
## 1     NA     NA              NA              NA    NA    NA     NA    NA
##   sd.tn sd.accuracy sd.informedness sd.markedness sd.mcc sd.loglik
## 1    NA          NA              NA            NA     NA        NA
##   sd.rsq.loglik sd.kappa
## 1            NA       NA

今のところ、Class Stratificationは識別モデルに対してだけ実装されています。 回帰モデルでもClass Stratificationが必要な人はご連絡ください。

閾値の計算方法を制御する

旧バージョンのcv.modelsは識別モデルの閾値の決定にpROCパッケージのcoords関数を使い、閾値を決定する指標にはYouden’s Jを使っていました。 Ver. 0.1.0からこれをOptimalCutpointsパッケージのoptimal.cutpoints関数に置き換えたため、Youden’s J以外の指標も閾値決定に使えるようになりました。 デフォルトでは旧バージョンと同じくYouden’s Jで閾値を決定しますが、cutpoint.options引数を設定することで他の指標も使うことができます。 cutpoint.optionsにはoptimal.cutpointsに渡す引数をリストに入れて指定します。 たとえば、Youden’s Jの代わりにSensitivityを最大化するように閾値を決定するには、以下の例のようにcutpoint.optionslist(methods = "MaxSe")を指定します。

cv <- cv.models(
    gbm(Species ~ ., data = iris, n.trees = 50, distribution = "multinomial"),
    n.cores = 1, n.trees = 50, positive.class = "setosa",
    cutpoint.options = list(methods = "MaxSe")
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##   method threshold sensitivity specificity ppv npv dlr.positive
## 1  MaxSe 0.9956124           1           1   1   1          Inf
##   dlr.negative fp fn auc tp tn accuracy informedness markedness mcc
## 1            0  0  0   1  5 10        1            1          1   1
##        loglik rsq.loglik kappa sd.threshold sd.sensitivity sd.specificity
## 1 -0.02259671  0.8983966     1  0.002632476              0              0
##   sd.ppv sd.npv sd.dlr.positive sd.dlr.negative sd.fp sd.fn sd.auc
## 1      0      0             NaN               0     0     0      0
##      sd.tp    sd.tn sd.accuracy sd.informedness sd.markedness sd.mcc
## 1 1.414214 1.414214           0               0             0      0
##    sd.loglik sd.rsq.loglik sd.kappa
## 1 0.00721783     0.0168212        0

optimal.cutpointsは同時に複数の閾値決定方法で閾値を計算することができます。 以下の例ではcutpoint.optionslist(methods = c("MaxSe", "MaxSp"))を指定し、Sensitivityを最大化するときの閾値、Specificityを最大化するときの閾値両方を計算しています(とはいえ、このモデルだと同じ閾値になるのですが・・・)。

cv <- cv.models(
    gbm(Species ~ ., data = iris, n.trees = 50, distribution = "multinomial"),
    n.cores = 1, n.trees = 50, positive.class = "setosa",
    cutpoint.options = list(methods = c("MaxSe", "MaxSp"))
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##   method threshold sensitivity specificity ppv npv dlr.positive
## 1  MaxSe 0.9959348           1           1   1   1          Inf
## 2  MaxSp 0.9959348           1           1   1   1          Inf
##   dlr.negative fp fn auc tp tn accuracy informedness markedness mcc
## 1            0  0  0   1  5 10        1            1          1   1
## 2            0  0  0   1  5 10        1            1          1   1
##       loglik rsq.loglik kappa sd.threshold sd.sensitivity sd.specificity
## 1 -0.0217214  0.8977191     1  0.001891336              0              0
## 2 -0.0217214  0.8977191     1  0.001891336              0              0
##   sd.ppv sd.npv sd.dlr.positive sd.dlr.negative sd.fp sd.fn sd.auc
## 1      0      0             NaN               0     0     0      0
## 2      0      0             NaN               0     0     0      0
##      sd.tp    sd.tn sd.accuracy sd.informedness sd.markedness sd.mcc
## 1 1.763834 1.763834           0               0             0      0
## 2 1.763834 1.763834           0               0             0      0
##     sd.loglik sd.rsq.loglik sd.kappa
## 1 0.006943235    0.01251464        0
## 2 0.006943235    0.01251464        0

cutpoint.optionsにはoptimal.cutpoints関数の引数のうち、methodsop.prevcontrolを指定することができ、そのほかの値は無視されます。これで問題がないだろうと思っているのですが、もし他の引数を使う必要がある場合には教えてください。 optimal.cutpointsの詳しい使い方はマニュアルや文献(Lopez-Raton et al. 2014)を見てください。

データ分割数の設定

クロスバリデーションのデータ分割数を変えるにはfoldsオプションを変更します。

cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50, folds = 5
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##          mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.09713131 0.3090713 0.9697169 0.9419013 0.8216823 0.9672297 0.02673631
##     sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.0448088  0.009822299  0.01831205  0.0383783    0.0103613

乱数を固定する

クロスバリデーションや一部の統計モデルの当てはめには乱数が使われいてるため、計算の度に結果が変わります。 これを固定したい場合にはcv.modelsseed引数に適当な値を指定します。 seedを指定した場合、並列計算のありなしにかかわらず、同じ結果を返します。

# 実行の度に値が変わる。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman  kendall q.squared     sd.mse
## 1 0.1077953 0.3213418 0.9701033 0.9318539 0.834346 0.9585987 0.04654643
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.07098322   0.01239665  0.03755285 0.06017613   0.01853208
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.1020297 0.3106249 0.9697572 0.9334639 0.8325937 0.9654339 0.05024967
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.07847039   0.01311519  0.05309581  0.0948489   0.01666312
# seedに値を指定すると、結果が固定される。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50, seed = 12345
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.1059846 0.3208402 0.9668684 0.9323122 0.8262999  0.961941 0.04231908
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.05817747   0.01486364  0.04670806 0.06900898   0.01723832
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50, seed = 12345
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.1059846 0.3208402 0.9668684 0.9323122 0.8262999  0.961941 0.04231908
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.05817747   0.01486364  0.04670806 0.06900898   0.01723832

グループを指定してクロスバリデーションを行う

cv.modelsはデフォルトではランダムにグループ分けを行い、モデルの性能評価を行いますが、空間自己相関などが存在し、データが独立ではない場合、ランダムなグループ分けはモデルの性能を過大評価することがあります(Roberts et al. 2017)。 このような場合、完全にランダムなデータ分割ではなく、グループを考慮したクロスバリデーションを行うことで、問題を軽減できる可能性があり(Roberts et al. 2017)、blockCVパッケージ(Valavi et al. 2018)でそのようなグループ作成を行えるようです。 外部で定義したグループを用いるにはgroup引数にグループ分けに用いるデータの観測数と同じ長さのベクトルデータを指定します。 グループの指定にはcharacterfactorintegerlogicalを用いることができます。

以下の例では異なる種を予測したときの性能を評価しています。

cv <- cv.models(
    randomForest(Petal.Length ~ ., data = iris, n.cores = 1, n.trees = 50),
    n.cores = 1, group = iris$Species
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##       mse     rmse r.squared  spearman   kendall q.squared   sd.mse
## 1 4.19397 1.749769 0.1980847 0.3585859 0.2723849 -106.9464 4.610268
##    sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 1.303234    0.3233384   0.3140936  0.2351445     177.1127

並列計算を制御する

cv.modelsはデフォルトで全てのCPUコアを使用して計算を行います。 計算に使われるコア数を制御したいときには、n.cores引数に適当な値を指定します。

# コアを2つ使って計算する。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50, n.cores = 2
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##         mse      rmse r.squared  spearman   kendall q.squared     sd.mse
## 1 0.1012589 0.3099502 0.9712767 0.9448699 0.8547116 0.9650142 0.04881779
##      sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.07593647   0.01071404  0.03180521 0.06488336   0.01608696

モデルのパラメーターを選択する

cv.modelsを使うとクロスバリデーションによる予測性能を用いて、モデルのハイパーパラメーターの選択を行うこともできます。 ハイパーパラメーターの選択を行うにはgrid引数に、候補となるパラメーターのベクトルを格納したリストを指定します。 以下の例ではgbmのハイパーパラメーターのinteraction.depthの候補に15n.minobsinnodeの候補に110を指定しています。 cv.modelsgridに指定された候補パラメーターの組み合わせ全てに対して、モデルの構築とクロスバリデーションによる予測性能評価を行います。

# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50,
    grid = list(interaction.depth = c(1, 5), n.minobsinnode = c(1, 10))
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##   interaction.depth n.minobsinnode        mse      rmse r.squared
## 1                 1              1 0.07988630 0.2798786 0.9742615
## 2                 5              1 0.06872487 0.2577735 0.9772989
## 3                 1             10 0.10508683 0.3198442 0.9680260
## 4                 5             10 0.09860740 0.3079757 0.9666660
##    spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1 0.9359951 0.8447668 0.9708169 0.02398398 0.04155642   0.01329978
## 2 0.9429622 0.8400153 0.9738134 0.02762220 0.05030685   0.01088985
## 3 0.9164073 0.8077601 0.9626880 0.03695090 0.05564301   0.01371834
## 4 0.9083837 0.7899425 0.9636894 0.04147805 0.06462149   0.01675626
##   sd.spearman sd.kendall sd.q.squared
## 1  0.03642433 0.06185616   0.01265028
## 2  0.03219425 0.06936421   0.01260916
## 3  0.09274789 0.11435953   0.01490257
## 4  0.09086103 0.11045777   0.01662154

また、gbmのように、predict関数の引数にも調整可能なパラメーターがあるモデルの場合、grid.predictに同様のリストを指定し、パラメーターの違いによる予測性能の違いを評価することができます。

# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    grid = list(interaction.depth = c(1, 5), n.minobsinnode = c(1, 10)),
    grid.predict = list(n.trees = c(10, 50, 80))
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##    interaction.depth n.minobsinnode n.trees        mse      rmse r.squared
## 1                  1              1      10 0.72135299 0.8408298 0.9224749
## 2                  1              1      50 0.07437223 0.2641699 0.9766791
## 3                  1              1      80 0.06657522 0.2497271 0.9779838
## 4                  5              1      10 0.46146821 0.6727457 0.9752794
## 5                  5              1      50 0.07149702 0.2639825 0.9776135
## 6                  5              1      80 0.07833486 0.2764058 0.9756564
## 7                  1             10      10 0.71841345 0.8328385 0.9225018
## 8                  1             10      50 0.10433898 0.3153330 0.9661261
## 9                  1             10      80 0.09594713 0.3013134 0.9675966
## 10                 5             10      10 0.49663532 0.6985129 0.9686535
## 11                 5             10      50 0.09840113 0.3053468 0.9719476
## 12                 5             10      80 0.09214228 0.2968738 0.9738044
##     spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1  0.9100195 0.8139757 0.7533619 0.22153823 0.12630736  0.031931726
## 2  0.9464613 0.8578524 0.9738740 0.03624250 0.07138713  0.012497343
## 3  0.9517171 0.8659651 0.9764205 0.03415546 0.06840736  0.012573607
## 4  0.9394519 0.8575398 0.8394801 0.13002987 0.09933937  0.016536960
## 5  0.9315475 0.8300317 0.9741809 0.02282192 0.04484861  0.008948665
## 6  0.9346320 0.8337089 0.9718290 0.02375459 0.04636464  0.009063234
## 7  0.9058810 0.8109394 0.7365153 0.29370950 0.16597668  0.043481842
## 8  0.9108231 0.8202898 0.9589747 0.04656443 0.07381725  0.023443762
## 9  0.8996299 0.7966864 0.9614480 0.04429300 0.07569933  0.024146913
## 10 0.9464612 0.8489458 0.8347409 0.14019311 0.09840448  0.019054120
## 11 0.9357587 0.8244511 0.9659731 0.04899660 0.07575139  0.015987939
## 12 0.9303874 0.8186353 0.9680151 0.04114785 0.06673509  0.014146147
##    sd.spearman sd.kendall sd.q.squared
## 1   0.03140692 0.04132981  0.061570485
## 2   0.03006538 0.05517784  0.013313957
## 3   0.02801302 0.04286534  0.013154050
## 4   0.03547433 0.05082152  0.023407873
## 5   0.04800136 0.07469330  0.009038385
## 6   0.04301943 0.07100622  0.009030727
## 7   0.05438383 0.06365259  0.083130132
## 8   0.09329241 0.13637365  0.024444145
## 9   0.09960391 0.15109024  0.025276808
## 10  0.02209911 0.04808539  0.036860816
## 11  0.02818431 0.04750567  0.018693113
## 12  0.03442903 0.05676997  0.016835996

最良モデルの取り出し

find.best.models関数を用いることで、cv.modelsの結果から最良モデルを取り出すことができます。 以下の例ではQ^2が最高になるモデルを選んでいます。

# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    grid = list(interaction.depth = c(1, 5), n.minobsinnode = c(1, 10)),
    grid.predict = list(n.trees = c(10, 50, 80))
)
print(cv)
## Result of cross validation
## Cross validation metrics:
##    interaction.depth n.minobsinnode n.trees        mse      rmse r.squared
## 1                  1              1      10 0.71167534 0.8270316 0.9374110
## 2                  1              1      50 0.07633178 0.2728128 0.9772633
## 3                  1              1      80 0.06970996 0.2611661 0.9776161
## 4                  5              1      10 0.45281760 0.6701436 0.9750218
## 5                  5              1      50 0.07109058 0.2625991 0.9739245
## 6                  5              1      80 0.07727443 0.2733543 0.9719832
## 7                  1             10      10 0.71146347 0.8261406 0.9370324
## 8                  1             10      50 0.10740502 0.3166357 0.9689644
## 9                  1             10      80 0.09877502 0.3048053 0.9706678
## 10                 5             10      10 0.49193687 0.6941843 0.9711012
## 11                 5             10      50 0.09008986 0.2964273 0.9736137
## 12                 5             10      80 0.08848874 0.2940106 0.9741172
##     spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1  0.9140896 0.8194126 0.7644424 0.31880848 0.17541726  0.022041155
## 2  0.9418795 0.8499907 0.9735839 0.02678038 0.04600656  0.008559696
## 3  0.9366348 0.8371464 0.9750647 0.02200164 0.04085541  0.007939684
## 4  0.9510719 0.8660040 0.8337415 0.09028281 0.06433510  0.013696117
## 5  0.9432535 0.8464335 0.9716884 0.02636445 0.04867451  0.019067865
## 6  0.9454451 0.8497483 0.9693665 0.03025065 0.05324833  0.020593886
## 7  0.9190058 0.8274564 0.7627697 0.31277873 0.17936674  0.033327878
## 8  0.9112350 0.7899700 0.9637852 0.05728259 0.08911181  0.013159216
## 9  0.9135241 0.7877474 0.9664255 0.04853941 0.08075179  0.011550374
## 10 0.9267630 0.8192063 0.8377699 0.14277131 0.10564650  0.010978775
## 11 0.9438938 0.8333159 0.9692335 0.02950186 0.04967337  0.008040985
## 12 0.9422706 0.8342727 0.9696525 0.02822796 0.04768542  0.008342816
##    sd.spearman sd.kendall sd.q.squared
## 1   0.04705451 0.05682279  0.047419462
## 2   0.02446964 0.04300466  0.007060169
## 3   0.02530149 0.04402734  0.008251505
## 4   0.01492247 0.03505610  0.057982144
## 5   0.03209873 0.05947417  0.019892648
## 6   0.02695628 0.05381158  0.021819174
## 7   0.05046289 0.06740698  0.053868218
## 8   0.05373625 0.06788535  0.015249590
## 9   0.04972679 0.06481529  0.013258042
## 10  0.04770969 0.08144966  0.026125456
## 11  0.02586059 0.05719960  0.009364542
## 12  0.02284420 0.04468042  0.009575278
# 最良モデルの取り出し。
best <- find.best.models(cv, "q.squared")
print(best)
## ------------------------------------------------
## Result of model selection by cross validation
## -> 1 best models was found:
## ------------------------------------------------
## Model index 1
## ------------------------------------------------ 
## Result of cross validation
## Cross validation metrics:
##   interaction.depth n.minobsinnode n.trees        mse      rmse r.squared
## 1                 1              1      80 0.06970996 0.2611661 0.9776161
##    spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1 0.9366348 0.8371464 0.9750647 0.02200164 0.04085541  0.007939684
##   sd.spearman sd.kendall sd.q.squared
## 1  0.02530149 0.04402734  0.008251505

find.best.models関数の結果はcv.best.modelsオブジェクトになります。 cv.best.modelsの実体はcv.resultオブジェクトのリストです。

# find.best.models()の結果のクラスを表示。
class(best)
## [1] "cv.best.models"
# find.best.models()の結果はcv.resultクラスのリスト。
class(best[[1]])
## [1] "cv.result"
print(best[[1]])
## Result of cross validation
## Cross validation metrics:
##   interaction.depth n.minobsinnode n.trees        mse      rmse r.squared
## 1                 1              1      80 0.06970996 0.2611661 0.9776161
##    spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1 0.9366348 0.8371464 0.9750647 0.02200164 0.04085541  0.007939684
##   sd.spearman sd.kendall sd.q.squared
## 1  0.02530149 0.04402734  0.008251505

cv.modelsオブジェクトからのデータの取り出し

以下の関数を使ってcv.modelsオブジェクトからデータを取り出すことができます。

cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    grid = list(interaction.depth = c(1, 5), n.minobsinnode = c(1, 10)),
    grid.predict = list(n.trees = c(10, 50, 80))
)

# 10番目のモデルの予測結果を取得。
# 結果はdata.frameで、response列は応答変数の値(生データ)、
# prediction列はモデルの予測値、index列はモデリングに使われたデータでの列番号。
fit <- extract.fit(cv, 10)
head(fit)
##   response prediction index
## 1      1.4   2.276516     1
## 2      1.4   2.255164     2
## 3      1.6   2.249139    30
## 4      1.4   2.276516    34
## 5      1.3   2.249139    43
## 6      4.9   4.449403    53
# 10番目のモデルの詳細情報を取得。
# 結果はcv.resultオブジェクト。
extract.result(cv, 10)
## Result of cross validation
## Cross validation metrics:
##   interaction.depth n.minobsinnode n.trees       mse      rmse r.squared
## 1                 5             10      10 0.4936772 0.6986842 0.9696327
##    spearman   kendall q.squared   sd.mse    sd.rmse sd.r.squared
## 1 0.9432281 0.8429034 0.8370312 0.110278 0.07829823   0.02237619
##   sd.spearman sd.kendall sd.q.squared
## 1  0.01946093 0.03992628   0.02108426
# 性能評価の表を取得
extract.metrics(cv)
##    interaction.depth n.minobsinnode n.trees        mse      rmse r.squared
## 1                  1              1      10 0.71055597 0.8267423 0.9374964
## 2                  1              1      50 0.08046765 0.2760962 0.9784017
## 3                  1              1      80 0.07417160 0.2652973 0.9793937
## 4                  5              1      10 0.45216030 0.6659131 0.9811965
## 5                  5              1      50 0.06951450 0.2617910 0.9803264
## 6                  5              1      80 0.07096561 0.2645191 0.9791059
## 7                  1             10      10 0.73633947 0.8437043 0.9080399
## 8                  1             10      50 0.10587326 0.3163376 0.9572463
## 9                  1             10      80 0.09896640 0.3065534 0.9582283
## 10                 5             10      10 0.49367720 0.6986842 0.9696327
## 11                 5             10      50 0.09058819 0.2875990 0.9723143
## 12                 5             10      80 0.08860219 0.2847839 0.9729779
##     spearman   kendall q.squared     sd.mse    sd.rmse sd.r.squared
## 1  0.9305355 0.8443569 0.7668644 0.27800977 0.17337557  0.023233970
## 2  0.9367463 0.8460317 0.9719726 0.03741907 0.06862551  0.011128382
## 3  0.9371733 0.8433806 0.9737136 0.03381360 0.06488393  0.010702370
## 4  0.9544721 0.8748789 0.8500226 0.13932830 0.09843247  0.006896546
## 5  0.9470792 0.8508433 0.9764101 0.01708627 0.03299778  0.005994176
## 6  0.9343438 0.8276356 0.9756209 0.01729721 0.03325448  0.007359089
## 7  0.8969229 0.8032909 0.6810873 0.27860921 0.16500002  0.117739581
## 8  0.9228504 0.8188642 0.9509120 0.05006332 0.08030339  0.041988400
## 9  0.9142913 0.8008107 0.9537898 0.04454427 0.07447134  0.042144557
## 10 0.9432281 0.8429034 0.8370312 0.11027801 0.07829823  0.022376187
## 11 0.9395509 0.8267150 0.9697106 0.05674628 0.09354147  0.018542407
## 12 0.9415195 0.8350511 0.9704691 0.05453334 0.09128904  0.017661454
##    sd.spearman sd.kendall sd.q.squared
## 1   0.02856329 0.04808475  0.042827590
## 2   0.02898569 0.05839199  0.013520327
## 3   0.03255064 0.06756512  0.013780382
## 4   0.02399544 0.04190647  0.027951917
## 5   0.02248642 0.04457215  0.005891470
## 6   0.03905071 0.06728755  0.007310669
## 7   0.05940689 0.06239035  0.298277308
## 8   0.04986328 0.06419038  0.050315486
## 9   0.06817481 0.09408433  0.046649363
## 10  0.01946093 0.03992628  0.021084265
## 11  0.02016400 0.04580704  0.018840181
## 12  0.01767506 0.03643516  0.017926971

extract.fit関数はcv.best.modelsオブジェクト、cv.resultオブジェクトにも使えます。

# ベストモデルの1番目のモデルから、予測結果を取得。
best <- find.best.models(cv, "q.squared")
fit <- extract.fit(best, 1)
head(fit)
##   response prediction index
## 1      1.3   1.500114     3
## 2      1.5   1.525980    40
## 3      1.5   1.465115    49
## 4      4.5   4.223293    67
## 5      4.1   3.656145    68
## 6      4.8   4.611026    71
# 10番目のモデルをcv.resultオブジェクトとして取り出し。
result <- extract.result(cv, 10)
# 10番目のモデルの予測結果を取得。
fit <- extract.fit(result)
head(fit)
##   response prediction index
## 1      1.4   2.276516     1
## 2      1.4   2.255164     2
## 3      1.6   2.249139    30
## 4      1.4   2.276516    34
## 5      1.3   2.249139    43
## 6      4.9   4.449403    53

簡易作図

plot関数を使って、予測値と応答変数の関係を作図することができます。 1点が生データの1点を表し、線は\(Y = X\)の線を表します。 ただし、識別モデルに対してはうまく動きません。

# lmモデルの予測値と応答変数の関係をプロット。
cv <- cv.models(
    lm(Petal.Length ~ ., data = iris)
)
plot(cv)

# gbmモデルのパラメーター組み合わせ候補を作成。
cv <- cv.models(
    gbm(
        Petal.Length ~ ., data = iris, weights = iris$Sepal.Width,
        distribution = "gaussian", n.cores = 1
    ),
    grid = list(interaction.depth = c(1, 5), n.minobsinnode = c(1, 10)),
    grid.predict = list(n.trees = c(10, 50, 80))
)
# 10番目の予測結果を作図。
plot(cv, 10)

技術的詳細

@準備中

モデルごとの注意事項

gbm

n.treesの指定

gbmを使う際にはgrid.predictに渡すn.treesの候補で使われる以上の値をgbmの呼び出しにも指定して下さい。 以下の例ではgrid.predictに指定されたn.treesの候補のなかで、gbm呼び出しのn.treesの値を超えたものは、どれもgbmの呼び出しに指定されたn.treesの値を使ったときと同じ結果になってしまっています。

# gbm呼び出しのn.treesは10だが、
# grid.predictではn.treesの候補として5、10、50、80を指定している。
cv <- cv.models(
    gbm(
        Petal.Length ~ ., data = iris, weights = iris$Sepal.Width,
        distribution = "gaussian", n.trees = 10
    ),
    grid.predict = list(n.trees = c(5, 10, 50, 80))
)
# そのような場合、gbmの呼び出しに使われたn.treesの値(10)以上の値を
# grid.predictに指定しても、予測性能はn.trees = 10の時と同じになる。
print(cv)
## Result of cross validation
## Cross validation metrics:
##   n.trees       mse      rmse r.squared  spearman   kendall q.squared
## 1       5 1.3757917 1.1672744 0.8703568 0.8406009 0.7338241 0.5324178
## 2      10 0.6899316 0.8230134 0.9254704 0.9177348 0.8269810 0.7669629
## 3      50 0.6899316 0.8230134 0.9254704 0.9177348 0.8269810 0.7669629
## 4      80 0.6899316 0.8230134 0.9254704 0.9177348 0.8269810 0.7669629
##      sd.mse   sd.rmse sd.r.squared sd.spearman sd.kendall sd.q.squared
## 1 0.2910870 0.1213905   0.10703870  0.06042366 0.06285473   0.06734191
## 2 0.1953957 0.1182299   0.05092181  0.04161940 0.05705307   0.05109554
## 3 0.1953957 0.1182299   0.05092181  0.04161940 0.05705307   0.05109554
## 4 0.1953957 0.1182299   0.05092181  0.04161940 0.05705307   0.05109554

weightsの指定

formulaを関数の呼び出し外で定義し、以下のような方法でweightsを指定した場合、gbmに対する性能評価がうまく動きません。

# これは動く。
cv <- cv.models(
    gbm(
        Petal.Length ~ ., data = iris, weights = iris$Sepal.Width,
        distribution = "gaussian", n.trees = 10
    ),
    n.trees = 10, n.cores = 1
)

# 原因はわからないが、これは動かない。
f <- Petal.Length ~ .
cv <- cv.models(
    gbm(
        f, data = iris, weights = iris$Sepal.Width,
        distribution = "gaussian", n.trees = 10
    ),
    n.trees = 10, n.cores = 1
)

原因は不明ですが、以下の例のようにiris$を除去してweightsを指定すれば動作させることが可能です。

f <- Petal.Length ~ .
cv <- cv.models(
    gbm(
        f, data = iris, weights = Sepal.Width,
        distribution = "gaussian", n.trees = 10
    ),
    n.trees = 10, n.cores = 1
)

data以外にもデータを用いる関数(gbmやglmのweightなど)

glmgbmweightsのように、一部の関数はdata以外にもデータを用いて推定を行います。 このような場合、call関数で作成したcallオブジェクトはcv.modelsに用いることができません。

# weightsを指定し、かつcall関数で作成したcallオブジェクトは
# cv.modelsに使えない。
call.gbm <- call(
    "gbm", Petal.Length ~ ., data = iris, weights = iris$Sepal.Width,
    distribution = "gaussian", n.trees = 10, n.cores = 1
)
# エラー
cv <- cv.models(call.gbm, n.trees = 10)

# GLMも同様。
call.glm <- call(
    "glm", Petal.Length ~ ., data = iris, weights = iris$Sepal.Width
)
cv <- cv.models(call.glm)

このような場合にはcall関数ではなく、substitute関数を使い、また、weightsの指定にはiris$を使わず、変数名を指定して下さい。

# substituteでcallを作成。
call.gbm <- substitute(
    gbm(
        Petal.Length ~ ., data = iris, weights = Sepal.Width,
        distribution = "gaussian", n.trees = 10, n.cores = 1
    )
)
# OK!
cv <- cv.models(call.gbm, n.trees = 10)

# GLMも同様。
call.glm <- substitute(
    glm(Petal.Length ~ ., data = iris, weights = Sepal.Width)
)
cv <- cv.models(call.glm)

オフセット項のあるモデル

一部のモデルのpredictメソッドははオフセット項があるときに予測値を正しく計算しませんが、cv.modelsは以下の条件に当てはまるときにはオフセット項を補正した予測値を使って計算を行います。

  • モデルが回帰モデル。
  • オフセット項は1つ。
  • オフセット項がformula = y ~ x + offset(param)もしくはoffset = param形式で指定されている。formula = y ~ x + offset(data$param)offset = data$paramは非対応。
  • 予測値をresponseスケールで計算。(そういう状況があるのかわかりませんが)linkスケールでは補正が行われません。
  • モデルがglmmMLとrangerではない(以下参照)。

モデルがリンク関数を用いる場合、cv.modelsは以下の例のようにオフセットに使われるデータがresponseスケール、オフセットの指定がlinkスケールであることを仮定します。 事前に変換済みのデータを用いた場合には補正がうまく行かないのでご注意下さい(変換済みであることを検出できないので、警告メッセージを表示できません)。

cv <- cv.models(
    glm(Petal.Length ~ Sepal.Width + offset(log(Petal.Width)), data = iris)
)

ranger::ranger

ranger::rangerpredictメソッドはオフセット項の補正を行うのか行わないのかよくわかりません。 このため、現在cv.modelsではrangerモデルのオフセットを補正していません。

glmmML::glmmML

glmmMLパッケージのglmmML関数はpredictメソッドを持ちません。 このため、cv.modelsは自前でpredictメソッドを実装してありますが、オフセット項の補正方法がわからないため、このメソッドは現在オフセット項の補正に対応していません。

ランダム効果のあるモデル

cv.modelsはランダム効果のあるモデルはランダム効果を無視した予測値(population-level prediction)を用いて、予測性能を計算しています。

既知の問題

引用文献

更新履歴