16章

評価にどのデータを用いるかを決める際に大きく

の2通りが存在する。

16.1

16.1.1

再代入法：評価指標をモデル構築の際に用いたデータに適用すること。 モデル適合度の推定値を得ることができるが、再代入法はモデルの予測性能を評価することはできない。

16.1.2

ランダム化によるモデル当てはめの評価には様々な手法がある。 原理はRandomization testと全く同じ。 また、ランダム化するのは応答変数のみでなく説明変数でも可能。

16.2

再サンプリングには様々な手法がある。 サンプル分割やk重交差検証法、1個抜き検証法、ブートストラップ、0.632+ブートストラップ、ジャックナイフ…など。

16.2.1

k重交差検証法（k fold cross-validation）：データをk個に分割し、k-1個のデータで訓練したモデルを残りの1個のデータで当てはめた時の指標を産出する。この操作をk回繰り返して指標の平均を産出する。

s_mammals_data <- read.csv("hsdm/data/tabular/species/summary_mammals_and_bioclim.csv", row.names=1)
library(boot)

set.seed(555)
cv.error.10=rep(0,10)
for (i in 1:10){
    glm.fit=glm(VulpesVulpes~poly(bio3+bio4+bio7+bio11+bio12,i),data=s_mammals_data)
    cv.error.10[i]=cv.glm(s_mammals_data,glm.fit,K=10)$delta[1]
}
cv.error.10
##  [1] 0.1508741 0.1399128 0.1386917 0.1374293 0.1327043 0.1327239 0.1316569
##  [8] 0.1315851 0.1318129 0.1315490

cv.glm()でk重交差検証法を実行できる。

delta: The first component is the raw cross-validation estimate of prediction error. The second component is the adjusted cross-validation estimate.

delta[1]で指定しているのは生のcross-validation推定値らしい。違いがそこまでわかっていない。

library(tidyverse)
library(tidymodels)
#   ref: https://stateofther.github.io/finistR2019/s-tidymodels.html


smd_split <- s_mammals_data %>% 
    mutate(VulpesVulpes=if_else(VulpesVulpes==0, 1, 0), 
           VulpesVulpes=as.factor(VulpesVulpes)) %>% initial_split(0.9)
smd_train <- training(smd_split)
smd_test <- testing(smd_split)

#   logit_fit <- logistic_reg(mode = "classification") %>%
#       set_engine(engine = "glm") %>% 
#       fit(VulpesVulpes ~ bio3+bio7+bio11+bio12, data = smd_train) %>%
#       step_poly(bio3, bio7, bio11, bio12)
#   
#   pred_logit_train <- predict(logit_fit, new_data = smd_train)
#   pred_logit_test <- predict(logit_fit, new_data = smd_test)
#   prob_logit_test <- predict(logit_fit, new_data = smd_test, type="prob")
#   
#   metrics(bind_cols(smd_test, pred_logit_test), truth = VulpesVulpes, estimate = .pred_class)
#   

multimetric <- metric_set(accuracy, bal_accuracy, sens, yardstick::spec, precision, recall, ppv, npv)
#   multimetric(bind_cols(forest_test, pred_logit_test), truth = VulpesVulpes, estimate = .pred_class)

cv_train <- vfold_cv(smd_train, v = 10, repeats = 9, strata = "VulpesVulpes")

logit_mod <- 
    logistic_reg(mode = "classification") %>%
    set_engine(engine = "glm")

## compute mod on kept part
cv_fit <- function(splits, mod, ...) {
    res_mod <-
        fit(mod, VulpesVulpes ~ bio3+bio7+bio11+bio12, data = analysis(splits), family = binomial)
    return(res_mod)
}

## get predictions on holdout sets
cv_pred <- function(splits, mod){
    # Save the 10%
    holdout <- assessment(splits)
    pred_assess <- bind_cols(truth = holdout$VulpesVulpes, predict(mod, new_data = holdout))
    return(pred_assess)
}

## get probs on holdout sets
cv_prob <- function(splits, mod){
    holdout <- assessment(splits)
    prob_assess <- bind_cols(truth = as.factor(holdout$VulpesVulpes), 
                             predict(mod, new_data = holdout, type = "prob"))
    return(prob_assess)
}

res_cv_train <- cv_train %>% 
    mutate(res_mod = map(splits, .f = cv_fit, logit_mod), ## fit model
           res_pred = map2(splits, res_mod, .f = cv_pred), ## predictions
           res_prob = map2(splits, res_mod, .f = cv_prob)  ## probabilities
    )

#   res_cv_train %>% 
#       mutate(metrics = map(res_pred, multimetric, truth = truth, estimate = .pred_class)) %>% 
#       unnest(metrics) %>% 
#       ggplot() + 
#       aes(x = id, y = .estimate) +
#       geom_point() + 
#       facet_wrap(~ .metric, scales = "free_y") 

res_cv_train %>% 
    mutate(roc = map(res_prob, roc_curve, truth = truth, .pred_1)) %>% 
    unnest(roc) %>% 
    ggplot() +
    aes(x = 1 - specificity, y = sensitivity, color = id2) +
    geom_path() +
    geom_abline(lty = 3) + facet_wrap(~id)

AUCがすごいことになっている。きちんと確認していないが、factorの01の順番に依存している可能性がある（岡本くんだん）。

vulpes_data<-s_mammals_data %>% select(bio3, bio4, bio7, bio11, bio12, VulpesVulpes) %>%
    mutate(VulpesVulpes=as.factor(VulpesVulpes))

myRF <- function(formula, train, test){
    model <- randomForest(formula, train)
    predict(model,test,type="prob")[,"pos"]
}

set.seed(555)

vulpes_RF_cv <- Daim(formula=VulpesVulpes~., model=myRF, data=vulpes_data, labpos="1", control=Daim.control(method="cv", k=10, k.runs=10), cutoff="cv")

vulpes_RF_cv
summary(vulpes_RF_cv)
auc(vulpes_RF_cv)$auc.loob
auc(vulpes_RF_cv)$auc.samples
vulpes_data<-s_mammals_data %>% select(bio3, bio4, bio7, bio11, bio12, VulpesVulpes) %>%
    mutate(VulpesVulpes=as.factor(VulpesVulpes))

smd_split <- vulpes_data %>% 
    as_tibble() %>%
    initial_split(0.9)
smd_train <- training(smd_split)
smd_test <- testing(smd_split)

16.2.2

1個抜き交差検証法（ジャックナイフ）：k-fold cross-validationのkをn（全データ数）にしたもの（と理解している）。 当然、k-foldよりも計算時間がかかる。

実行結果は省略。

16.2.3

繰り返し分割交差検証法：k-fold validationをR回(R>>k)繰り返す手法。リサンプリングありのBootstrapに近いイメージ。

生息適地モデリング研究で一般的によく使われている手法。また、前述の2つの検証法に比べて不確実性や安定性も評価できる。

set.seed(555)
train=sample(2488,1244)

attach(s_mammals_data)
## The following objects are masked from s_mammals_data (pos = 3):
## 
##     bio11, bio12, bio3, bio4, bio7, ConnochaetesGnou, GuloGulo,
##     PantheraOnca, PteropusGiganteus, TenrecEcaudatus, VulpesVulpes,
##     X_WGS84, Y_WGS84
## The following objects are masked from s_mammals_data (pos = 4):
## 
##     bio11, bio12, bio3, bio4, bio7, ConnochaetesGnou, GuloGulo,
##     PantheraOnca, PteropusGiganteus, TenrecEcaudatus, VulpesVulpes,
##     X_WGS84, Y_WGS84
## The following objects are masked from s_mammals_data (pos = 10):
## 
##     bio11, bio12, bio3, bio4, bio7, ConnochaetesGnou, GuloGulo,
##     PantheraOnca, PteropusGiganteus, TenrecEcaudatus, VulpesVulpes,
##     X_WGS84, Y_WGS84
glm.fit=glm(VulpesVulpes~bio3+bio7+bio11+bio12,family="binomial",data=s_mammals_data,subset=train)

mean((VulpesVulpes-predict(glm.fit,s_mammals_data))[-train]^2)
## [1] 10.62033
glm.fit2=glm(VulpesVulpes~poly(bio3+bio7+bio11+bio12,2),family="binomial",data=s_mammals_data,subset=train)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
mean((VulpesVulpes-predict(glm.fit2,s_mammals_data))[-train]^2)
## [1] 9.820958
glm.fit3=glm(VulpesVulpes~poly(bio3+bio7+bio11+bio12,3),family="binomial",data=s_mammals_data,subset=train)
mean((VulpesVulpes-predict(glm.fit3,s_mammals_data))[-train]^2)
## [1] 4.454756

なぜか教科書と結果が異なる…。わからん。 3次の多項式GLMが一番当てはまりが良いことは一緒。

16.3

完全に独立なデータを用いてテストを行う場合、「完全に独立なデータ」とは何かを気にする必要がある。

完全に独立なデータには4つの状況に分けることができる。

  1. 同じ時期で同じ地域
  2. 同じ時期だが異なる地域
  3. 異なる時期だが同じ時期
  4. 時期も地域も異なる場合

4.は2.3.の合わせ技でなんとかなるので割愛。 1.の問題は学習データセットと空間的に独立でないこと。 2.と3.については生物学的な変化に留意する必要がある。 具体的には、2.では遺伝的な分化とそれによるニッチの多様性が生じている可能性。 3.では未知の要因が個体群動態に影響を与えている可能性を考慮する必要がある。