AS4-0 Gruop4 課堂筆記

packages = c(
  "dplyr","ggplot2","d3heatmap","googleVis","devtools","plotly", "xgboost",
  "magrittr","caTools","ROCR","corrplot", "rpart", "rpart.plot",
  "doParallel", "caret", "glmnet", "Matrix", "e1071", "randomForest"
  )
existing = as.character(installed.packages()[,1])
for(pkg in packages[!(packages %in% existing)]) install.packages(pkg)

rm(list=ls(all=T))
options(digits=5, scipen=12)
library(dplyr)
library(caTools)
library(rpart)
library(rpart.plot)
library(randomForest)

---

1. Suprime Court Decision

1.1 Prepare Data

D = read.csv("data/stevens.csv")

library(caTools)
set.seed(3000)
spl = sample.split(D$Reverse, SplitRatio = 0.7)
TR = subset(D, spl)
TS = subset(D, !spl)

1.2 CART (Classification & Regression Tree) - rpart::rpart()

library(rpart)
rpart1 = rpart(Reverse ~ ., TR[,3:9],          # simplify the formula
               method="class", minbucket=25)

1.3 Plot Decision Tree

library(rpart.plot)
prp(rpart1)
rpart.plot(rpart1,cex=0.6)

1.4 Make Prediction

pred = predict(rpart1, TS, type = "class")  # predict classes
x = table(actual = TS$Reverse, pred); x
sum(diag(x))/sum(x)                         # .65882

1.5 ROC & AUC

library(ROCR)
PredictROC = predict(rpart1, TS)              # predict prob.
head(PredictROC)
perf = prediction(PredictROC[,2], TS$Reverse)
perf = performance(perf, "tpr", "fpr")
plot(perf)

pred = predict(rpart1, TS)[,2]     # prob. of Reverse = 1         
colAUC(pred, TS$Reverse, T)        # AUC = 0.6927

1.6 DPP of Decision Tree

rpart.plot(rpart1,cex=0.75,varlen=3,faclen=2,box.palette="GnRd")

source("DPP.R")
auc = DPP(predict(rpart1)[,2], TR$Reverse, 1)     # AUC = 0.7884

【QUIZ-1】比較以上prp()、rpart.plot()和DPP()這幾張圖形

• 一棵決策樹的形狀， 和由它所產生的預測機率分布之間，有什麼關係呢 ? 以上DPP 圖中的分布點數、位置和高度，分別會對應到決策樹的什麼特徵呢?

• 決策樹的分類是利用現有的變數，將資料做最好的區隔變數來適配，於是藉由設定參數讓決策樹演算出最好的模型，每個分類下的資料點會計算出其機率，再對應其類別。將所有資料點依演算法做出的模型做切割，每個資料點會到不同的分類(不同分類有不同的機率與類別)，而這分類的過程也就是決策樹的分支。/DPP的橫軸為機率；縱軸為數量，將決策樹分類後的每個bucket有其機率則其會分布在DPP所對應的機率上，而因此基本上bucket數與DPP上的長條數會一樣(本題因兩個機率相同則合併計算)，另外決策樹的bucket會顯示%數為所含總資料數的比例，反應在DPP的高度。

1.7 The effect of minbucket

t5 = rpart(Reverse ~ ., TR[,3:9], method="class", minbucket=5)
prp(t5)

t100 = rpart(Reverse ~ ., TR[,3:9], method="class", minbucket=100)
prp(t100)

【QUIZ-2】Based on the 2 plots above, what is the effect of minbucket?

• Minbucket 代表在決策樹做分類時，每個bucket所需包含的最少數量，這個設定可以避免過度適配原有資料而導致決策樹過於複雜流失更多的ACC。

1.8 Random Forest Method - randomForest::randomForest()

library(randomForest)
TR$Reverse = as.factor(TR$Reverse)  # Convert outcome to factor
TS$Reverse = as.factor(TS$Reverse)

# Build model
rf1 = randomForest(Reverse ~ ., TR[,3:9], ntree=200, nodesize=25)

# Make predictions
pred = predict(rf1, TS)
x = table(TS$Reverse, pred)
sum(diag(x)) / sum(x)          # 0.6824

pred = predict(rf1, TS, type='prob')[,2]
table(TS$Reverse, pred > 0.5) %>% {sum(diag(.))/sum(.)}

1.9 Cross Validation & Parameter Tuning - caret package

library(caret)
cv1 = train(
  Reverse ~ ., TR[,3:9], method = "rpart", 
  trControl = trainControl(method="cv", number=10), # 10 fold CV
  tuneGrid=expand.grid(cp = seq(0.01,0.5,0.01))     # parameter combination
  )
cv1; plot(cv1)

【QUIZ-3】Which cp which we should we choose? Why?

• 0.19， cp = cost of complexity 多設一個節點，損失下降如果比設定的cp低，就不再往下分支了；cp越高，代表model越不複雜。 從交叉驗證選擇最佳的cp的過程中發現本題cp=0.06~0.19的ACC皆屬一致，在ACC一致的情況下會傾向選擇複雜度較低的cp值，故選擇0.19。

1.10 The Final Model

rpart2 = rpart(Reverse ~ ., TR[,3:9], method="class", cp=0.19)
pred = predict(rpart2, TS, type='prob')[,2]
table(TS$Reverse, pred > 0.5) %>% {sum(diag(.)) / sum(.)}  # 0.7235

rpart.plot(rpart2)

【學習重點】

• Classification & Regression Tree rpart
  • method=“classs” : 用來分類
  • predict時，type可為“class”或“prob”
• A Bag of Tree - Random Forest randomForest
  • 如果要分類的話要把y變數轉換為factor型態
  • 可以調整mtry, ntree, nodesize, cp
• Overfitting : Curse of Complexity
  • 樹太過複雜導致
  • 透過參數調校解決問題
• Method Parameter : Cost of Compexity
  • 每多設一個節點會增加複雜度的成本，所以透過設定cp參數來限制樹的複雜度
  • 可以利用CV選擇最佳cp
• Parameter Tuning by CV (caret package)
  • 利用Cross-Validation來選擇最佳的cp
  • 把原始training data分成k份，每次輪流選擇一份當作testing data其他k-1份當作training data，建立多個模型來找最好的參數

---

2 D2HawkEye

2.1 Prepare & Examine Data

rm(list=ls(all=TRUE))
D = read.csv("data/ClaimsData.csv")

# Percentage of patients in each cost bucket
table(D$bucket2009)/nrow(D)

# Split the data
library(caTools)
set.seed(88)
spl = sample.split(D$bucket2009, SplitRatio = 0.6)
TR = subset(D, spl)
TS = subset(D, !spl)

# sapply(): Apply Function & Aggregate Output
sapply(list(D, TR, TS), function(x) 
  table(x$bucket2009) %>% prop.table) %>% t

# Data Exploration 
mean(TR$age)                 # 72.64
mean(TR$diabetes)            # .3809 

2.2 Smart Baseline

# Baseline Accuracy
acc.base = table(TS$bucket2009, TS$bucket2008) %>% {sum(diag(.)) / sum(.)}
acc.base   # acc.base = .6838

# Penalty Matrix
Penalty = matrix(c(
  0,1,2,3,4,
  2,0,1,2,3,
  4,2,0,1,2,
  6,4,2,0,1,
  8,6,4,2,0), byrow=TRUE, nrow=5)
Penalty

# Penalty Error of Baseline Method
cm = as.matrix(table(TS$bucket2009, TS$bucket2008))
pen.base = sum(cm*Penalty) / nrow(TS)
pen.base  # pen.base: .7386

# Dumb baseline - Using the most frequent class
acc.dumb = (table(TS$bucket2009)/nrow(TS))[1]                 
pen.dumb = sum(table(TS$bucket2009) * c(0,2,4,6,8))/nrow(TS)  
c(acc.dumb, pen.dumb)   # 0.6713 1.0443

2.3 CART Model

library(rpart)
library(rpart.plot)
rpart1 = rpart(bucket2009 ~ ., data=TR[c(1:14, 16)], 
               method="class", cp=0.00005)
prp(rpart1)

# Make predictions
pred = predict(rpart1, newdata = TS, type = "class")
acc.rpart1 = table(TS$bucket2009, pred) %>% { sum(diag(.)) /nrow(TS)}
acc.rpart1    # acc.rpart1 = .71267 <- .68381
cm = as.matrix(table(TS$bucket2009, pred))
pen.rpart1 = sum(cm*Penalty)/nrow(TS)
pen.rpart1    # pen.rpart1: .75789 <- .7386

2.4 CART Model with Customized Loss Matrix

rpart2 = rpart(bucket2009 ~ ., data=TR[c(1:14, 16)], 
               method="class", cp=0.00005, 
               parms=list(loss=Penalty) )    # customized loss function # # 原本決策樹選擇ACC最高的結果，但將penalty矩陣加入參數，將會選擇penalty最低的模型。
               

pred = predict(rpart2, TS, type = "class")
acc.rpart2 = table(TS$bucket2009, pred) %>% {sum(diag(.))/sum(.)}
acc.rpart2    # acc.rpart2 = .64727   
cm = as.matrix(table(TS$bucket2009, pred))
pen.rpart2 = sum(cm*Penalty)/nrow(TS)    
pen.rpart2    # pen.rpart2: .64182 

# Summary
data.frame(
  accuracy = c(acc.dumb, acc.base, acc.rpart1, acc.rpart2),
  penalty = c(pen.dumb, pen.base, pen.rpart1, pen.rpart2),
  row.names = c("dumb","base","CART","CART w/LossMx"))

【學習重點】

• 5 classes: 5X5 confussion/penalty(payoff) matrices
• dumb vs. smart baseline
• model1: accuracy up, penalty up
• model2: accuracy down, penalty down
• accuracy != profitability
• “customized” model optimize criteria in classification methods (not always available)
• This CV/PT process is over simplified
• we will cover a better process next time

---

3 Boston House Price

rm(list=ls(all=TRUE))
D = read.csv("data/boston.csv")

3.1 Examine Data

# Plot observations
plot(D$LON, D$LAT)

# Tracts alongside the Charles River
with(D, points(LON[CHAS==1], LAT[CHAS==1], col="blue", pch=19))
points(D$LON[D$CHAS==1], D$LAT[D$CHAS==1], col="blue", pch=19)

# Plot MIT
with(D, points(LON[TRACT==3531], LAT[TRACT==3531], col="red", pch=19) )

# Plot polution
summary(D$NOX)
points(D$LON[D$NOX>=0.55], D$LAT[D$NOX>=0.55], col="green", pch=20)

# Plot prices
plot(D$LON, D$LAT)
summary(D$MEDV)
points(D$LON[D$MEDV>=21.2], D$LAT[D$MEDV>=21.2], col="red", pch=20)

3.2 Linear Regression using LAT and LON

latlonlm = lm(MEDV ~ LAT + LON, data=D)
summary(latlonlm)

3.3 Visualize regression output

plot(D$LON, D$LAT)
points(D$LON[D$MEDV>=21.2], D$LAT[D$MEDV>=21.2], col="red", pch=20)

# latlonlm$fitted.values
points(D$LON[latlonlm$fitted.values >= 21.2], 
       D$LAT[latlonlm$fitted.values >= 21.2], 
       col="blue", pch="x")

library(rpart)
library(rpart.plot)
# CART model
latlontree = rpart(MEDV ~ LAT + LON, data=D)
prp(latlontree)

# Visualize output
plot(D$LON, D$LAT)
points(D$LON[D$MEDV>=21.2], D$LAT[D$MEDV>=21.2], 
       col="red", pch=20)

fittedvalues = predict(latlontree)
points(D$LON[fittedvalues>21.2], 
       D$LAT[fittedvalues>=21.2], col="blue", pch="x")

3.4 Simplify tree by increasing minbucket

latlontree = rpart(MEDV ~ LAT + LON, data=D, minbucket=50)
rpart.plot(latlontree)

# Visualize Output
plot(D$LON,D$LAT)
abline(v=-71.07)
abline(h=42.21)
abline(h=42.17)
points(D$LON[D$MEDV>=21.2], 
       D$LAT[D$MEDV>=21.2], col="red", pch=20)

3.5 Use all the variables

# Split the data
library(caTools)
set.seed(123)
split = sample.split(D$MEDV, SplitRatio = 0.7)
train = subset(D, split==TRUE)
test = subset(D, split==FALSE)

# Create linear regression
linreg = lm(MEDV ~ LAT + LON + CRIM + ZN + INDUS + CHAS + NOX + RM + 
              AGE + DIS + RAD + TAX + PTRATIO, data=train)
summary(linreg)

# Make predictions
linreg.pred = predict(linreg, newdata=test)
linreg.sse = sum((linreg.pred - test$MEDV)^2)
linreg.sse

3.6 Create a CART model

tree = rpart(MEDV ~ LAT + LON + CRIM + ZN + INDUS + CHAS + NOX + RM + 
               AGE + DIS + RAD + TAX + PTRATIO, data=train)
prp(tree)

# Make predictions
tree.pred = predict(tree, newdata=test)
tree.sse = sum((tree.pred - test$MEDV)^2)
tree.sse

3.7 Load libraries for cross-validation

library(caret)

# Number of folds
tr.control = trainControl(method = "cv", number = 10)

# cp values
cp.grid = expand.grid( .cp = (0:10)*0.001)

# Cross-validation
tr = train(MEDV ~ LAT + LON + CRIM + ZN + INDUS + CHAS + NOX + RM + 
             AGE + DIS + RAD + TAX + PTRATIO, 
           data = train, method = "rpart", 
           trControl = tr.control, tuneGrid = cp.grid)
tr

# Extract tree
best.tree = tr$finalModel
prp(best.tree)

# Make predictions
best.tree.pred = predict(best.tree, newdata=test)
best.tree.sse = sum((best.tree.pred - test$MEDV)^2)
best.tree.sse

3.8 Ensembling

en.pred = (best.tree.pred + linreg.pred)/2
sum((en.pred - test$MEDV)^2)

#       METHOD   SSE
#           lm  3037
#         rpart 4329
#      rpart.cv 3660
#   lm+rpart.cv 2599  

【學習重點】

• drawing map with spatial data
• tree in spatial data
• tree vs. linear regression
• the concept of ensembling

---