機器學習概論簡介

本週開始教導機器學習的部分，與生物統計學的數學模型有所不同了，這邊每一個演算法都可能有自己的運算邏輯！
至於哪個演算法比較好呢? 依據不同的資料型態，這可能很難說，原則上決定的方式還是以「資料科學實驗」的流程為主，在一個「驗證組」中進行模型的選擇，並在「測試組」。
我們這裡介紹一個機器學習界常講的說法：無免費午餐定理(No Free Lunch Theorem)，主要在說對於各種問題都沒有單一的最佳優化算法。
除了上週我們開始接觸到「調參數」之外，我們還需要學習更多的機器學習算法，從而在你的任務中找到最佳的預測模型。

– 以後你會發現在一個特定任務的「資料科學實驗」中，你可能會使用10種以上的算法，而每種算法至少有100組參數在比較，從而在驗證組中找到最佳的模型最後再應用至測試組。

第一節：線性判別分析(1)

線性判別分析(Linear Discriminant Analysis)是一種監督學習的降維技術，他的運算過程是試圖找出一個超平面，從而將不同屬性的樣本分離：

F01

他的思想可以用一句話概括，就是「找到一個投影方式，讓投影後的組內差異最小化，組間差異最大化」。

– 請至這裡下載範例資料

dat <- read.csv("ECG_train.csv", header = TRUE, fileEncoding = 'CP950', stringsAsFactors = FALSE, na.strings = "")

第一節：線性判別分析(2)

「投影」是線性代數裡面的一個詞彙，簡單來說是說矩陣的乘法可以幫助我們把「原始資料」「投影」到新的「座標系」上，我們來簡單做一個小示範。

– 這是一組資料，我們可以透過下面這個「旋轉投影矩陣\(R\)」，強迫二維座標點轉軸\(\theta\)度

\[R = \begin{pmatrix} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \end{pmatrix}\]

X_mat <- cbind(x = c(1.7, 2.5, -2.1), y = c(3.1, -2.8, -0.7))

rotate_degree <- 0.5
R_mat <- matrix(c(cos(rotate_degree), -sin(rotate_degree), sin(rotate_degree), cos(rotate_degree)), nrow = 2)
XR_mat <- X_mat %*% R_mat

par(mfrow = c(1, 2))
plot(X_mat[,1], X_mat[,2], xlim = c(-5, 5), ylim = c(-5, 5), pch = 19, cex = 2, col = 2:4, main = 'before rotate')
plot(XR_mat[,1], XR_mat[,2], xlim = c(-5, 5), ylim = c(-5, 5), pch = 19, cex = 2, col = 2:4, main = 'after rotate')

– 你可以換一下【rotate_degree】看看效果

第一節：線性判別分析(3)

有這樣的了解後，我們來試著直接使用套件「MASS」中的函數「lda()」進行投影變換。

– 我們的目標是找到一個矩陣\(R\)，讓我們能用「Rate」和「QTc」預測「AMI」：

library(MASS)

dat[,'AMI'] <- factor(dat[,'AMI'])
lda_fit <- lda(AMI ~ Rate + QTc, data = dat)
lda_fit$scaling

##              LD1         LD2
## Rate -0.02371907  0.03725259
## QTc  -0.01476256 -0.01480724

– 找到矩陣後我們來圖像化表達一下線性判別分析的效果：

X_mat <- cbind(dat[,'Rate'], dat[,'QTc'])
R_mat <- lda_fit$scaling
XR_mat <- X_mat %*% R_mat
col_idx <- as.integer(dat[,'AMI'])
col_list <- c('#FF000030', '#00A00030', '#0000FF30')

par(mfrow = c(1, 2))
plot(X_mat[,1], X_mat[,2], pch = 19, cex = 0.5, col = col_list[col_idx], main = 'before rotate')
plot(XR_mat[,1], XR_mat[,2], pch = 19, cex = 0.5, col = col_list[col_idx], main = 'after rotate')

沒什麼用對吧，畢竟這是線性的轉換，其實他跟softmat regression沒什麼差別。

第一節：線性判別分析(4)

如果我們直接取得他預測的class，你會發現結果怪怪的：

result_list <- predict(lda_fit, dat)
table(dat[,'AMI'], result_list$class)

##          
##           not-AMI NSTEMI STEMI
##   not-AMI    1055      0     1
##   NSTEMI      168      0     0
##   STEMI       255      0     0

原因又是一樣的，就是我們這個樣本其實沒有很好的分離能力，而無論結果怎麼樣都最有可能是盛行率最大的那個值。

– 這是對不同類別平均數的距離：

head(result_list$posterior)

##     not-AMI     NSTEMI     STEMI
## 1 0.5799250 0.14117064 0.2789043
## 2 0.8068244 0.07610796 0.1170676
## 3 0.6695866 0.13224111 0.1981723
## 4 0.6559637 0.14644156 0.1975948
## 5 0.7454229 0.09747472 0.1571024
## 6 0.7566320 0.09973263 0.1436353

我們來使用ROC曲線與原始資料做出比較：

library(pROC)

roc_curve <- roc(response = (dat[,'AMI'] == 'STEMI'), predictor = result_list$posterior[,3])
plot(roc_curve)

練習1：與迴歸分析做出比較

讓我們使用「Rate」和「QTc」預測「LVD」，我們一樣先將樣本分成「訓練組」、「驗證組」與「測試組」。

subdat <- dat[!(dat[,"LVD"] %in% NA) & !(dat[,"Rate"] %in% NA) & !(dat[,"QTc"] %in% NA),]
X <- subdat[,c('Rate', 'QTc')]
Y <- subdat[,"LVD"]

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- X[train_idx,]
valid_X <- X[valid_idx,]
test_X <- X[test_idx,]

train_Y <- Y[train_idx]
valid_Y <- Y[valid_idx]
test_Y <- Y[test_idx]

請你比較線性判別分析與邏輯斯回歸的差異！

練習1答案(1)

這是邏輯斯回歸的結果：

glm_fit <- glm(train_Y ~ ., data = train_X, family = 'binomial')

valid_pred <- predict(glm_fit, valid_X)
test_pred <- predict(glm_fit, test_X)

library(pROC)

roc_valid <- roc(valid_Y ~ valid_pred)
best_pos <- which.max(roc_valid$sensitivities + roc_valid$specificities)
best_cut <- roc_valid$thresholds[best_pos]

tab_test <- table(test_pred >= best_cut, test_Y)
sens <- tab_test[2,2] / sum(tab_test[,2])
spec <- tab_test[1,1] / sum(tab_test[,1])

roc_test <- roc(test_Y ~ test_pred)
plot(roc_test)

points(spec, sens, pch = 19)
text(0.5, 0.5, paste0('AUC(valid) = ', formatC(roc_valid$auc, digits = 3, format = 'f'),
                      '\nSens = ', formatC(sens, digits = 3, format = 'f'),
                      '\nSpec = ', formatC(spec, digits = 3, format = 'f'),
                      '\nAUC(test) = ', formatC(roc_test$auc, digits = 3, format = 'f')), col = 'red')

練習1答案(2)

這是線性判別分析的結果：

lda_fit <- lda(train_Y ~ ., data = train_X)

valid_list <- predict(lda_fit, valid_X)
test_list <- predict(lda_fit, test_X)

valid_pred <- valid_list$posterior[,2]
test_pred <- test_list$posterior[,2]

library(pROC)

roc_valid <- roc(valid_Y ~ valid_pred)
best_pos <- which.max(roc_valid$sensitivities + roc_valid$specificities)
best_cut <- roc_valid$thresholds[best_pos]

tab_test <- table(test_pred >= best_cut, test_Y)
sens <- tab_test[2,2] / sum(tab_test[,2])
spec <- tab_test[1,1] / sum(tab_test[,1])

roc_test <- roc(test_Y ~ test_pred)
plot(roc_test)

points(spec, sens, pch = 19)
text(0.5, 0.5, paste0('AUC(valid) = ', formatC(roc_valid$auc, digits = 3, format = 'f'),
                      '\nSens = ', formatC(sens, digits = 3, format = 'f'),
                      '\nSpec = ', formatC(spec, digits = 3, format = 'f'),
                      '\nAUC(test) = ', formatC(roc_test$auc, digits = 3, format = 'f')), col = 'red')

結果完全一樣，所以拜託以後論文中不要同時比較邏輯斯回歸與線性判別分析！

第二節：樸素貝葉斯分類器(1)

樸素貝葉斯分類器(Naive Bayes Classification)是基於機率理論，他直接假設所有的隨機變數之間具獨立特性，因此可以直接利用條件機率相乘的方法，計算出聯合機率分布。

– 我們假設具有樣本\(X\)具有條件\(X = \begin{pmatrix} x_1 & x_2 & x_3 \end{pmatrix}\)，我們想要知道發生事件(\(Y\))的可能性\(L(Y|X)\)為何，那可以得到下列式子：

\[ \begin{align} L(Y|X) & = P(Y)P(x_1|Y)P(x_2|Y)P(x_3|Y) \end{align} \]

– 同樣的，我們也能計算出不發生事件(\(Y'\))的可能性\(L(Y'|X)\)：

\[ \begin{align} L(Y'|X) & = P(Y')P(x_1|Y')P(x_2|Y')P(x_3|Y') \end{align} \]

比較\(L(Y|X)\)與\(L(Y'|X)\)的大小，我們就能猜測是否發生事件。

– 那使用獨立假設有甚麼好處呢?好處在於我們應該有足夠多的樣本計算\(P(x_i|Y)\)，但如果條件足夠多，樣本中可能就不存在同時滿足\(X = \begin{pmatrix} x_1 & x_2 & x_3 \end{pmatrix}\)的樣本，從而無法做出準確的估計。

– 但獨立假設的缺點也很明顯，又是一個線性的預測工具，我們不可能透過樸素貝葉斯分類器找出交互作用存在。

第二節：樸素貝葉斯分類器(2)

還是有點難理解，我們用下面這個例子來進行解釋，我們想要預測樣本\(X\)具有下列特徵「女性、年齡介於31~40之間、上班族、月收入中等者」，會不會辦信用卡：

F02

\[ \begin{align} L(\text{辦卡:會}|X) & = P(\text{辦卡:會})P(\text{性別:女}|\text{辦卡:會})P(\text{年齡:31~40}|\text{辦卡:會})P(\text{職業:上班族}|\text{辦卡:會})P(\text{月收入:中}|\text{辦卡:會}) \\ & = \frac{6}{10} \times \frac{4}{6} \times \frac{2}{6} \times \frac{3}{6} \times \frac{1}{6} = 0.0111 \\\\ L(\text{辦卡:不會}|X) & = P(\text{辦卡:不會})P(\text{性別:女}|\text{辦卡:不會})P(\text{年齡:31~40}|\text{辦卡:不會})P(\text{職業:上班族}|\text{辦卡:不會})P(\text{月收入:中}|\text{辦卡:不會}) \\ & = \frac{4}{10} \times \frac{1}{4} \times \frac{3}{4} \times \frac{2}{4} \times \frac{1}{4} = 0.0188 \end{align} \]

很明顯因為\(L(\text{辦卡:不會}|X) > L(\text{辦卡:會}|X)\)，所以我們可以推論這個樣本\(X\)較有可能不會辦信用卡。

第二節：樸素貝葉斯分類器(3)

有這樣的了解後，我們來試著直接使用套件「e1071」中的函數「naiveBayes()」進行分析：

– 讓我們能用「GENDER」和「Rate」預測「AMI」：

library(e1071)

dat[,'GENDER'] <- factor(dat[,'GENDER'])
dat[,'AMI'] <- factor(dat[,'AMI'])

fit_Bayes <- naiveBayes(AMI ~ GENDER + Rate, data = dat)
fit_Bayes

## 
## Naive Bayes Classifier for Discrete Predictors
## 
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
## 
## A-priori probabilities:
## Y
##   not-AMI    NSTEMI     STEMI 
## 0.7147674 0.1132839 0.1719488 
## 
## Conditional probabilities:
##          GENDER
## Y            female      male
##   not-AMI 0.4915094 0.5084906
##   NSTEMI  0.1904762 0.8095238
##   STEMI   0.1411765 0.8588235
## 
##          Rate
## Y             [,1]     [,2]
##   not-AMI 90.12830 23.72952
##   NSTEMI  80.92262 22.54098
##   STEMI   83.04706 22.14905

裡面有一個變數是連續變項耶，不用緊張，我們可以用常態分布的機率密度函數求出機率。

第二節：樸素貝葉斯分類器(4)

我們來嘗試一下，計算當\(Rate = 100\)同時\(GENDER = male\)時，這個樣本的\(AMI\)是哪個比較有可能：

\[ \begin{align} L(\text{not-AMI}|\text{Rate = 100 & GENDER = male}) & = P(\text{not-AMI})P(\text{Rate = 100}|\text{not-AMI})P(\text{GENDER = male}|\text{not-AMI}) \\ L(\text{NSTEMI}|\text{Rate = 100 & GENDER = male}) & = P(\text{NSTEMI})P(\text{Rate = 100}|\text{NSTEMI})P(\text{GENDER = male}|\text{NSTEMI}) \\ L(\text{STEMI}|\text{Rate = 100 & GENDER = male}) & = P(\text{STEMI})P(\text{Rate = 100}|\text{STEMI})P(\text{GENDER = male}|\text{STEMI}) \end{align} \]

用R語言進行計算：

L1 <- prop.table(fit_Bayes$apriori)[1] * dnorm(100, mean = fit_Bayes$tables$Rate[1,1], sd = fit_Bayes$tables$Rate[1,2]) * fit_Bayes$tables$GENDER[1,2]
L2 <- prop.table(fit_Bayes$apriori)[2] * dnorm(100, mean = fit_Bayes$tables$Rate[2,1], sd = fit_Bayes$tables$Rate[2,2]) * fit_Bayes$tables$GENDER[2,2]
L3 <- prop.table(fit_Bayes$apriori)[3] * dnorm(100, mean = fit_Bayes$tables$Rate[3,1], sd = fit_Bayes$tables$Rate[3,2]) * fit_Bayes$tables$GENDER[3,2]
L_vec <- c(L1, L2, L3)
L_vec / sum(L_vec)

##   not-AMI    NSTEMI     STEMI 
## 0.6424388 0.1300580 0.2275032

看來最有可能的答案是\(AMI = not-AMI\)，我們來看看結果：

predict(fit_Bayes, newdata = data.frame(Rate = 100, GENDER = factor('male', levels = c('female', 'male'))), type = 'raw')

##        not-AMI   NSTEMI     STEMI
## [1,] 0.6424388 0.130058 0.2275032

練習2：比較樸素貝葉斯分類器與邏輯斯回歸

承練習1，讓我們再使用「Rate」和「QTc」預測「LVD」，我們一樣先將樣本分成「訓練組」、「驗證組」與「測試組」。

subdat <- dat[!(dat[,"LVD"] %in% NA) & !(dat[,"Rate"] %in% NA) & !(dat[,"QTc"] %in% NA),]
X <- subdat[,c('Rate', 'QTc')]
Y <- subdat[,"LVD"]

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- X[train_idx,]
valid_X <- X[valid_idx,]
test_X <- X[test_idx,]

train_Y <- Y[train_idx]
valid_Y <- Y[valid_idx]
test_Y <- Y[test_idx]

請你比較樸素貝葉斯分類器與邏輯斯回歸哪個更好！

練習2答案

這是樸素貝葉斯分類器的結果：

fit_Bayes <- naiveBayes(train_Y ~ ., data = train_X, family = 'binomial')

valid_pred <- predict(fit_Bayes, valid_X, type = 'raw')[,2]
test_pred <- predict(fit_Bayes, test_X, type = 'raw')[,2]

library(pROC)

roc_valid <- roc(valid_Y ~ valid_pred)
best_pos <- which.max(roc_valid$sensitivities + roc_valid$specificities)
best_cut <- roc_valid$thresholds[best_pos]

tab_test <- table(test_pred >= best_cut, test_Y)
sens <- tab_test[2,2] / sum(tab_test[,2])
spec <- tab_test[1,1] / sum(tab_test[,1])

roc_test <- roc(test_Y ~ test_pred)
plot(roc_test)

points(spec, sens, pch = 19)
text(0.5, 0.5, paste0('AUC(valid) = ', formatC(roc_valid$auc, digits = 3, format = 'f'),
                      '\nSens = ', formatC(sens, digits = 3, format = 'f'),
                      '\nSpec = ', formatC(spec, digits = 3, format = 'f'),
                      '\nAUC(test) = ', formatC(roc_test$auc, digits = 3, format = 'f')), col = 'red')

在驗證組中有好一點點，因此應該會選到樸素貝葉斯分類器，並且測試組中AUC也確實較高。

– 不過這不是重點，重點是結果終於不一樣了！論文中同時出現樸素貝葉斯分類器與邏輯斯回歸是可以接受的。

第三節：k-近鄰演算法(1)

有沒有覺得機器學習的模型還滿簡單的?統計方法好像還比較難?

– 我們再來教個更簡單的：k-近鄰演算法(k-nearest neighbor algorithm)

這幾乎不用什麼數學了，完全就是直觀的邏輯，那就是每次有新的樣本進來的時候，我們都觀察他在空間中最接近的k個點在哪，然後將這k個點的結果「平均」起來。
我們來簡單做一下分析，假定要使用「Rate」和「QTc」預測「AMI」，假定有個新的樣本\(Rate = 150\)同時\(QTc=600\)，這是原始圖像：

subdat <- dat[!(dat[,"AMI"] %in% NA) & !(dat[,"Rate"] %in% NA) & !(dat[,"QTc"] %in% NA),]
X_mat <- cbind(subdat[,'Rate'], subdat[,'QTc'])
col_idx <- as.integer(subdat[,'AMI'])
col_list <- c('#FF000030', '#00A00030', '#0000FF30')

plot(X_mat[,1], X_mat[,2], pch = 19, cex = 0.5, col = col_list[col_idx], main = 'All data')
points(150, 600, pch = 19, cex = 1)

第三節：k-近鄰演算法(2)

我們先計算這個新的點與所有已知點的距離：

new_vec <- c(150, 600)
new_X <- matrix(new_vec, nrow = nrow(X_mat), ncol = 2, byrow = TRUE)
distance <-  apply((X_mat - new_X)^2, 1, sum)

假定我們要找最接近的20個點來預測新的點，我們可以透過這個方法：

k_param <- 20
limit_val <- sort(distance)[k_param]
used_id <- which(distance <= limit_val)

par(mfrow = c(1, 2))

plot(X_mat[-used_id,1], X_mat[-used_id,2], pch = 19, cex = 0.5, col = '#00000020', main = 'Original scale')
points(X_mat[used_id,1], X_mat[used_id,2], col = col_list[col_idx[used_id]], pch = 19, cex = 1)
points(150, 600, pch = 15, cex = 1)

plot(X_mat[-used_id,1], X_mat[-used_id,2], xlim = c(100, 200), ylim = c(500, 700), pch = 19, cex = 0.5, col = '#00000020', main = 'Interested zoom')
points(X_mat[used_id,1], X_mat[used_id,2], col = col_list[col_idx[used_id]], pch = 19, cex = 1)
points(150, 600, pch = 15, cex = 1)

所以答案就是較有可能是\(AMI = not-AMI\)：

table(subdat[used_id,'AMI'])

## 
## not-AMI  NSTEMI   STEMI 
##      16       1       3

第三節：k-近鄰演算法(3)

你應該有注意到，這樣計算距離對\(QTc\)有較高的加權，所以為了以示公平，通常會先對兩個軸都進行標準化後再計算距離：

norm_X_mat <- scale(X_mat)
new_vec <- c(150, 600)
norm_vec <- (new_vec - attr(norm_X_mat, 'scaled:center')) / attr(norm_X_mat, 'scaled:scale')
new_norm_X <- matrix(norm_vec, nrow = nrow(X_mat), ncol = 2, byrow = TRUE)
distance <-  apply((norm_X_mat - new_norm_X)^2, 1, sum)

在畫圖的時候我們使用原來的scale：

k_param <- 20
limit_val <- sort(distance)[k_param]
used_id <- which(distance <= limit_val)

par(mfrow = c(1, 2))

plot(X_mat[-used_id,1], X_mat[-used_id,2], pch = 19, cex = 0.5, col = '#00000020', main = 'Original scale')
points(X_mat[used_id,1], X_mat[used_id,2], col = col_list[col_idx[used_id]], pch = 19, cex = 1)
points(150, 600, pch = 15, cex = 1)

plot(X_mat[-used_id,1], X_mat[-used_id,2], xlim = c(100, 200), ylim = c(500, 700), pch = 19, cex = 0.5, col = '#00000020', main = 'Interested zoom')
points(X_mat[used_id,1], X_mat[used_id,2], col = col_list[col_idx[used_id]], pch = 19, cex = 1)
points(150, 600, pch = 15, cex = 1)

結果有一點點不一樣：

table(subdat[used_id,'AMI'])

## 
## not-AMI  NSTEMI   STEMI 
##      15       1       4

第三節：k-近鄰演算法(4)

我們來直接使用套件「class」的函數「knn()」進行分析：

– 需要注意的是，由於k-近鄰演算法每次都要跟原始資料比，所以是沒有辦法精簡出摘要數據的，每次我們都需要跟原始資料做比較。

– 我們直接將樣本分成「訓練組」、「驗證組」與「測試組」。

subdat <- dat[!(dat[,"AMI"] %in% NA) & !(dat[,"Rate"] %in% NA) & !(dat[,"QTc"] %in% NA),]
X <- subdat[,c('Rate', 'QTc')]
Y <- subdat[,"AMI"]

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- X[train_idx,]
valid_X <- X[valid_idx,]
test_X <- X[test_idx,]

train_Y <- Y[train_idx]
valid_Y <- Y[valid_idx]
test_Y <- Y[test_idx]

直接看看在驗證組中的表現如何，並比較沒有標準化與有標準化的結果：

library(class)

pred_valid.1 <- knn(train = train_X, test = valid_X, cl = train_Y, k = 20, prob = TRUE)
pred_valid.2 <- knn(train = scale(train_X), test = scale(valid_X), cl = train_Y, k = 20, prob = TRUE)

table(valid_Y, pred_valid.1)

##          pred_valid.1
## valid_Y   not-AMI NSTEMI STEMI
##   not-AMI     212      2     6
##   NSTEMI       24      0     2
##   STEMI        48      1     0

table(valid_Y, pred_valid.2)

##          pred_valid.2
## valid_Y   not-AMI NSTEMI STEMI
##   not-AMI     209      2     9
##   NSTEMI       24      0     2
##   STEMI        48      1     0

練習3：使用k-近鄰演算法分析

承練習2，讓我們再使用「Rate」和「QTc」預測「LVD」，我們一樣先將樣本分成「訓練組」、「驗證組」與「測試組」。

subdat <- dat[!(dat[,"LVD"] %in% NA) & !(dat[,"Rate"] %in% NA) & !(dat[,"QTc"] %in% NA),]
X <- subdat[,c('Rate', 'QTc')]
Y <- subdat[,"LVD"]

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_X <- X[train_idx,]
valid_X <- X[valid_idx,]
test_X <- X[test_idx,]

train_Y <- Y[train_idx]
valid_Y <- Y[valid_idx]
test_Y <- Y[test_idx]

請你比較k-近鄰演算法分析、樸素貝葉斯分類器與邏輯斯回歸哪個更好！

練習3答案

這是k-近鄰演算法(k = 20)的結果(你可以試試其他數值k)：

valid_pred <- knn(train = scale(train_X), test = scale(valid_X), cl = train_Y, k = 20, prob = TRUE)
test_pred <- knn(train = scale(train_X), test = scale(test_X), cl = train_Y, k = 20, prob = TRUE)

library(pROC)

roc_valid <- roc(valid_Y ~ attr(valid_pred, 'prob'), direction = '>')
best_pos <- which.max(roc_valid$sensitivities + roc_valid$specificities)
best_cut <- roc_valid$thresholds[best_pos]

tab_test <- table(attr(test_pred, 'prob') >= best_cut, test_Y)
sens <- tab_test[2,2] / sum(tab_test[,2])
spec <- tab_test[1,1] / sum(tab_test[,1])

roc_test <- roc(test_Y ~ attr(test_pred, 'prob'), direction = '>')
plot(roc_test)

text(0.5, 0.5, paste0('AUC(valid) = ', formatC(roc_valid$auc, digits = 3, format = 'f'),
                      '\nSens = ', formatC(sens, digits = 3, format = 'f'),
                      '\nSpec = ', formatC(spec, digits = 3, format = 'f'),
                      '\nAUC(test) = ', formatC(roc_test$auc, digits = 3, format = 'f')), col = 'red')

同樣的，如果依照驗證組的AUC來選擇的話樸素貝葉斯分類器是最佳選擇，不過我們又多了一種演算法能嘗試了。

課程小結

今天一口氣教了3個機器學習的方法，但你應該有注意到這三個算法都是「分類器」，所以是沒有辦法處理「連續變項」的預測，更別提是「存活分析」的資料了

– 你應該有開始覺得機器學習並沒有這麼困難，在過去3節課我們才分別教會你3個統計方法，而今天一口氣就學會3個機器學習方法，有沒有覺得很有信心！

k-近鄰演算法是我們第一個接觸到帶有非線性預測成分在裡面的算法，其他幾個仍然屬於線性預測的範疇。

– 透過大量的實驗應該有發現，線性預測也不見得不好，事實上大多數時候線性預測帶來的效果並不一定比較差。

從泛用性來看，如果想要用線性模型，統計方法似乎更為理想，可以面對更多不同的資料型態做出預測。下週我們將會把linear regression、logistic regression、softmax regression以及Cox proportional hazards model進行一個總整理，擴展他們的公式並讓他們能對更複雜的資料作出精準的預測。

機器學習及演算法