情報科学3
Rによる分類の例
矢内 勇生
2021-09-30
Rによる機械学習の例として、分類を行う方法を紹介する。各手法の説明は、今後の授業で行う。今回は、（どれくらいうまくいくかは別として）Rで分類が行えるということを確認する。
準備
必要なパッケージを読み込む。パッケージの読み込みにはpacman::p_load() を使う。pacmanパッケージがインストール済みでない場合は、まずインストールする。
if (!require(pacman)) install.packages("pacman")
pacman::p_load(tidyverse,
palmerpenguins, # データソース
caret, # 混同行列, kNN
MASS, # LDA
nnet, # 多項ロジスティック回帰
e1071, # SVM
rpart, # 決定木
randomForest, # ランダムフォレスト
scales) # expand_range()次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。LinuxにはIPAexフォント がインストールされていることを想定している（IPAex はインストールすれば maxOSやWindows でも使える）。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
my_font <- "Yu Gothic"
} else {
my_font <- "Japan1"
}
} else if (capabilities("aqua")) {
my_font <- "HiraginoSans-W3"
} else {
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))palmerpenguins データ
今回は、palmerpenguinsパッケージのpenguinsデータを使う。
data(penguins)データセットの詳細は?penguins で確認できる。
いくつかの方法で、データの中身を確認してみよう。 まず、データセットのサイズ（次元, 行数と列数）を確認する。
dim(penguins)## [1] 344 8行数（観測数, サンプルサイズ）。
nrow(penguins)## [1] 344列数（変数の数）。
ncol(penguins)## [1] 8変数名（列名）。
names(penguins)## [1] "species" "island" "bill_length_mm"
## [4] "bill_depth_mm" "flipper_length_mm" "body_mass_g"
## [7] "sex" "year"データセットの基本情報を確認する。データセットを読み込んだ直後に、毎回これを実行することを推奨。ただし、R Markdown ファイルには書かない（あるいはチャンクオプションで eval = FALSE を指定する）ほうがいいかもしれない。
glimpse(penguins)## Rows: 344
## Columns: 8
## $ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex <fct> male, female, female, NA, female, male, female, male…
## $ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…同上（tidyverse を読み込まなくても使える）。
str(penguins)## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)
## $ species : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ island : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bill_length_mm : num [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
## $ bill_depth_mm : num [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
## $ flipper_length_mm: int [1:344] 181 186 195 NA 193 190 181 195 193 190 ...
## $ body_mass_g : int [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
## $ sex : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
## $ year : int [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...データセットの先頭部（最初の\(n\) 行。デフォルトはn=6）を表示。
head(penguins)## # A tibble: 6 × 8
## species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex
## <fct> <fct> <dbl> <dbl> <int> <int> <fct>
## 1 Adelie Torge… 39.1 18.7 181 3750 male
## 2 Adelie Torge… 39.5 17.4 186 3800 fema…
## 3 Adelie Torge… 40.3 18 195 3250 fema…
## 4 Adelie Torge… NA NA NA NA <NA>
## 5 Adelie Torge… 36.7 19.3 193 3450 fema…
## 6 Adelie Torge… 39.3 20.6 190 3650 male
## # … with 1 more variable: year <int>データセットの最後部（最後の\(n\) 行。デフォルトはn=6）を表示。
tail(penguins, n = 3)## # A tibble: 3 × 8
## species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex
## <fct> <fct> <dbl> <dbl> <int> <int> <fct>
## 1 Chinst… Dream 49.6 18.2 193 3775 male
## 2 Chinst… Dream 50.8 19 210 4100 male
## 3 Chinst… Dream 50.2 18.7 198 3775 fema…
## # … with 1 more variable: year <int>基本的な記述統計を表示。
summary(penguins)## species island bill_length_mm bill_depth_mm
## Adelie :152 Biscoe :168 Min. :32.10 Min. :13.10
## Chinstrap: 68 Dream :124 1st Qu.:39.23 1st Qu.:15.60
## Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30
## Mean :43.92 Mean :17.15
## 3rd Qu.:48.50 3rd Qu.:18.70
## Max. :59.60 Max. :21.50
## NA's :2 NA's :2
## flipper_length_mm body_mass_g sex year
## Min. :172.0 Min. :2700 female:165 Min. :2007
## 1st Qu.:190.0 1st Qu.:3550 male :168 1st Qu.:2007
## Median :197.0 Median :4050 NA's : 11 Median :2008
## Mean :200.9 Mean :4202 Mean :2008
## 3rd Qu.:213.0 3rd Qu.:4750 3rd Qu.:2009
## Max. :231.0 Max. :6300 Max. :2009
## NA's :2 NA's :2問題の設定
今回は、くちばしの長さ (bill_length_mm) と くちばしの厚み (bill_depth_mm) からペンギンの種類 (species) を予測（分類）する。
with(penguins, table(species))## species
## Adelie Chinstrap Gentoo
## 152 68 124種類が英語で記録されているので日本語に変え、利用する3つの変数のなかに欠測値があるものを（とりあえず）取り除く。
penguins_s1 <- penguins %>%
mutate(species = case_when(
species == "Adelie" ~ 1L,
species == "Chinstrap" ~ 2L,
TRUE ~ 3L,
)) %>%
mutate(species = factor(species,
levels = 1:3,
labels = c("アデリー", "ヒゲ", "ジェンツー"))) %>%
dplyr::select(species, bill_length_mm, bill_depth_mm) %>%
na.omit()予測に使う2つの変数（機械学習ではこれを特徴量 (feature)と呼ぶ）と、実際の種類（つまり、分類問題の正解。機械学習における分類問題ではこれをラベル (label)と呼ぶ）の関係を図示する。
p1 <- ggplot(penguins_s1,
aes(x = bill_length_mm,
y = bill_depth_mm)) +
geom_point(aes(color = species,
shape = species)) +
scale_color_brewer(palette = "Accent") +
labs(x = "くちばしの長さ (mm)",
y = "くちばしの厚み (mm)",
color = "ペンギンの種類",
shape = "ペンギンの種類")
plot(p1)
予測に使うRのformulaを用意しておこう。 species を bill_length_mm と bill_depth_mm から予測する（species を bill_length_mm と bill_depth_mm に回帰する）。
pgnf <- formula(species ~ bill_length_mm + bill_depth_mm)分類の境界線 (decision boundaries; discrimination lines) を求める範囲のグリッドデータを用意しておく。
bill_len_lim <- range(penguins_s1$bill_length_mm) %>%
expand_range(mul = 0.05)
bill_dep_lim <- range(penguins_s1$bill_depth_mm) %>%
expand_range(mul = 0.05)
newdf <- expand_grid(
bill_length_mm = seq(bill_len_lim[1], bill_len_lim[2], length.out = 500),
bill_depth_mm = seq(bill_dep_lim[1], bill_dep_lim[2], length.out = 500)
)
glimpse(newdf, n = 5)## Rows: 250,000
## Columns: 2
## $ bill_length_mm <dbl> 30.725, 30.725, 30.725, 30.725, 30.725, 30.725, 30.725,…
## $ bill_depth_mm <dbl> 12.68000, 12.69852, 12.71703, 12.73555, 12.75407, 12.77…これからいくつかの方法で分類を行い、そのたびに上で作った散布図に、分類の境界線を上書きした図を作る。繰り返し行う作業は関数化してしまったほうが効率が良い。 そこで、分類結果を図示するための関数decision_region()を自作する。
decision_region <- function(pred_class, title = NULL) {
# pred_class には予測（分類）結果を数値で表したベクトルを渡す
df <- newdf %>%
mutate(predicted = pred_class)
p <- p1 +
geom_raster(data = df,
aes(fill = as.factor(predicted)),
alpha = 0.2,
show.legend = FALSE) +
geom_contour(data = df,
aes(z = predicted),
color = "gray") +
scale_fill_brewer(palette = "Accent") +
scale_x_continuous(limits = bill_len_lim, expand = c(0, 0)) +
scale_y_continuous(limits = bill_dep_lim, expand = c(0, 0)) +
ggtitle(title)
return(p)
}LDA（線形判別分析）
LDA (Linear Discriminant Analysis) で分類してみよう。MASS::lda()を使う。
classifier_lda <- lda(pgnf, data = penguins_s1)予測（分類）がどれだけうまくできたか確認するために、混同行列 (confusion matrix) を表示する。 predict() で予測を行い、caret::confusionMatrix() で混同行列を作る。
pred_lda <- predict(classifier_lda)
confusionMatrix(pred_lda$class,
reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 5 0
## ヒゲ 2 59 2
## ジェンツー 0 4 121
##
## Overall Statistics
##
## Accuracy : 0.962
## 95% CI : (0.9359, 0.9796)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.94
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.8676 0.9837
## Specificity 0.9738 0.9854 0.9817
## Pos Pred Value 0.9675 0.9365 0.9680
## Neg Pred Value 0.9894 0.9677 0.9908
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1725 0.3538
## Detection Prevalence 0.4503 0.1842 0.3655
## Balanced Accuracy 0.9803 0.9265 0.9827Prediction が予測結果、Reference が正解である。Accuracy に予測の正解率が示されている。他の部分の読み方については、今後の授業で説明する。
分類の境界線 (decision boundaries; discrimination lines) を図示する。
p_lda <- predict(classifier_lda, newdata = newdf)$class %>%
as.numeric() %>%
decision_region(title = "LDAによる分類")
plot(p_lda)
ロジスティック回帰
ロジスティック回帰による分類を行う。分類するカテゴリが3つあるので、多項ロジスティック回帰 (multinomial logistic regression) を利用する。nnet::multinom()を使う。
classifier_logit <- multinom(pgnf, data = penguins_s1)## # weights: 12 (6 variable)
## initial value 375.725403
## iter 10 value 26.783497
## iter 20 value 24.005561
## iter 30 value 23.973372
## iter 40 value 23.956920
## iter 50 value 23.951745
## iter 60 value 23.950452
## final value 23.948361
## converged混同行列。
predict(classifier_logit) %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 148 3 0
## ヒゲ 3 61 2
## ジェンツー 0 4 121
##
## Overall Statistics
##
## Accuracy : 0.9649
## 95% CI : (0.9395, 0.9817)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9448
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9801 0.8971 0.9837
## Specificity 0.9843 0.9818 0.9817
## Pos Pred Value 0.9801 0.9242 0.9680
## Neg Pred Value 0.9843 0.9746 0.9908
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4327 0.1784 0.3538
## Detection Prevalence 0.4415 0.1930 0.3655
## Balanced Accuracy 0.9822 0.9394 0.9827分類の境界線を図示する。
p_logit <- predict(classifier_logit, newdata = newdf) %>%
as.numeric() %>%
decision_region(title = "多項ロジスティック回帰による分類")
plot(p_logit)
SVM（サポートベクタマシン）
SVM (support vector machine) による分類を行う。 e1071::svm() を使う。kernel引数で利用するカーネルを指定する。
線形カーネル
線形カーネル (linear kernel) を使ってSVMで分類を行う。
svm_linear <- svm(pgnf,
data = penguins_s1,
kernel = "linear")混同行列。
predict(svm_linear) %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 5 0
## ヒゲ 2 60 3
## ジェンツー 0 3 120
##
## Overall Statistics
##
## Accuracy : 0.962
## 95% CI : (0.9359, 0.9796)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.94
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.8824 0.9756
## Specificity 0.9738 0.9818 0.9863
## Pos Pred Value 0.9675 0.9231 0.9756
## Neg Pred Value 0.9894 0.9711 0.9863
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1754 0.3509
## Detection Prevalence 0.4503 0.1901 0.3596
## Balanced Accuracy 0.9803 0.9321 0.9810e1071パッケージには、分類境界を可視化する関数が用意されている。
par(family = my_font)
plot(svm_linear, data = penguins_s1)
横軸と縦軸が入れ替わっている。
ggplot2 で描き直す。
p_svm_linear <- predict(svm_linear, newdata = newdf) %>%
as.numeric() %>%
decision_region(title = "SVM（線形カーネル）による分類")
plot(p_svm_linear)
ガウスカーネル
ガウスカーネル (Gaussian kernel or radial basis function; RBF) を使う。 gamma引数で、関数の複雑さ（円の曲率）を調整する。大きな値にするほど、複雑な関数を使った分類になる。 まず、gamma = 0.5でやってみる。
svm_rbf1 <- svm(pgnf,
data = penguins_s1,
kernel = "radial",
gamma = 0.5) # 関数の複雑さを決めるパラメタ混同行列。
predict(svm_rbf1) %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 5 0
## ヒゲ 2 60 3
## ジェンツー 0 3 120
##
## Overall Statistics
##
## Accuracy : 0.962
## 95% CI : (0.9359, 0.9796)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.94
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.8824 0.9756
## Specificity 0.9738 0.9818 0.9863
## Pos Pred Value 0.9675 0.9231 0.9756
## Neg Pred Value 0.9894 0.9711 0.9863
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1754 0.3509
## Detection Prevalence 0.4503 0.1901 0.3596
## Balanced Accuracy 0.9803 0.9321 0.9810分類境界を図示する。
p_svm_rbf1 <- predict(svm_rbf1, newdata = newdf) %>%
as.numeric() %>%
decision_region(title = "SVM（ガウシアンカーネル）による分類: gamma = 0.5")
plot(p_svm_rbf1)
境界線が直線ではなくなったことがわかる。
次に、gamma = 10 にしてみる。
svm_rbf2 <- svm(pgnf,
data = penguins_s1,
kernel = "radial",
gamma = 10) # 関数の複雑さを決めるパラメタ混同行列。
predict(svm_rbf2) %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 2 0
## ヒゲ 2 65 2
## ジェンツー 0 1 121
##
## Overall Statistics
##
## Accuracy : 0.9795
## 95% CI : (0.9583, 0.9917)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9679
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.9559 0.9837
## Specificity 0.9895 0.9854 0.9954
## Pos Pred Value 0.9868 0.9420 0.9918
## Neg Pred Value 0.9895 0.9890 0.9909
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1901 0.3538
## Detection Prevalence 0.4415 0.2018 0.3567
## Balanced Accuracy 0.9881 0.9706 0.9896分類境界を図示する。
p_svm_rbf2 <- predict(svm_rbf2, newdata = newdf) %>%
as.numeric() %>%
decision_region(title = "SVM（ガウシアンカーネル）による分類: gamma = 10")
plot(p_svm_rbf2)
先程よりも複雑な境界線が引かれている。その結果、観測データの分類がうまくできている一方で、グラフの右下や右上が緑色（アデリーペンギン）と予測（分類）されており、関数を複雑にし過ぎると弊害がありそうだということがわかる。
決定木
決定木 (decision tree) による分類を行う。 rpart::rpart() を使う。
classifier_tree <- rpart(pgnf, data = penguins_s1)混同行列。
predict(classifier_tree, type = "class") %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 148 6 4
## ヒゲ 3 62 8
## ジェンツー 0 0 111
##
## Overall Statistics
##
## Accuracy : 0.9386
## 95% CI : (0.9077, 0.9616)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9036
##
## Mcnemar's Test P-Value : 0.004637
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9801 0.9118 0.9024
## Specificity 0.9476 0.9599 1.0000
## Pos Pred Value 0.9367 0.8493 1.0000
## Neg Pred Value 0.9837 0.9777 0.9481
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4327 0.1813 0.3246
## Detection Prevalence 0.4620 0.2135 0.3246
## Balanced Accuracy 0.9639 0.9358 0.9512分類境界を図示する。
p_tree <- predict(classifier_tree,
newdata = newdf,
type = "class") %>%
as.numeric() %>%
decision_region(title = "決定木による分類")
plot(p_tree)
境界線がそれぞれの特徴量（説明変数）の軸と垂直（平行）に引かれていることがわかる。
ランダムフォレスト
ランダムフォレスト (random forest) による分類を行う。randomForest::randomForest() を使う。ntree 引数で、いくつの決定木を使うかを決める。
classifier_rf <- randomForest(pgnf,
data = penguins_s1,
ntree = 50)混同行列。
predict(classifier_rf) %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 3 0
## ヒゲ 2 63 4
## ジェンツー 0 2 119
##
## Overall Statistics
##
## Accuracy : 0.9678
## 95% CI : (0.9432, 0.9838)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9495
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.9265 0.9675
## Specificity 0.9843 0.9781 0.9909
## Pos Pred Value 0.9803 0.9130 0.9835
## Neg Pred Value 0.9895 0.9817 0.9819
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1842 0.3480
## Detection Prevalence 0.4444 0.2018 0.3538
## Balanced Accuracy 0.9855 0.9523 0.9792分類境界を図示する。
p_rf <- predict(classifier_rf,
newdata = newdf,
type = "class") %>%
as.numeric() %>%
decision_region(title = "ランダムフォレストによる分類")
plot(p_rf)
決定木と比較すると、複雑な境界線になることがわかる。
\(k\)近傍法による分類
\(k\)近傍 (\(k\)-nearest neighbor; kNN) 法による分類を行う。 caret::knn3() を使う。
classifier_kNN <- knn3(pgnf,
data = penguins_s1,
k = 5)混同行列。
predict(classifier_kNN, newdata = penguins_s1, type = "class") %>%
confusionMatrix(reference = penguins_s1$species)## Confusion Matrix and Statistics
##
## Reference
## Prediction アデリー ヒゲ ジェンツー
## アデリー 149 2 0
## ヒゲ 2 64 3
## ジェンツー 0 2 120
##
## Overall Statistics
##
## Accuracy : 0.9737
## 95% CI : (0.9506, 0.9879)
## No Information Rate : 0.4415
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9587
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: アデリー Class: ヒゲ Class: ジェンツー
## Sensitivity 0.9868 0.9412 0.9756
## Specificity 0.9895 0.9818 0.9909
## Pos Pred Value 0.9868 0.9275 0.9836
## Neg Pred Value 0.9895 0.9853 0.9864
## Prevalence 0.4415 0.1988 0.3596
## Detection Rate 0.4357 0.1871 0.3509
## Detection Prevalence 0.4415 0.2018 0.3567
## Balanced Accuracy 0.9881 0.9615 0.9832分類境界を図示する。
p_kNN <- predict(classifier_kNN,
newdata = newdf,
type = "class") %>%
as.numeric() %>%
decision_region(title = "k近傍法による分類 (k=5)")
plot(p_kNN)
このように、Rを使うとさまざまな方法で分類を行うことができる（もちろんPythonでもできるが、機械学習の基本的な問題については、Rだけでも対処できる。RもPythonも使えるのが理想）。方法ごとに異なるパッケージを利用することになるが、それぞれの方法に固有のパラメタ（引数）の指定が必要だという点を除けば、基本的な使い方はほぼ同じである。
ここで実行した分類を、引数（たとえば、svm() のgammaや knn3() のk）の値を変えて実行してみよう。どのような結果が得られるだろうか。