In Marchen/cv.models: Model Selection and Hyper-Parameter Tuning by Cross Validation

```{R, preparation, echo = FALSE, message = FALSE, include = FALSE} library(cv.models) library(gbm) set.seed(12345)

# はじめに

このパッケージはクロスバリデーションを用いて統計・機械学習モデルの予測能力の推定や、ハイパーパラメーターの調整を行うパッケージです。
現在のバージョンはまだ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](https://github.com/Marchen/cv.models)にご連絡下さい。

# インストール

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

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

基本的な使い方

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

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

```{R, basic_usage} library(cv.models) data(iris) cv <- cv.models(glm(Petal.Length ~ ., data = iris)) print(cv)

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

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

```{R, basic_usage2}
# substitute()関数を使う場合。
call.glm <- substitute(glm(Petal.Length ~ ., data = iris))
cv <- cv.models(call.glm)
print(cv)

# call()関数を使う場合。
call.glm <- call("glm", Petal.Length ~ ., data = iris)
cv <- cv.models(call.glm)
print(cv)

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

```{R, wrong_usage} model <- glm(Petal.Length ~ ., data = iris) cv <- cv.models(model) print(cv)

## 計算される指標の一覧

`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$<br>参考のために計算していますが、相関係数の２乗は切片と傾きが再調整されるので使用をおすすめしません。予測値と実測値をプロットした点が、y = xの線上に乗るのが理想的なモデルですが、例えば予測値と実測値をプロットした点がy = -2x + 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)$, <br>Consonni et al. (2009)を参照。<br>現在のバージョンではデータの平均（$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)$, <br>Lopez-Raton et al. (2014)を参照。
**Diagnostic likelihood ratio (DRL)-**    |"dlr.negative"  |$DRL.negative = (1 - Sensitivity) / Specificity$, <br>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)))$, <br>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`の引数に追加します。

```{R, predict_args}
library(gbm)
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50
)
print(cv)

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

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

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

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

```{R, positive_class} cv <- cv.models( gbm(Species ~ ., data = iris, n.trees = 50, distribution = "multinomial"), n.cores = 1, n.trees = 50, positive.class = "versicolor" ) print(cv)

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

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

```{R, aggregate_method}
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    aggregate.method = "join", n.trees = 50
)
print(cv)

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

Class Stratificationを使う

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

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

```{R, class_stratification} library(randomForest)

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)

今のところ、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")`を指定します。

```{R, cutoff_example1}
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)

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

```{R, cutoff_example2} 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)

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

## データ分割数の設定

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

```{R, folds}
cv <- cv.models(
    gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1),
    n.trees = 50, folds = 5
)
print(cv)

乱数を固定する

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

```{R, set_seed}

実行の度に値が変わる。

cv <- cv.models( gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1), n.trees = 50 ) print(cv) cv <- cv.models( gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1), n.trees = 50 ) print(cv)

seedに値を指定すると、結果が固定される。

cv <- cv.models( gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1), n.trees = 50, seed = 12345 ) print(cv) cv <- cv.models( gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1), n.trees = 50, seed = 12345 ) print(cv)

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

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

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

```{R, grouped_cross_validation}
cv <- cv.models(
    randomForest(Petal.Length ~ ., data = iris, n.cores = 1, n.trees = 50),
    n.cores = 1, group = iris$Species
)
print(cv)

並列計算を制御する

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

```{R, n_cores}

コアを２つ使って計算する。

cv <- cv.models( gbm(Petal.Length ~ ., data = iris, distribution = "gaussian", n.cores = 1), n.trees = 50, n.cores = 2 ) print(cv)

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

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

```{R, hyperparameter1}
# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
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)

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

```{R, hyperparmeter2}

ハイパーパラメータが違うと予測性能がどう変わるかを評価する。

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)

## 最良モデルの取り出し

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

```{R, bestmodel}
# ハイパーパラメータが違うと予測性能がどう変わるかを評価する。
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)
# 最良モデルの取り出し。
best <- find.best.models(cv, "q.squared")
print(best)

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

```{R, bestmodel2}

find.best.models()の結果のクラスを表示。

class(best)

find.best.models()の結果はcv.resultクラスのリスト。

class(best[[1]]) print(best[[1]])

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

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

```{R, extract_data}
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)

# 10番目のモデルの詳細情報を取得。
# 結果はcv.resultオブジェクト。
extract.result(cv, 10)

# 性能評価の表を取得
extract.metrics(cv)

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

```{R, extract_data2}

ベストモデルの１番目のモデルから、予測結果を取得。

best <- find.best.models(cv, "q.squared") fit <- extract.fit(best, 1) head(fit)

10番目のモデルをcv.resultオブジェクトとして取り出し。

result <- extract.result(cv, 10)

10番目のモデルの予測結果を取得。

fit <- extract.fit(result) head(fit)

## 簡易作図

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

```{R, plot}
# 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の値を使ったときと同じ結果になってしまっています。

```{R, gbm_wrong_ntree}

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)

### weightsの指定

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

```{R, gbm_weights_1, eval = FALSE}
# これは動く。
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を指定すれば動作させることが可能です。

```{R, gbm_weights_2} 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など）

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

```{R, call_with_other_data_error, eval = FALSE}
# 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$を使わず、変数名を指定して下さい。

```{R, call_with_other_data_ok}

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`は非対応。
* 予測値をresponseスケールで計算。（そういう状況があるのかわかりませんが）linkスケールでは補正が行われません。
* モデルがglmmMLとrangerではない（以下参照）。

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

```{R}
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が高いモデルがよいモデルだと考えてモデルを選んでいますが、本当にそれでよいのか自信がありません。もし、間違いに気付いた人は教えてください。

引用文献

Cohen, J. 1960. A Coefficient of Agreement for Nominal Scales. Educational and Psychological Measurement 20: 37-46.
Consonni, V., D. Ballabio, and R. Todeschini. 2009. Comments on the Definition of the Q2 Parameter for QSAR Validation. Journal of Chemical Information and Modeling 49:1669-1678.
Lawson CR, Hodgson JA, Wilson RJ, et al. 2014. Prevalence, thresholds and the performance of presence-absence models. Methods in Ecology and Evolution 5:54-64.
Lopez-Raton, M., M. X. Rodriguez-Alvarez, C. Cadarso-Suarez, and F. Gude-Sampedro. 2014. OptimalCutpoints: An R Package for Selecting Optimal Cutpoints in Diagnostic Tests. Journal of Statistical Software 61:36.
Roberts, D. R., V. Bahn, S. Ciuti, M. S. Boyce, J. Elith, G. Guillera-Arroita, S. Hauenstein, J. J. Lahoz-Monfort, B. Schröder, W. Thuiller, D. I. Warton, B. A. Wintle, F. Hartig, and C. F. Dormann. 2017. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography 40:913-929.
Valavi, R., J. Elith, J. J. Lahoz-Monfort, and G. Guillera-Arroita. 2018. blockCV: an R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods in Ecology and Evolution 00:1-8.

更新履歴

2019-07-09: 0.1.4
回帰モデルの性能指標にSpearman's ρとKendall's τを追加。
Travis-CIによる継続的インテグレーションの開始。
バグ修正。
2018-11-10: 0.1.3
group引数によるユーザー定義グループによるクロスバリデーションのサポート。
2018-09-14: 0.1.2
model.adapterパッケージver.0.1.0に対応。
Class Stratificationに対応。
2018-01-31: 0.1.1
call()関数やモデルオブジェクトを使ってcv.modelsを呼び出した場合、クロスバリデーションが行われないバグを修正した。
2018-01-28: 0.1.0
find.best.modelsがcv.modelsに指定した乱数の種子を用いて結果を固定するようにした。
pROCパッケージへの依存をOptimalCutpointsパッケージに置き換え、Youden's J以外の指標でも閾値を選択できるようにした。この変更により1-npvなどの指標が計算されなくなり、Diagnostic Likelihood Ratioが計算されるようになった。
Cohen's Kappaを識別モデルの性能評価指標に追加した。
バグ修正：seedを指定したとき、並列計算ありなしで結果が変わるバグを修正。
バグ修正：print.cv.best.models()のメッセージの修正。
2017-09-10: 0.0.4
試験的に対数尤度と尤度に基づいたR^2を識別モデルの性能評価指標に追加。
バグ修正：CPUコアの割り当てがおかしかったバグを修正。
2017-07-07: 0.0.3
バグ修正：まれに計算に失敗する問題を修正。
2017-07-07: 0.0.2
バグ修正：説明変数と応答変数の対応がおかしくなっていたバグを修正。
2017-07-05: 0.0.1
バージョン付け開始。

Marchen/cv.models documentation built on Sept. 2, 2020, 2:04 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Marchen/cv.models
Model Selection and Hyper-Parameter Tuning by Cross Validation

In Marchen/cv.models: Model Selection and Hyper-Parameter Tuning by Cross Validation

基本的な使い方

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

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

Class Stratificationを使う

setosaが少ない偏ったデータを作成。

クロスバリデーション。

乱数を固定する

実行の度に値が変わる。

seedに値を指定すると、結果が固定される。

並列計算を制御する

コアを２つ使って計算する。

ハイパーパラメータが違うと予測性能がどう変わるかを評価する。

find.best.models()の結果のクラスを表示。

find.best.models()の結果はcv.resultクラスのリスト。

ベストモデルの１番目のモデルから、予測結果を取得。

10番目のモデルをcv.resultオブジェクトとして取り出し。

10番目のモデルの予測結果を取得。

技術的詳細

モデルごとの注意事項

gbm

n.treesの指定

gbm呼び出しのn.treesは10だが、

grid.predictではn.treesの候補として5、10、50、80を指定している。

そのような場合、gbmの呼び出しに使われたn.treesの値（10）以上の値を

grid.predictに指定しても、予測性能はn.trees = 10の時と同じになる。

substituteでcallを作成。

OK!

GLMも同様。

`ranger::ranger`

`glmmML::glmmML`

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

既知の問題

引用文献

更新履歴

R Package Documentation

Browse R Packages

We want your feedback!

Marchen/cv.models Model Selection and Hyper-Parameter Tuning by Cross Validation

In Marchen/cv.models: Model Selection and Hyper-Parameter Tuning by Cross Validation

基本的な使い方

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

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

Class Stratificationを使う

setosaが少ない偏ったデータを作成。

クロスバリデーション。

乱数を固定する

実行の度に値が変わる。

seedに値を指定すると、結果が固定される。

並列計算を制御する

コアを２つ使って計算する。

ハイパーパラメータが違うと予測性能がどう変わるかを評価する。

find.best.models()の結果のクラスを表示。

find.best.models()の結果はcv.resultクラスのリスト。

ベストモデルの１番目のモデルから、予測結果を取得。

10番目のモデルをcv.resultオブジェクトとして取り出し。

10番目のモデルの予測結果を取得。

技術的詳細

モデルごとの注意事項

gbm

n.treesの指定

gbm呼び出しのn.treesは10だが、

grid.predictではn.treesの候補として5、10、50、80を指定している。

そのような場合、gbmの呼び出しに使われたn.treesの値（10）以上の値を

grid.predictに指定しても、予測性能はn.trees = 10の時と同じになる。

substituteでcallを作成。

OK!

GLMも同様。

ranger::ranger

glmmML::glmmML

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

既知の問題

引用文献

更新履歴

R Package Documentation

Browse R Packages

We want your feedback!

Marchen/cv.models
Model Selection and Hyper-Parameter Tuning by Cross Validation

`ranger::ranger`

`glmmML::glmmML`