必要なパッケージ
require(here)
require(plotly)
require(MASS)
require(kernlab)
require(ranger)背景
科学研究などのさまざまな場面において、あるサンプル（標本）について測定された複数の特徴量から、そのサンプルが属するクラスを特定するモデルが必要とされています。このようなモデルは、クラスと測定された特徴との関係を記述するために利用できるだけでなく、新しいサンプルのクラスを予測するためにも役立ちます。
例えば、異なる品種について測定された複数の形質がある場合、その品種をクラスとみなし、形質の観測値に基づいてクラスを特定するモデルを構築することで、品種間の形質変動を記述することができます。また、リモートセンシングのデータを用いて、各地点で栽培されている作物種を特定するようなケースもあります。本講義では、このような場合に有効なクラス同定モデルの構築方法を紹介します。
分類モデルを構築するためには、さまざまな手法が存在します。本講義では、線形判別分析、サポートベクターマシン、ランダムフォレストについて取り上げます。分類モデルでは、モデルに関わる特徴量がサンプル数よりはるかに多い場合、「large p small n」問題が発生することがありますが、サポートベクターマシンやランダムフォレストは、このような問題がある場合でも有効に利用できる手法です。
線形判別分析(LDA)
この講義では、Hao et al. (2015, PLoS ONE 10: e0137748) のリモートセンシングデータを基に、主成分分析を行います。この論文では、Landsatなどの人工衛星から取得したリモートセンシング画像を用いて、対象地点で栽培されている作物を予測する機械学習モデルを構築することが研究の目的とされています。リモートセンシングは1年間に11回行われ、そこからNDVI（正規化差分植生指数）が算出されました。NDVIは以下のように算出されます。 \[ NDVI=\frac{\rho(NIR)-\rho(R)}{\rho(NIR)+\rho(R)} \] 植生は近赤外線（NIR）を反射し、赤色（レッド：R）を吸収するため、植生があるところではNDVIが高くなります。ここでは、Haoら（2015）のデータの一部を用いて、11個のNDVIデータセットからコムギを除く4種類の作物（ワタ、ブドウ、トウモロコシ、スイカ）を識別するモデルを構築してみましょう。モデル構築には1000サンプルのデータ（訓練データ）を、モデル評価には500サンプルのデータ（検証データ）を使用します。
まず、モデル構築のためのデータをロードします。
ndvi <- read.csv(here("data", "ndvi_train.csv"),
row.names = 1, stringsAsFactors = T) # Read data
head(ndvi) # check first 6 lines## Class DOY96 DOY112 DOY130 DOY144 DOY159 DOY192 DOY208 DOY228 DOY256 DOY270
## 21 cotton 0.1278 0.0668 0.0553 0.0972 0.1220 0.7251 0.7916 0.7432 0.6334 0.4983
## 22 cotton 0.1278 0.0665 0.0553 0.0868 0.1218 0.7027 0.7585 0.7323 0.6186 0.4871
## 27 cotton 0.1157 0.0841 0.1234 0.1253 0.1519 0.2402 0.2969 0.3398 0.2512 0.3306
## 47 cotton 0.1402 0.0731 0.0370 0.1084 0.1758 0.7599 0.8256 0.7902 0.6509 0.5038
## 53 cotton 0.1108 0.0537 0.0510 0.0972 0.1455 0.6849 0.8143 0.8284 0.6767 0.5568
## 55 cotton 0.1426 0.0731 0.0379 0.1085 0.1814 0.7408 0.8140 0.8133 0.6240 0.4971
## DOY289
## 21 0.3435
## 22 0.3374
## 27 0.1993
## 47 0.3319
## 53 0.3317
## 55 0.3590次に、主成分分析を行って、データのもつ変異を要約してみましょう。読み込んだデータから、NDVIデータ（2列目から最後の列まで）を変数Xとして抽出します（Xは行列）。また、作物種を変数clとして抽出します。
X <- as.matrix(ndvi[, 2:ncol(ndvi)]) # Extract wavelength data
cl <- ndvi$Class
res.pca <- prcomp(X) # principal component analysis次に、主成分分析の結果を描いてみよう。
pairs(res.pca$x[,1:3], col = cl)
plotly パッケージの plot_ly 関数を使用して 3 次元散布図を描画します。
plotly::plot_ly(x = res.pca$x[,1], y = res.pca$x[,2],
z = res.pca$x[,3], color = cl,
type = "scatter3d", mode = "markers", size = 0.1)主成分得点をもとに散布図を描くと、第1、第2主成分のグラフでは、スイカが他の作物種から離れた場所に位置していることがわかる。一方、ブドウ、ワタ、トウモロコシは、第1〜3主成分の空間で互いに重なり合っていることがわかります。
次に、線形判別分析を用いて分類モデルを構築してみる。線形判別分析には、MASSパッケージに含まれるlda関数を使用します。lda関数は、どの変数で何を判別するかを指定する必要があります。ここでは、“Class ~ .”として指定します。これは、ndviというデータの残りの変数でClass（因子でなければならない）を識別するモデルを作成することを意味します。
res.lda <- MASS::lda(Class ~ . , data = ndvi) # linear discriminant analysis
res.lda # Display results## Call:
## lda(Class ~ ., data = ndvi)
##
## Prior probabilities of groups:
## cotton grape maize watermelon
## 0.586 0.170 0.065 0.179
##
## Group means:
## DOY96 DOY112 DOY130 DOY144 DOY159 DOY192
## cotton 0.09183481 0.07079317 0.06931945 0.1138387 0.2074708 0.7767010
## grape 0.10465529 0.08929941 0.07177353 0.1158300 0.1834359 0.5488500
## maize 0.10179385 0.07642615 0.07822000 0.1518662 0.3258077 0.8426938
## watermelon 0.09710670 0.09306425 0.07633296 0.1052223 0.1310408 0.6249503
## DOY208 DOY228 DOY256 DOY270 DOY289
## cotton 0.8182379 0.8005734 0.6696370 0.5292558 0.3142102
## grape 0.6512294 0.6459153 0.5694894 0.5027753 0.2954747
## maize 0.8459585 0.8108923 0.5997385 0.4333092 0.2346062
## watermelon 0.7282106 0.5984654 0.3159866 0.1956598 0.1454872
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3
## DOY96 1.239276 8.9295380 -14.2410016
## DOY112 22.052720 -9.6572911 -9.4822330
## DOY130 3.153688 9.3666079 12.9301262
## DOY144 -3.195624 -0.7050797 13.0531728
## DOY159 -1.873597 -2.9388017 -24.1224798
## DOY192 -6.897472 8.9389458 0.9706753
## DOY208 19.383390 -5.0019306 6.9063329
## DOY228 -3.897492 10.8862230 -1.1247100
## DOY256 -11.118489 -2.4111915 -3.3169339
## DOY270 -2.130809 -3.4693096 3.7571186
## DOY289 -4.117050 -4.0018919 6.1395059
##
## Proportion of trace:
## LD1 LD2 LD3
## 0.7801 0.1437 0.0763解析結果は、各サンプルが各クラス（グループ）に属する確率（事前確率）、各クラスのNDVI値の平均値、線形判別関数の係数、そして各線形判別関数がグループ間分散のどの程度を説明しているかを示しています。具体的には、線形判別関数の係数と各線形判別関数で説明される群間分散の割合が表示されます。
次に、判別得点を計算します。
fit.lda <- predict(res.lda) # calculate discriminant scores
head(fit.lda$x)## LD1 LD2 LD3
## 21 -0.7454548 0.01170553 1.54443730
## 22 -0.9495328 -0.03185156 1.14782997
## 27 0.2183818 -3.90433866 -1.60240434
## 47 -0.7081727 0.36198457 0.02702861
## 53 -1.2852901 0.07155552 1.31392312
## 55 -0.6945178 0.49373205 -0.02349397主成分得点と同様に、散布図を描いてみます。
pairs(fit.lda$x[,1:3], col = cl)
3次元で散布図を描いてみましょう。
plotly::plot_ly(x = fit.lda$x[,1], y = fit.lda$x[,2],
z = fit.lda$x[,3], color = cl,
type = "scatter3d", mode = "markers", size = 0.1)なお、軸（判別関数）はグループをうまく分離するように配置されているので、主成分分析の場合よりも、グループ間の重なりがはっきりしています。例えば、トウモロコシとワタ、ブドウとワタなどは、LD2やLD3を使って何とか分離できそうです。
構築したモデルを用いて、分類精度を確認してみましょう。
cl.est <- fit.lda$class # Estimation class for all samples
tab <- table(cl, cl.est) # assignment to tab
tab## cl.est
## cl cotton grape maize watermelon
## cotton 544 39 2 1
## grape 53 106 3 8
## maize 19 0 46 0
## watermelon 0 3 0 176diag(tab) # Extract only the diagonal components## cotton grape maize watermelon
## 544 106 46 176sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.872各サンプルについて同定（推定）されたクラスは，predict 関数を用いて取得することができます．真のクラス（cl）と同定（推定）クラス（cl.est）の分割表を作成し，比較することにより，正しく（あるいは誤って）同定された標本数を確認することができます．表の対角成分は正しく同定されたケースなので、diag関数で対角成分のみを取り出し、それらを合計し、総サンプル数で割ることで、正しく同定された割合を算出します。正答率は87.2%であり、特にブドウとワタを間違えるケースが多い（170サンプル中53サンプル）ことがわかります。
なお、ここで算出された精度は、あくまで訓練データ（モデル構築に使用したデータ）にモデルを「当てはめた」ものです。このような場合、前回の講義と同様に「過剰適合」（over-fitting）が起こっている可能性があります。
構築するモデルが解析されたデータを記述するためのものであれば、訓練データで分類精度を評価しても構わないかもしれません。しかし、多くの場合、そのようなモデルを構築する目的は、新しいデータのクラスを予測することです。そのような場合には、交差検証などの手法により、モデル構築に使用しなかったデータ（以下、「検証データ」と呼ぶ）の分類精度を評価する必要があります。 ここでは、Hao et al. (2015)のデータから、モデル評価用の検証データを別途用意した。検証データを読み込んでみましょう。
ndvi.test <- read.csv(here("data", "ndvi_test.csv"),
row.names = 1, stringsAsFactors = T) # Read data
head(ndvi.test) # check first 6 lines## Class DOY96 DOY112 DOY130 DOY144 DOY159 DOY192 DOY208 DOY228 DOY256 DOY270
## 1 cotton 0.1305 0.0731 0.0637 0.1143 0.2143 0.7696 0.8281 0.7887 0.6716 0.4999
## 3 cotton 0.0860 0.0844 0.0851 0.1421 0.1327 0.3117 0.3341 0.3971 0.2237 0.3605
## 6 cotton 0.1216 0.0678 0.0646 0.1027 0.2037 0.7829 0.8274 0.7668 0.7010 0.5201
## 8 cotton 0.1134 0.0680 0.0528 0.1138 0.1762 0.6998 0.7722 0.7633 0.6218 0.5164
## 14 cotton 0.1352 0.0731 0.0553 0.0975 0.1254 0.7062 0.7467 0.7514 0.6363 0.4876
## 24 cotton 0.1459 0.0668 0.0532 0.0919 0.1286 0.7096 0.7538 0.6604 0.5642 0.4306
## DOY289
## 1 0.3400
## 3 0.2358
## 6 0.3285
## 8 0.3274
## 14 0.3491
## 24 0.2938cl.test <- ndvi.test$Class # Extract class (crop species) data次に、構築した識別モデルres.ldaを用いて、評価用データの識別を行ってみましょう。
pred.lda <- predict(res.lda, ndvi.test)
cl.pred <- pred.lda$class # Predict to unknown data
head(cl.pred, 10)## [1] cotton grape cotton cotton cotton cotton grape cotton cotton cotton
## Levels: cotton grape maize watermelonhead(cl.test, 10)## [1] cotton cotton cotton cotton cotton cotton cotton cotton cotton cotton
## Levels: cotton grape maize watermelon識別された（予測された）クラスと実際のクラスを分割表を使って比較してみましょう。先ほどと同じ手順で分類の正答率も計算してみます。
table(cl.test, cl.pred) # Create binary division table## cl.pred
## cl.test cotton grape maize watermelon
## cotton 265 24 1 0
## grape 23 56 3 4
## maize 15 0 22 0
## watermelon 0 0 0 87tab <- table(cl.test, cl.pred) # assignment to tab
sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.86ワタとぶどうの区別、ワタとトウモロコシの区別が難しいことがわかります。正解率は86%で、訓練データで計算した正解率より若干低くなっています。
検証データ内のサンプルについて、判別スコアを計算して、プロットしてみましょう。
pairs(pred.lda$x[, 1:3], col = cl.test,
pch = as.numeric(cl.pred))
色の違いは真のクラス、記号の違いは予測されたクラスを表しています。誤って分類されたサンプルは、クラス間の境界に位置することが多いことが分かります。
サポートベクターマシン(SVM)
線形判別分析は、元の変数Xの線形結合を用いてクラス識別を行う手法です。しかし、この方法は、クラス間の境界が線形分離可能（直線で分離できる）でなければ、あまり正確な結果を得ることができません。こうした場合、サポートベクターマシン（SVM）を使用することが有効です。SVMは、クラス間の境界が直線的に分離できない複雑なパターンの場合に、カーネル法という方法を用いて、より適切なモデルを構築することができます。
では、実際にSVMを使用してその精度を確認してみましょう。解析には、kernlabパッケージに含まれるksvmという関数を使用します。SVMの種類は、サポートベクターマシンを使用する目的に応じて選択する必要があります。クラス識別モデルを構築する場合は、まず「C-svc」を選択すると良いでしょう。SVMは多変量回帰や外れ値検出にも使用することができ、これらの場合にはそれぞれ「eps-svr」や「one-svc」を選択します。このほかにも、さまざまなタイプを指定することが可能です。
res.svm <- kernlab::ksvm(Class ~ . ,
data = ndvi, type = "C-svc") # Execute SVM
res.svm## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc (classification)
## parameter : cost C = 1
##
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 0.136294283540592
##
## Number of Support Vectors : 334
##
## Objective Function Value : -113.6496 -54.0899 -20.9744 -21.7727 -19.1527 -7.8739
## Training error : 0.038では、特定（推定）されたクラスと実際のクラスを比較してみましょう。
cl.est <- predict(res.svm) # Estimation class for all samples
tab <- table(cl, cl.est) # assignment to tab
tab## cl.est
## cl cotton grape maize watermelon
## cotton 580 6 0 0
## grape 17 152 1 0
## maize 13 0 52 0
## watermelon 0 1 0 178sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.962線形判別分析と比較して、分類の正答率（96.2％）が向上しています。また、ブドウとワタを間違える率も低くなっています（170サンプル中17サンプル）。
この訓練データを用いて、10分割交差検証によりモデルの分類精度を評価します。
res.svm <- kernlab::ksvm(Class ~ . , data = ndvi,
type = "C-svc",
cross = 10) # Execute SVM with cross-validation
res.svm## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc (classification)
## parameter : cost C = 1
##
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 0.137295216424973
##
## Number of Support Vectors : 334
##
## Objective Function Value : -113.4212 -54.001 -20.9704 -21.7312 -19.1381 -7.8824
## Training error : 0.038
## Cross validation error : 0.056cross(res.svm)## [1] 0.0561 - cross(res.svm)## [1] 0.944では、テストデータに対して、識別の精度を確認してみましょう。
cl.pred <- predict(res.svm, ndvi.test) # Predict to unknown data
tab <- table(cl.test, cl.pred) # assignment to tab
tab## cl.pred
## cl.test cotton grape maize watermelon
## cotton 283 7 0 0
## grape 14 70 0 2
## maize 13 1 23 0
## watermelon 0 0 0 87sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.926また、テストデータによる検証でも、高い分類精度（93％）を示しています。
SVMでは、様々なカーネル行列（サンプル間の関係を定義する行列。具体的には、特徴空間におけるサンプル間の内積行列である）を用いて分類モデルを作成することができます。
例えば、線形カーネルは、以下のように適用できます。
res.svm.lin <- kernlab::ksvm(Class ~ . , data = ndvi,
type = "C-svc",
kernel = "vanilladot", # linear kernel
cross = 10) ## Setting default kernel parametersres.svm.lin## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc (classification)
## parameter : cost C = 1
##
## Linear (vanilla) kernel function.
##
## Number of Support Vectors : 269
##
## Objective Function Value : -175.02 -52.5992 -10.4078 -12.954 -12.7449 -0.7291
## Training error : 0.091
## Cross validation error : 0.106線形カーネルで精度を確認します。
cl.pred.lin <- predict(res.svm.lin, ndvi.test) # Estimation class for all samples
tab <- table(cl.test, cl.pred.lin) # assignment to tab
tab## cl.pred.lin
## cl.test cotton grape maize watermelon
## cotton 266 23 1 0
## grape 13 71 1 1
## maize 13 2 22 0
## watermelon 0 0 0 87sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.892デフォルトのカーネル（放射基底関数カーネル、Radial basis function: RBF、rbfdot）よりも精度が低くなっています。
次に、3次多項式カーネルを使用してみます。
res.svm.pol <- kernlab::ksvm(Class ~ . , data = ndvi,
type = "C-svc",
kernel = "polydot",
kpar = list(degree = 3), # polynomial kernel
cross = 10)
res.svm.pol## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc (classification)
## parameter : cost C = 1
##
## Polynomial kernel function.
## Hyperparameters : degree = 3 scale = 1 offset = 1
##
## Number of Support Vectors : 189
##
## Objective Function Value : -2.5758 -2.1475 -0.0294 -0.0911 -0.0293 -0.0161
## Training error : 0
## Cross validation error : 0.084多項式カーネルで精度を確認します。
cl.pred.pol <- predict(res.svm.pol, ndvi.test) # Estimation class for all samples
tab <- table(cl.test, cl.pred.pol) # assignment to tab
tab## cl.pred.pol
## cl.test cotton grape maize watermelon
## cotton 269 15 6 0
## grape 9 77 0 0
## maize 5 1 31 0
## watermelon 0 0 0 87sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 0.928精度はデフォルトのカーネル（放射基底関数カーネル rbfdot）と同等ですが、2つのカーネルでは正解・不正解のパターンが異なっています。
table(cl.pred, cl.pred.pol, cl.test)## , , cl.test = cotton
##
## cl.pred.pol
## cl.pred cotton grape maize watermelon
## cotton 265 12 6 0
## grape 4 3 0 0
## maize 0 0 0 0
## watermelon 0 0 0 0
##
## , , cl.test = grape
##
## cl.pred.pol
## cl.pred cotton grape maize watermelon
## cotton 5 9 0 0
## grape 4 66 0 0
## maize 0 0 0 0
## watermelon 0 2 0 0
##
## , , cl.test = maize
##
## cl.pred.pol
## cl.pred cotton grape maize watermelon
## cotton 4 0 9 0
## grape 0 1 0 0
## maize 1 0 22 0
## watermelon 0 0 0 0
##
## , , cl.test = watermelon
##
## cl.pred.pol
## cl.pred cotton grape maize watermelon
## cotton 0 0 0 0
## grape 0 0 0 0
## maize 0 0 0 0
## watermelon 0 0 0 87Random Forst
機械学習手法ランダムフォレストを使って分類モデルを構築してみましょう。ランダムフォレストは、識別に用いる変数とクラスの間に非線形な関係がある場合にも用いることができます。
ランダムフォレストでは、rangerパッケージに含まれるranger関数を使用します。
res.rfo <- ranger::ranger(Class ~ . , data = ndvi) # Run RF
# Estimated class vs. actual
fit.rfo <- predict(res.rfo, ndvi) # Estimation class for all samples予測結果は次のようにして取り出すことができます。
cl.est <- fit.rfo$predictions
tab <- table(cl, cl.est) # assignment to tab
tab## cl.est
## cl cotton grape maize watermelon
## cotton 586 0 0 0
## grape 0 170 0 0
## maize 0 0 65 0
## watermelon 0 0 0 179sum(diag(tab)) / sum(tab) # rate of correct prediction## [1] 1このモデルは学習データに完全に適合しています。
モデルの精度は、OOB（out of bag）誤差で評価できます。
accuracy <- 1 - res.rfo$prediction.error
accuracy## [1] 0.935検証データにモデルを当てはめ、予測精度を確認します。
pred.rfo <- predict(res.rfo, ndvi.test) # Prediction of the classes for all samples
cl.pred <- pred.rfo$predictions
tab <- table(cl.test, cl.pred) # assignment to tab
tab## cl.pred
## cl.test cotton grape maize watermelon
## cotton 282 7 1 0
## grape 17 67 0 2
## maize 10 0 27 0
## watermelon 0 0 0 87sum(diag(tab)) / sum(tab) # rate of correct prediction## [1] 0.926テストデータに対する予測では、SVMと同程度の正答率で分類されていることがわかります。SVMと比較すると、特にトウモロコシの識別がより正確であるように見えます。
RFでは、変数の重要度を確認することができます。importance オプションを “none” から “impurity” に設定することで、識別における変数の重要度が算出されます。
res.rfo.wim <- ranger::ranger(Class ~ . , data = ndvi,
importance = "impurity")
var.imp <- ranger::importance(res.rfo.wim)
var.imp## DOY96 DOY112 DOY130 DOY144 DOY159 DOY192 DOY208 DOY228
## 26.52716 34.89458 19.88409 42.63709 61.51464 53.34131 45.55458 69.21589
## DOY256 DOY270 DOY289
## 100.17419 82.18946 54.16634barplot(var.imp, las = 2)
夏以降（DOY228以降）のデータが分類・予測に重要であることがわかります。
ブルーベリーデータへの適用
では、Gilbertら（2015, PLoS ONE 10:e0138494）のデータから5品種を抜き出したものを使用して、ブルーベリーの揮発性成分から品種を予測するモデルを構築してみましょう。
まず、モデル構築のためのデータをロードします。
bb <- read.csv(here("data", "blueberry_5var_all.csv"), stringsAsFactors = T) # load data
head(bb) # check the first six lines## gID Year Location Harvest Genotype Panelists Overall.Liking
## 1 2012 C1 Emerald 2012 C 1 Emerald 95 23.5
## 2 2012 C1 Endura 2012 C 1 Endura 95 17.5
## 3 2012 C1 Farthing 2012 C 1 Farthing 95 24.2
## 4 2012 C1 Meadowlark 2012 C 1 Meadowlark 95 22.3
## 5 2012 C1 Primadonna 2012 C 1 Primadonna 95 24.0
## 6 2012 C2 Emerald 2012 C 2 Emerald 90 20.0
## Texture Sweetness Sourness Flavor SSC TA SSC.TA pH Fructose Glucose
## 1 24.9 22.4 10.7 25.0 10.4 0.5110 20.4 3.43 48.74 49.39
## 2 25.1 19.2 22.6 27.8 10.5 0.6199 16.9 3.36 47.79 46.49
## 3 23.9 23.6 11.8 26.5 10.6 0.3064 34.6 3.55 51.31 51.66
## 4 28.7 21.5 11.2 24.4 9.2 0.2631 34.8 3.78 42.57 45.11
## 5 23.6 25.5 10.4 25.8 11.3 0.2346 48.2 3.95 56.83 56.35
## 6 23.0 20.2 13.6 22.3 10.6 0.3992 26.6 3.54 48.12 49.15
## Sucrose Total.Sugar X590.86.3 X616.25.1 X107.87.9 X110.62.3 X105.37.3
## 1 0.13 98.26 1.05 7.20 3.79 4.79 0
## 2 0.19 94.47 1.40 8.31 3.56 4.32 0
## 3 0.18 103.15 0.82 5.27 21.69 4.53 0
## 4 0.12 87.79 0.85 5.25 5.87 4.20 0
## 5 0.21 113.39 1.82 7.33 6.88 7.47 0
## 6 0.22 97.48 0.77 5.06 3.69 5.95 0
## X623.42.7 X123.51.3 X1576.87.0 X71.41.0 X1576.95.0 X108.88.3 X556.24.1
## 1 0.11 0.33 4.46 0.94 1.92 0.03 0.09
## 2 0.43 0.55 3.73 0.70 2.32 0.08 8.70
## 3 0.13 0.59 3.05 0.87 1.21 0.00 7.65
## 4 0.20 0.24 2.12 1.23 1.52 0.02 43.87
## 5 0.36 0.50 4.93 1.53 1.67 0.16 216.42
## 6 0.08 0.20 3.54 2.39 1.27 0.01 0.15
## X66.25.1 X6728.26.3 X928.95.0 X111.27.3 X123.92.2 X110.43.0 X111.71.7
## 1 1196.56 4345.51 355.57 190.77 0 0.72 0.90
## 2 907.13 5165.84 401.18 254.44 0 1.35 0.53
## 3 1058.29 3146.05 358.01 166.66 0 1.61 0.35
## 4 1380.10 4123.45 390.25 245.92 0 2.62 0.49
## 5 1766.28 4999.20 290.46 226.99 0 0.78 1.02
## 6 1166.00 3426.35 172.85 121.11 0 1.51 0.75
## X1191.16.8 X106.70.7 X7785.70.8 X928.68.7 X142.62.1 X110.93.0 X111.13.7
## 1 0.00 0.08 0.00 0.32 0.42 2.42 0.01
## 2 0.22 0.26 0.05 0.43 0.33 2.54 0.16
## 3 0.28 0.05 0.00 0.07 0.37 1.53 0.00
## 4 0.00 0.09 0.00 0.46 0.33 12.18 0.15
## 5 0.00 0.25 0.00 0.31 0.65 5.21 0.12
## 6 0.00 0.05 0.00 0.36 0.23 1.76 0.00
## X123.66.0 X3681.71.8 X142.92.7 X2497.18.9 X104.76.7 X5989.27.5 X470.82.6
## 1 1.54 0 4.27 0.00 0 0.79 1.87
## 2 1.42 0 10.75 20.05 0 4.94 10.63
## 3 0.04 0 1.92 5.78 0 0.43 3.29
## 4 3.00 0 6.63 0.00 0 4.22 1.26
## 5 1.96 0 7.46 0.00 0 2.51 1.92
## 6 0.64 0 0.00 1.55 0 0.37 0.72
## X122.78.1 X821.55.6 X78.70.6 X124.19.6 X2639.63.6 X53398.83.7 X112.40.3
## 1 0.13 0.50 36.42 0.74 0.21 0.71 0.42
## 2 0.42 4.79 16.18 0.69 0.18 3.38 0.05
## 3 0.07 0.32 10.86 0.26 0.29 1.41 0.01
## 4 0.20 1.62 21.46 0.25 0.18 0.46 0.25
## 5 0.07 0.34 10.41 0.89 0.20 0.69 0.17
## 6 0.05 0.24 18.33 0.03 0.04 0.18 0.31
## X119.36.8 X106.26.3 X141.27.5 X112.12.9 X629.50.5 X3879.26.3 X689.67.8
## 1 1.42 0.70 1.63 0.38 0.11 0.00 2.32
## 2 0.38 0.61 1.94 1.39 0.11 0.14 3.30
## 3 0.11 0.39 0.52 0.11 0.05 0.00 0.66
## 4 0.00 0.18 0.03 0.85 0.09 0.42 14.76
## 5 0.00 0.55 1.67 0.05 0.06 0.00 1.52
## 6 0.29 0.40 0.72 0.27 0.05 0.08 1.77
## X1139.30.6 X5989.33.3 X43219.68.7 X582.16.1
## 1 1.78 0.14 1.27 0.00
## 2 5.46 1.65 0.66 0.31
## 3 2.08 0.15 0.20 0.22
## 4 4.74 0.81 1.50 0.12
## 5 1.65 0.83 0.41 0.17
## 6 1.17 0.07 0.37 0.03まず、主成分分析を行って、データのもつ変異を要約してみます。読み込んだデータから、揮発性成分のデータ（20列目から最後の列まで）を変数xとして抽出します（xは行列）。また、品種名をculとして抽出します。
x <- as.matrix(bb[, 20:ncol(bb)]) # volatile components
x <- scale(x) # scale the variables
cl <- bb$Genotype # variety name
res.pca <- prcomp(x)plotly パッケージの plot_ly 関数を使用して 3 次元散布図を描画します。
plotly::plot_ly(x = res.pca$x[,1], y = res.pca$x[,2],
z = res.pca$x[,3], color = cl,
type = "scatter3d", mode = "markers", size = 0.1)主成分得点をもとに散布図を描くと、第1主成分（x）、第2主成分（y）のグラフでは、Meadowlarkが他の品種から少しわかれて位置していることがわかります。一方、Endura、Farthing、Primadonnaの3品種は、主に3主成分（z）で互いに分かれてくることがわかります。また、これら3品種とEmeraldは、Emeraldの分布が偏っているものの、互いに重なり合っていることがわかります。
次に、線形判別分析を用いて分類モデルを構築してみましょう。まずは、解析のためのデータを準備してみましょう。
bb.vc <- data.frame(y = cl, x)次に、線形判別分析をおこなってみます。
res.lda <- MASS::lda(y ~ . , data = bb.vc) # linear discriminant analysis判別得点を計算し、散布図を描いてみます。
fit.lda <- predict(res.lda) # calculate discriminant scores
pairs(fit.lda$x[,1:3], col = cl)
非常にきれいに分離されることがわかります。
次に、3次元で散布図を描きます。
plotly::plot_ly(x = fit.lda$x[,1], y = fit.lda$x[,2],
z = fit.lda$x[,3], color = cl,
type = "scatter3d", mode = "markers", size = 0.1)なお、軸（判別関数）はグループをうまく分離するように配置されており、グループ間がよく分離していることがわかります。
最後に、線形判別モデルの分類精度を確認してみましょう。
cl.est <- fit.lda$class # Estimation class for all samples
tab <- table(cl, cl.est) # assignment to tab
tab## cl.est
## cl Emerald Endura Farthing Meadowlark Primadonna
## Emerald 27 0 0 0 0
## Endura 0 27 0 0 0
## Farthing 0 0 27 0 0
## Meadowlark 0 0 0 27 0
## Primadonna 0 0 0 0 25diag(tab) # Extract only the diagonal components## Emerald Endura Farthing Meadowlark Primadonna
## 27 27 27 27 25sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 1分類の正解率は100%です。もちろん、これはデータに当てはめただけで、予測精度ではありません。
最後に、1個抜き交差検証を行って、精度を確認してみましょう。
cl.pred <- rep(NA, length(cl))
for(i in 1:nrow(bb.vc)) {
bb.train <- bb.vc[-i, ]
res.lda <- MASS::lda(y ~ . , data = bb.train)
fit.lda <- predict(res.lda, bb.vc[i, ])
cl.pred[i] <- fit.lda$class
}tab <- table(cl, cl.pred) # assignment to tab
tab## cl.pred
## cl 1 2 3 4 5
## Emerald 27 0 0 0 0
## Endura 0 27 0 0 0
## Farthing 0 0 27 0 0
## Meadowlark 0 0 0 27 0
## Primadonna 0 0 0 0 25diag(tab) # Extract only the diagonal components## [1] 27 27 27 27 25sum(diag(tab)) / sum(tab) # Calculate percentage of correct answers## [1] 11個抜き交差検証でも100％の精度が得られていることから、データによっては線形判別分析でも非常に正確に予測できる場合があることがわかります。どの手法が適しているかはデータに大きく依存し、例えば、非線形の方法が常に線形の方法よりも精度が高いとは限りません。このように、すべての手法があらゆる問題に対して同様に優れているわけではないことを「No Free Lunch」の定理と呼びます。万能な手法は存在せず、実際の問題解決においては、データの特性や問題の構造に基づいて最適な手法を選定することが重要です。