このパッケージはクロスバリデーションを用いて統計・機械学習モデルの予測能力の推定や、ハイパーパラメーターの調整を行うパッケージです。 現在のバージョンはまだ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., glm
、gbm
など)にしか対応していません。
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)を参照。 |
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
を使う(整数型の二項データを仮定)。TRUE
を使う。levels(response)[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になります)。
応答変数のクラスごとのサンプル数に偏りがある場合、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.options
にlist(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.options
にlist(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
関数の引数のうち、methods
、op.prev
、 control
を指定することができ、そのほかの値は無視されます。これで問題がないだろうと思っているのですが、もし他の引数を使う必要がある場合には教えてください。 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.models
のseed
引数に適当な値を指定します。 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
引数にグループ分けに用いるデータの観測数と同じ長さのベクトルデータを指定します。 グループの指定にはcharacter
、factor
、integer
、logical
を用いることができます。
以下の例では異なる種を予測したときの性能を評価しています。
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
の候補に1
と5
、n.minobsinnode
の候補に1
と10
を指定しています。 cv.models
はgrid
に指定された候補パラメーターの組み合わせ全てに対して、モデルの構築とクロスバリデーションによる予測性能評価を行います。
# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
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 <- 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
を使う際には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
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
)
glm
やgbm
のweights
のように、一部の関数は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
は以下の条件に当てはまるときにはオフセット項を補正した予測値を使って計算を行います。
formula = y ~ x + offset(param)
もしくはoffset = param
形式で指定されている。formula = y ~ x + offset(data$param)
やoffset = data$param
は非対応。モデルがリンク関数を用いる場合、cv.models
は以下の例のようにオフセットに使われるデータがresponseスケール、オフセットの指定がlinkスケールであることを仮定します。 事前に変換済みのデータを用いた場合には補正がうまく行かないのでご注意下さい(変換済みであることを検出できないので、警告メッセージを表示できません)。
cv <- cv.models(
glm(Petal.Length ~ Sepal.Width + offset(log(Petal.Width)), data = iris)
)
ranger::ranger
ranger::ranger
のpredict
メソッドはオフセット項の補正を行うのか行わないのかよくわかりません。 このため、現在cv.models
ではrangerモデルのオフセットを補正していません。
glmmML::glmmML
glmmML
パッケージのglmmML
関数はpredict
メソッドを持ちません。 このため、cv.models
は自前でpredict
メソッドを実装してありますが、オフセット項の補正方法がわからないため、このメソッドは現在オフセット項の補正に対応していません。
cv.models
はランダム効果のあるモデルはランダム効果を無視した予測値(population-level prediction)を用いて、予測性能を計算しています。
オブジェクトの構造がよくわからない
cv.models
オブジェクトはまぁいいとして、cv.best.models
オブジェクトとか、cv.result
オブジェクトとか、なんか混乱している気がします。 なので将来的に再設計・リファクタリングしてもうちょっと整理したいと思っています。
テストしてない
とりあえずテストコードをほとんど作ってないので、うまく動いてないことがある可能性があります。
AUCがいつも0.5以上になる
対象のモデルがランダムより当てはまりが悪い場合、つまり、あり/なしに対して反対の予測をする可能性が高い場合、pROC::coordsはデフォルトであり/なしを入れ替えて計算を行うので、AUCは0.5以下にはならず、MCCも0以下にはなりません。 多くの文献がAUCは0.5~1の間に収まると仮定しているので、これをどう扱うのが正しいのか、まだ自分の中で結論が出ていません。 ということで、とりあえず対応は保留してあります。 近い将来pROCパッケージををOptimalCutpointsパッケージに入れ替えようと思っているので、その頃にまた考えます。
OptimalCutpointsパッケージで同様の問題が起こるのかはわかりません。 何かわかった方は教えてください。
MSEの定義
論文によってはMSEの定義が MSE = sum((response - predict) ^ 2) / (sample size - number of parameter) になっていることがあります。どの定義を使ったらいいのか、そのうち調べます。
Diagnostic Likelihood Ratioを使ったベストモデルの選択
ベストモデル選択の時、Diagnostic Likelihood Ratioが高いモデルがよいモデルだと考えてモデルを選んでいますが、本当にそれでよいのか自信がありません。もし、間違いに気付いた人は教えてください。
group
引数によるユーザー定義グループによるクロスバリデーションのサポート。find.best.models
がcv.models
に指定した乱数の種子を用いて結果を固定するようにした。1-npv
などの指標が計算されなくなり、Diagnostic Likelihood Ratioが計算されるようになった。