AS4-0 課堂筆記

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=4, scipen=12)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

---

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)

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call prp with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

rpart.plot(rpart1,cex=0.6)

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call rpart.plot with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

1.4 Make Prediction

pred = predict(rpart1, TS, type = "class")  # predict classes
x = table(actual = TS$Reverse, pred); x

      pred
actual  0  1
     0 41 36
     1 22 71

sum(diag(x))/sum(x)                         # .65882

[1] 0.6588

1.5 ROC & AUC

library(ROCR)

Loading required package: gplots

Attaching package: 'gplots'

The following object is masked from 'package:stats':

    lowess

PredictROC = predict(rpart1, TS)              # predict prob.
head(PredictROC)

        0      1
1  0.3036 0.6964
3  0.3036 0.6964
4  0.4000 0.6000
6  0.4000 0.6000
8  0.4000 0.6000
21 0.3036 0.6964

perf = prediction(PredictROC[,2], TS$Reverse)
perf = performance(perf, "tpr", "fpr")
plot(perf)

library(ROCR)
PredictROC = predict(rpart1, TS)              # predict prob.
head(PredictROC)

         0       1
1  0.30357 0.69643
3  0.30357 0.69643
4  0.40000 0.60000
6  0.40000 0.60000
8  0.40000 0.60000
21 0.30357 0.69643

perf = prediction(PredictROC[,2], TS$Reverse)
perf = performance(perf, "tpr", "fpr")
plot(perf)

1.6 DPP of Decision Tree

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

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call rpart.plot with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

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

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

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

1. 柱狀圖偏左和偏右的形狀，假如此決策樹預測yes的機率較no高，此決策樹的形狀就會偏向右邊，反之若是no的機率較高，便會偏向左邊。
2. (1)minbucket的值設的越小，決策樹最後的分類結果就會越多，DDP柱狀圖也會因此越多。 (2)紅色柱狀和綠色柱狀重疊越多（橘色部分），預測越不準 (3)決策樹上的機率剛好能對應到旁邊數量單位的比率，DDP能將決策樹呈的結果以柱狀圖呈現，方便理解。

1.7 The effect of minbucket

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

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call prp with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

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

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call prp with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

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

• If we set meanbucket too small, it may become overfit to the training set, but if we set it too large, the model will become too simple, either of them cause poor model.

1.8 Random Forest Method - randomForest::randomForest()

library(randomForest)

randomForest 4.6-14
Type rfNews() to see new features/changes/bug fixes.

Attaching package: 'randomForest'

The following object is masked from 'package:dplyr':

    combine

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

[1] 0.6824

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

[1] 0.6882

1.9 Cross Validation & Parameter Tuning - caret package

library(caret)

Loading required package: lattice
Loading required package: ggplot2

Attaching package: 'ggplot2'

The following object is masked from 'package:randomForest':

    margin

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

CART 

396 samples
  6 predictor
  2 classes: '0', '1' 

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 357, 356, 356, 356, 357, 356, ... 
Resampling results across tuning parameters:

  cp    Accuracy  Kappa  
  0.01  0.6136    0.20703
  0.02  0.6137    0.20625
  0.03  0.6238    0.23127
  0.04  0.6314    0.25635
  0.05  0.6441    0.28402
  0.06  0.6441    0.28402
  0.07  0.6441    0.28402
  0.08  0.6441    0.28402
  0.09  0.6441    0.28402
  0.10  0.6441    0.28402
  0.11  0.6441    0.28402
  0.12  0.6441    0.28402
  0.13  0.6441    0.28402
  0.14  0.6441    0.28402
  0.15  0.6441    0.28402
  0.16  0.6441    0.28402
  0.17  0.6441    0.28402
  0.18  0.6441    0.28402
  0.19  0.6210    0.23063
  0.20  0.5830    0.13765
  0.21  0.5830    0.13765
  0.22  0.5530    0.04839
  0.23  0.5479    0.03226
  0.24  0.5454    0.01768
  0.25  0.5454    0.00000
  0.26  0.5454    0.00000
  0.27  0.5454    0.00000
  0.28  0.5454    0.00000
  0.29  0.5454    0.00000
  0.30  0.5454    0.00000
  0.31  0.5454    0.00000
  0.32  0.5454    0.00000
  0.33  0.5454    0.00000
  0.34  0.5454    0.00000
  0.35  0.5454    0.00000
  0.36  0.5454    0.00000
  0.37  0.5454    0.00000
  0.38  0.5454    0.00000
  0.39  0.5454    0.00000
  0.40  0.5454    0.00000
  0.41  0.5454    0.00000
  0.42  0.5454    0.00000
  0.43  0.5454    0.00000
  0.44  0.5454    0.00000
  0.45  0.5454    0.00000
  0.46  0.5454    0.00000
  0.47  0.5454    0.00000
  0.48  0.5454    0.00000
  0.49  0.5454    0.00000
  0.50  0.5454    0.00000

Accuracy was used to select the optimal model using the largest value.
The final value used for the model was cp = 0.18.

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

• cp=0.19, 他有最高的accuracy，設在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

[1] 0.7235

rpart.plot(rpart2)

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call rpart.plot with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

【討論：】

• Classification & Regression Tree rpart classfication呈現類別, regression呈現機率
• Tree Bagging - Random Forest randomForest
  多個決策術得分類器，隨機森林準確度更高，但會犧牲對圖的解釋性，因為個training set切的不一樣所以每個分配的indepedent variable也會不同
• Overfitting : Curse of Complexity overfitting代表在追求預設trainging set data的準確性而製作model的同時，可能將meanbucket,cp等參數設的太過符合training data，導致在預測其他data時反而失真。
• Model Parameter : Cost of Compexity 參數的多寡會影響模型複雜度，而這會影響到建立與使用一個模型的時間成本
• Parameter Tuning by CV (caret package) 透過ＣＶ找到最佳參數，反覆驗證能提高準確性

---

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)

       1        2        3        4        5 
0.671268 0.190170 0.089466 0.043325 0.005771 

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

# Check Percentages Among Subsets
library(dplyr)
table(TR$bucket2009) %>% prop.table

       1        2        3        4        5 
0.671266 0.190169 0.089468 0.043326 0.005771 

table(TS$bucket2009) %>% prop.table

      1       2       3       4       5 
0.67127 0.19017 0.08946 0.04332 0.00577 

table(D$bucket2009) %>% prop.table

       1        2        3        4        5 
0.671268 0.190170 0.089466 0.043325 0.005771 

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

          1      2       3       4        5
[1,] 0.6713 0.1902 0.08947 0.04332 0.005771
[2,] 0.6713 0.1902 0.08947 0.04333 0.005771
[3,] 0.6713 0.1902 0.08946 0.04332 0.005770

# Data Exploration 
mean(TR$age)                 # 72.64

[1] 72.64

mean(TR$diabetes)            # .3809 

[1] 0.3809

2.2 Baseline Model

# Baseline Accuracy
acc.base = table(TS$bucket2009, TS$bucket2008) %>% 
  {sum(diag(.)) / sum(.)}    # .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

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    2    0    1    2    3
[3,]    4    2    0    1    2
[4,]    6    4    2    0    1
[5,]    8    6    4    2    0

# 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

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    2    0    1    2    3
[3,]    4    2    0    1    2
[4,]    6    4    2    0    1
[5,]    8    6    4    2    0

# if we use the most popular class as baseline
acc.dumb = (table(TS$bucket2009)/nrow(TS))[1]                 # acc: .6713 
pen.dumb = sum(table(TS$bucket2009) * c(0,2,4,6,8))/nrow(TS)  # pen: 1.044

CART Model

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

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call prp with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

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

CART Model with Customized Loss Matrix

rpart2 = rpart(bucket2009 ~ ., data=TR[c(1:14, 16)], 
               method="class", cp=0.00005, 
               parms=list(loss=Penalty))
pred = predict(rpart2, TS, type = "class")
acc.rpart2 = table(TS$bucket2009, pred) %>% 
  {sum(diag(.))/sum(.)}                  # acc.rpart2: .64727   
cm = as.matrix(table(TS$bucket2009, pred))
pen.rpart2 = sum(cm*Penalty)/nrow(TS)    # pen.rpart2: .64182 

rpart2 = rpart(bucket2009 ~ ., data=TR[c(1:14, 16)], 
               method="class", cp=0.00005, 
               parms=list(loss=Penalty)     # customized loss function 將penalty納入模型考慮【
               )

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

pred = predict(rpart2, TS, type = "class")
acc.rpart2 = table(TS$bucket2009, pred) %>% {sum(diag(.))/sum(.)}
acc.rpart2    # acc.rpart2 = .64727   

[1] 0.64727

cm = as.matrix(table(TS$bucket2009, pred))
pen.rpart2 = sum(cm*Penalty)/nrow(TS)    
pen.rpart2    # pen.rpart2: .64182 

[1] 0.64182

【討論】

• 5 classes: 5X5 confussion/penalty(payoff) matrices penalty成本矩陣
• dumb vs. smart baseline dumb 以最多的類別作為預測基準 （最多是1,全預測為1) smart 以之前的資料作為預測基準 (去年你是2，今年也預測你2)
• model1: accuracy up, penalty up accuracy最高時，penalty不一定對低
• model2: accuracy down, penalty down 加入懲罰矩陣當參數後，會選擇penalty最低的結果，但accuracy不一定是最高的
• accuracy != profitability
• “customized” model optimize criteria in classification methods (not always available) penalty矩陣是我們自己訂的，不一定在每個model都適用
• This CV/PT process is over simplified cp最好經過反覆驗證之後再決定
• 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))

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

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.385   0.449   0.538   0.555   0.624   0.871 

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)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0    17.0    21.2    22.5    25.0    50.0 

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

Linear Regression using LAT and LON

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

Call:
lm(formula = MEDV ~ LAT + LON, data = D)

Residuals:
   Min     1Q Median     3Q    Max 
-16.46  -5.59  -1.30   3.69  28.13 

Coefficients:
            Estimate Std. Error t value          Pr(>|t|)    
(Intercept) -3178.47     484.94   -6.55 0.000000000138619 ***
LAT             8.05       6.33    1.27               0.2    
LON           -40.27       5.18   -7.77 0.000000000000045 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.69 on 503 degrees of freedom
Multiple R-squared:  0.107, Adjusted R-squared:  0.104 
F-statistic: 30.2 on 2 and 503 DF,  p-value: 0.000000000000416

Visualize regression output

# 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)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.385   0.449   0.538   0.555   0.624   0.871 

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)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0    17.0    21.2    22.5    25.0    50.0 

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

Load CART packages

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")

Simplify tree by increasing minbucket

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

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")

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)

Call:
lm(formula = MEDV ~ LAT + LON + CRIM + ZN + INDUS + CHAS + NOX + 
    RM + AGE + DIS + RAD + TAX + PTRATIO, data = train)

Residuals:
   Min     1Q Median     3Q    Max 
-14.51  -2.71  -0.68   1.79  36.88 

Coefficients:
              Estimate Std. Error t value    Pr(>|t|)    
(Intercept) -252.29051  436.71008   -0.58       0.564    
LAT            1.54396    5.19196    0.30       0.766    
LON           -2.98711    4.78583   -0.62       0.533    
CRIM          -0.18080    0.04391   -4.12 0.000047725 ***
ZN             0.03250    0.01877    1.73       0.084 .  
INDUS         -0.04297    0.08473   -0.51       0.612    
CHAS           2.90418    1.22005    2.38       0.018 *  
NOX          -21.61279    5.41371   -3.99 0.000079778 ***
RM             6.28440    0.48270   13.02     < 2e-16 ***
AGE           -0.04430    0.01785   -2.48       0.014 *  
DIS           -1.57736    0.28418   -5.55 0.000000056 ***
RAD            0.24509    0.09728    2.52       0.012 *  
TAX           -0.01112    0.00545   -2.04       0.042 *  
PTRATIO       -0.98349    0.19389   -5.07 0.000000638 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.6 on 350 degrees of freedom
Multiple R-squared:  0.665, Adjusted R-squared:  0.653 
F-statistic: 53.4 on 13 and 350 DF,  p-value: <2e-16

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

[1] 3037

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

[1] 4329

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

CART 

364 samples
 13 predictor

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 327, 326, 328, 329, 328, 327, ... 
Resampling results across tuning parameters:

  cp     RMSE   Rsquared  MAE  
  0.000  4.729  0.7598    3.068
  0.001  4.745  0.7584    3.105
  0.002  4.756  0.7569    3.091
  0.003  4.807  0.7491    3.151
  0.004  4.876  0.7418    3.230
  0.005  4.917  0.7368    3.281
  0.006  4.954  0.7315    3.322
  0.007  4.899  0.7360    3.290
  0.008  4.933  0.7297    3.323
  0.009  4.906  0.7320    3.308
  0.010  4.906  0.7320    3.308

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was cp = 0.

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

Bad 'data' field in model 'call' field.
         To make this warning go away:
             Call prp with roundint=FALSE,
             or rebuild the rpart model with model=TRUE.

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

[1] 3660

Ensembling

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

[1] 2599

#       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 linear regression 可以看每個參數對結果的邊際效用 tree 是分析非線性關係
• the concept of ensembling 把兩個簡單的模型再做一個新的模型，可以綜效兩個模型觀察到的特性

---