library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
yardstick::spec(),
yardstick::precision(),
yardstick::recall(),
)
library(tidyverse)
library(MASS)
library(broom) # 解析結果を tibble 形式に集約
library(tidymodels)統計データ解析
第9講 判別分析 - 分析の評価
準備
以下で利用する共通パッケージを読み込む．
判別結果の評価
問題
前回と同様に東京の気候データの線形判別を行い，以下を確認しなさい．
9月と10月の気温と湿度のデータを抽出する．
tw_data <- read_csv("data/tokyo_weather.csv") tw_subset <- tw_data |> filter(month %in% c(9,10)) |> select(temp, humid, month) |> mutate(month = as_factor(month)) # 月を因子化線形判別関数を構成する．
tw_lda <- lda(formula = month ~ temp + humid, data = tw_subset)混同行列を用いて構成した判別関数の評価を行う．1
tw_lda_fitted <- predict(tw_lda) # あてはめ値を求める tw_subset |> mutate(fitted = tw_lda_fitted[["class"]]) # あてはめ値を加えたデータフレームROC曲線を描画しAUCを求める．2
tw_subset |> bind_cols(tw_lda_fitted[["posterior"]])
1 関数 lda() の返値の class は vector 型なので，関数 mutate() を用いてデータフレームの列として加えることができる．
2 関数 lda() の返値の posterior は matrix 型なので，関数 bind_cols() を用いてデータフレームとして結合することができる．
解答例
データを整理する．
tw_data <- read_csv("data/tokyo_weather.csv")
tw_subset <- tw_data |>
filter(month %in% c(9,10)) |> # 9,10月のデータ
select(temp, humid, month) |> # 気温・湿度・月を選択
mutate(month = as_factor(month)) # 月を因子化判別関数を作成する．
tw_lda <- lda(formula = month ~ temp + humid,
data = tw_subset)混同行列を用いて判別結果を評価する．
tw_lda_fitted <- predict(tw_lda)
tw_lda_cm <- tw_subset |>
bind_cols(fitted = tw_lda_fitted[["class"]]) |> # 関数 mutate() でも良い
conf_mat(truth = month, estimate = fitted)
tw_lda_cm # 簡便な表示 Truth
Prediction 9 10
9 28 1
10 2 30summary(tw_lda_cm) # 詳細な表示# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy binary 0.951
2 kap binary 0.902
3 sens binary 0.933
4 spec binary 0.968
5 ppv binary 0.966
6 npv binary 0.938
7 mcc binary 0.902
8 j_index binary 0.901
9 bal_accuracy binary 0.951
10 detection_prevalence binary 0.475
11 precision binary 0.966
12 recall binary 0.933
13 f_meas binary 0.949関数 ggplot2::autoplot() を用いて視覚化する．
autoplot(tw_lda_cm, type = "mosaic") # モザイクプロット
autoplot(tw_lda_cm, type = "heatmap") # 行列表示
ROC曲線の描画を行う．
tw_subset |>
bind_cols(tw_lda_fitted[["posterior"]]) |>
roc_curve(truth = month, `9`) |> # 9月(`9`列)への判別を陽性とする
autoplot()
AUCを計算する．
tw_subset |>
bind_cols(tw_lda_fitted[["posterior"]]) |>
roc_auc(truth = month, `9`)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.991以下は12ヶ月分のデータを用いた例で，これを参考に説明変数は適宜選択して検討してほしい．
データを整理する．
tw_subset12 <- tw_data |>
select(temp, solar, wind, humid, month) |>
mutate(month = as_factor(month))Jan, Feb などの文字を扱いたい場合は関数 as_factor() を用いるかわりに month(month, label = TRUE) とすればよい． その場合，列名の指定なども ‘Jan:Dec’ などと変わるので注意する．
判別関数を作成する．
tw_lda12 <- lda(month ~ ., # 右辺の . は month 以外の全てを説明変数として指定
data = tw_subset12)判別結果を評価する．
tw_lda12_fitted <- predict(tw_lda12)
tw_lda12_cm <- tw_subset12 |>
bind_cols(fitted = tw_lda12_fitted[["class"]]) |>
conf_mat(truth = month, estimate = fitted)
tw_lda12_cm # 表示 Truth
Prediction 1 2 3 4 5 6 7 8 9 10 11 12
1 23 9 0 0 0 0 0 0 0 0 0 4
2 4 11 1 0 0 0 0 0 0 0 0 1
3 1 2 17 5 2 0 0 0 0 2 4 3
4 0 0 4 11 3 0 0 0 0 0 0 0
5 0 0 3 6 18 3 0 0 0 0 0 0
6 0 0 0 0 4 19 3 0 3 1 1 0
7 0 0 0 0 1 3 18 7 5 0 0 0
8 0 0 0 0 0 0 2 18 6 0 0 0
9 0 0 0 0 0 3 8 6 15 1 1 0
10 0 0 0 1 2 2 0 0 1 25 8 1
11 0 0 4 7 1 0 0 0 0 2 7 4
12 3 6 2 0 0 0 0 0 0 0 9 18summary(tw_lda12_cm) # 詳細表示# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.548
2 kap multiclass 0.507
3 sens macro 0.545
4 spec macro 0.959
5 ppv macro 0.552
6 npv macro 0.959
7 mcc multiclass 0.508
8 j_index macro 0.504
9 bal_accuracy macro 0.752
10 detection_prevalence macro 0.0833
11 precision macro 0.552
12 recall macro 0.545
13 f_meas macro 0.541 autoplot(tw_lda12_cm, type = "mosaic") # モザイクプロット
autoplot(tw_lda12_cm, type = "heatmap") # 行列表示
tw_subset12 |>
bind_cols(tw_lda12_fitted[["posterior"]]) |>
roc_curve(truth = month, `1`:`12`) |>
autoplot()
tw_subset12 |>
bind_cols(tw_lda12_fitted[["posterior"]]) |>
roc_auc(truth = month, `1`:`12`)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc hand_till 0.935訓練・試験データによる評価
問題
Wine Quality Data Set を用いて，以下を確認しなさい．
8:2の比率で訓練データと試験データに分割する ．
wq_split <- initial_split(wq_data, prop = 0.8, strata = grade) training(wq_split)) testing(wq_split)訓練データを用いて線形・2次判別関数を構成する．
訓練データを用いて評価を行う．
試験データを用いて評価を行う．
解答例
データを整理する．
wq_data <-
read_delim("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv",
delim = ";") |>
mutate(grade = factor(case_when( # quality を A,B,C,D に割り当てる
quality >= 7 ~ "A",
quality >= 6 ~ "B",
quality >= 5 ~ "C",
.default = "D"
)))訓練・試験のためにデータを分割する．
set.seed(987987) # 適宜シード値は設定する
wq_split <- initial_split(wq_data, prop = 0.8,
strata = grade)オプション strata (層別) を指定すると分割したデータの grade の比率が保たれる．
線形および2次判別関数を作成する．ただし，grade のもとになっている quality は除く．
wq_formula <- grade ~ . - quality
wq_lda <- lda(formula = wq_formula, data = training(wq_split))
wq_qda <- qda(formula = wq_formula, data = training(wq_split))訓練データを用いて判別結果を評価する．
wq_lda_train_cm <- training(wq_split) |>
bind_cols(fitted = predict(wq_lda)[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
summary(wq_lda_train_cm) # 線形判別の評価# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.615
2 kap multiclass 0.380
3 sens macro 0.484
4 spec macro 0.844
5 ppv macro 0.534
6 npv macro 0.848
7 mcc multiclass 0.381
8 j_index macro 0.328
9 bal_accuracy macro 0.664
10 detection_prevalence macro 0.25
11 precision macro 0.534
12 recall macro 0.484
13 f_meas macro 0.501autoplot(wq_lda_train_cm, type = "heatmap") +
labs(title = "Linear Discriminant (training data)")
wq_qda_train_cm <- training(wq_split) |>
bind_cols(fitted = predict(wq_qda)[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
summary(wq_qda_train_cm) # 2次判別の評価# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.626
2 kap multiclass 0.411
3 sens macro 0.549
4 spec macro 0.852
5 ppv macro 0.568
6 npv macro 0.852
7 mcc multiclass 0.411
8 j_index macro 0.401
9 bal_accuracy macro 0.700
10 detection_prevalence macro 0.25
11 precision macro 0.568
12 recall macro 0.549
13 f_meas macro 0.556autoplot(wq_qda_train_cm, type = "heatmap") +
labs(title = "Quadratic Discriminant (training data)")
試験データを用いて判別結果を評価する．
wq_lda_test_cm <- testing(wq_split) |>
bind_cols(fitted = predict(wq_lda,
newdata = testing(wq_split))[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
summary(wq_lda_test_cm) # 線形判別の評価# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.607
2 kap multiclass 0.379
3 sens macro 0.485
4 spec macro 0.843
5 ppv macro 0.569
6 npv macro 0.847
7 mcc multiclass 0.382
8 j_index macro 0.329
9 bal_accuracy macro 0.664
10 detection_prevalence macro 0.25
11 precision macro 0.569
12 recall macro 0.485
13 f_meas macro 0.508autoplot(wq_lda_test_cm, type = "heatmap") +
labs(title = "Linear Discriminant (test data)")
wq_qda_test_cm <- testing(wq_split) |>
bind_cols(fitted = predict(wq_qda,
newdata = testing(wq_split))[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
summary(wq_qda_test_cm) # 2次判別の評価# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.604
2 kap multiclass 0.383
3 sens macro 0.502
4 spec macro 0.844
5 ppv macro 0.541
6 npv macro 0.847
7 mcc multiclass 0.383
8 j_index macro 0.346
9 bal_accuracy macro 0.673
10 detection_prevalence macro 0.25
11 precision macro 0.541
12 recall macro 0.502
13 f_meas macro 0.513autoplot(wq_qda_test_cm, type = "heatmap") +
labs(title = "Quadratic Discriminant (test data)")
LOO交叉検証法による予測誤差の評価
問題
データセット MASS::biopsy を用いて2次判別の分析を行いなさい．
bio_data <- as_tibble(biopsy) |>
na.omit() |> # NA を除く
select(-ID) # IDを除く全てのデータを用いて訓練誤差を評価する．
LOO交叉検証法を用いて予測誤差を評価する．
bio_qda_loo <- qda(class ~ ., data = bio_data, CV = TRUE)
解答例
データを整理する．
bio_data <- as_tibble(biopsy) |>
na.omit() |> # NA を除く
select(-ID) # IDを除く散布図で項目間の関係を確認する．
bio_data |>
GGally::ggpairs(diag = list(mapping = aes(colour = class),
lower = list(mapping = aes(colour = class))))
以下は2次判別の LOO CV による評価の例である． 線形判別も同様に行うことができる．
bio_formula <- class ~ .
bio_qda <- qda(formula = bio_formula, data = bio_data)
bio_qda_loo <- qda(formula = bio_formula, data = bio_data, CV = TRUE)訓練誤差での評価を確認し，可視化する．
bio_qda_cm <- bio_data |>
bind_cols(fitted = predict(bio_qda)[["class"]]) |>
conf_mat(truth = class, estimate = fitted)
summary(bio_qda_cm) # 2次判別によるあてはめ値の評価(訓練誤差)# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy binary 0.959
2 kap binary 0.911
3 sens binary 0.950
4 spec binary 0.975
5 ppv binary 0.986
6 npv binary 0.914
7 mcc binary 0.912
8 j_index binary 0.925
9 bal_accuracy binary 0.963
10 detection_prevalence binary 0.627
11 precision binary 0.986
12 recall binary 0.950
13 f_meas binary 0.968autoplot(bio_qda_cm, type = "heatmap") +
labs(title = "Training Error")
LOO CV での評価を確認し，可視化する．
bio_qda_loo_cm <- bio_data |>
bind_cols(fitted = bio_qda_loo[["class"]]) |>
conf_mat(truth = class, estimate = fitted)
summary(bio_qda_loo_cm) # LOO CVによる予測の評価(予測誤差)# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy binary 0.950
2 kap binary 0.893
3 sens binary 0.937
4 spec binary 0.975
5 ppv binary 0.986
6 npv binary 0.893
7 mcc binary 0.895
8 j_index binary 0.912
9 bal_accuracy binary 0.956
10 detection_prevalence binary 0.618
11 precision binary 0.986
12 recall binary 0.937
13 f_meas binary 0.961autoplot(bio_qda_loo_cm, type = "heatmap") +
labs(title = "Test Error (LOO CV)")
あてはめ値による評価は LOO CV より若干良くなっており，あてはめ値では精度を過剰に評価する可能性があることが示唆される．
k-重交叉検証法による予測誤差の評価
問題
Wine Quality Data Set を用いて 線形判別と2次判別の分析を行いなさい．
LOO交叉検証法を用いて予測誤差を評価する．
wq_lda_loo <- lda(formula = wq_formula, data = training(wq_split), CV = TRUE) wq_qda_loo <- qda(formula = wq_formula, data = training(wq_split), CV = TRUE)k-重交叉検証法を用いて予測誤差を評価する．
解答例
前の問題で既に整理してある wq_data/wq_split を用いる． wq_lda と比較する．
wq_lda_loo <- lda(formula = wq_formula, data = training(wq_split),
CV = TRUE)
training(wq_split) |> # 訓練誤差による評価
bind_cols(fitted = predict(wq_lda)[["class"]]) |>
conf_mat(truth = grade, estimate = fitted) |>
autoplot(type = "heatmap") +
labs(title = "Training Error (LDA)")
wq_lda_loo_cm <- training(wq_split) |> # LOO CV 予測誤差による評価
bind_cols(fitted = wq_lda_loo[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
wq_lda_loo_cm |> autoplot(type = "heatmap") +
labs(title = "LOO CV (LDA)")
線形判別の過学習は微小であることがわかる．
2次判別の LOO CV を wq_qda と比較する．
wq_qda_loo <- qda(formula = wq_formula, data = training(wq_split),
CV = TRUE)
training(wq_split) |> # 訓練誤差による評価
bind_cols(fitted = predict(wq_qda)[["class"]]) |>
conf_mat(truth = grade, estimate = fitted) |>
autoplot(type = "heatmap") +
labs(title = "Training Error (QDA)")
wq_qda_loo_cm <- training(wq_split) |> # LOO CV 予測誤差による評価
bind_cols(fitted = wq_qda_loo[["class"]]) |>
conf_mat(truth = grade, estimate = fitted)
wq_qda_loo_cm |> autoplot(type = "heatmap") +
labs(title = "LOO CV (QDA)")
2次判別は若干過学習していることが示唆される．
LOO CV により線形・2次判別の予測誤差を比較する．
summary(wq_lda_loo_cm) # 線形# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.600
2 kap multiclass 0.357
3 sens macro 0.465
4 spec macro 0.839
5 ppv macro 0.507
6 npv macro 0.842
7 mcc multiclass 0.357
8 j_index macro 0.304
9 bal_accuracy macro 0.652
10 detection_prevalence macro 0.25
11 precision macro 0.507
12 recall macro 0.465
13 f_meas macro 0.479summary(wq_qda_loo_cm) # 2次# A tibble: 13 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.576
2 kap multiclass 0.330
3 sens macro 0.467
4 spec macro 0.832
5 ppv macro 0.481
6 npv macro 0.833
7 mcc multiclass 0.330
8 j_index macro 0.300
9 bal_accuracy macro 0.650
10 detection_prevalence macro 0.25
11 precision macro 0.481
12 recall macro 0.467
13 f_meas macro 0.473予測誤差の観点からは線形判別の方が良さそうであることが伺える．
package::tidymodels を用いた k-重交叉検証は以下のようになる．
まず lda/qda を tidymodels 用に宣言する．
library(discrim) # 以下の判別モデルを設定するために必要
tidy_qda <- discrim_quad() |> set_engine("MASS") |> set_mode("classification")
tidy_lda <- discrim_linear() |> set_engine("MASS") |> set_mode("classification")交叉検証用にデータ分割を行う．
wq_folds <- vfold_cv(training(wq_split),
v = 10) # 10-fold を指定 (既定値)評価指標を設定する．
wq_metrics <- metric_set(accuracy, # 精度
sens, # 感度 (真陽性率)
spec, # 特異度 (真陰性率)
precision, # 適合率
recall, # 再現率(sensと同じ)
roc_auc, # AUC
kap, # kappa
f_meas) # f値共通の処理を定義する．
wq_workflow <- workflow() |>
add_formula(wq_formula) 線形判別を行い，評価する．
wq_lda_cv <- wq_workflow |>
add_model(tidy_lda) |>
fit_resamples(resamples = wq_folds,
metrics = wq_metrics)
wq_lda_cv |> collect_metrics()# A tibble: 8 × 6
.metric .estimator mean n std_err .config
<chr> <chr> <dbl> <int> <dbl> <chr>
1 accuracy multiclass 0.597 10 0.00753 Preprocessor1_Model1
2 f_meas macro 0.467 10 0.0221 Preprocessor1_Model1
3 kap multiclass 0.349 10 0.0130 Preprocessor1_Model1
4 precision macro 0.511 10 0.0339 Preprocessor1_Model1
5 recall macro 0.456 10 0.0187 Preprocessor1_Model1
6 roc_auc hand_till 0.784 10 0.0109 Preprocessor1_Model1
7 sens macro 0.456 10 0.0187 Preprocessor1_Model1
8 spec macro 0.837 10 0.00291 Preprocessor1_Model12次判別を行い，評価する．
wq_qda_cv <- wq_workflow |>
add_model(tidy_qda) |>
fit_resamples(resamples = wq_folds,
metrics = wq_metrics)
wq_qda_cv |> collect_metrics()# A tibble: 8 × 6
.metric .estimator mean n std_err .config
<chr> <chr> <dbl> <int> <dbl> <chr>
1 accuracy multiclass 0.568 10 0.00958 Preprocessor1_Model1
2 f_meas macro 0.445 10 0.0217 Preprocessor1_Model1
3 kap multiclass 0.317 10 0.0171 Preprocessor1_Model1
4 precision macro 0.460 10 0.0256 Preprocessor1_Model1
5 recall macro 0.447 10 0.0201 Preprocessor1_Model1
6 roc_auc hand_till 0.732 10 0.0192 Preprocessor1_Model1
7 sens macro 0.447 10 0.0201 Preprocessor1_Model1
8 spec macro 0.830 10 0.00401 Preprocessor1_Model1