22-08-25 chap16_DecisionTree

chap16_DecisionTree

 

install.packages('rpart') # rpart() : 분류모델 생성 
install.packages("rpart.plot") # prp(), rpart.plot() : rpart 시각화
install.packages('rattle') # fancyRpartPlot() : node 번호 시각화 

library(rpart) 
library(rpart.plot) 
library(rattle) 

####################################
# 암 진단 분류 분석 : 이항분류
####################################
# "wdbc_data.csv" : 유방암 진단결과 데이터 셋 분류

# 1. 데이터셋 가져오기 
wdbc <- read.csv('C:/itwill/3_Rwork/data/wdbc_data.csv')
str(wdbc) # 569 obs. of  32 variables:

# 2. 데이터 탐색 및 전처리 
wdbc_df <- wdbc[-1] # id 칼럼 제외
head(wdbc_df)
wdbc_df$diagnosis # 진단결과 : B -> '양성', M -> '악성'
table(wdbc_df$diagnosis)

# 목표변수(y변수)를 factor형 변환 : 0, 1 dummy 변수  
wdbc_df$diagnosis <- as.factor(wdbc$diagnosis) # 0, 1
wdbc_df$diagnosis[1:10] 


# 3. 훈련셋과 검정셋 생성 : 7 : 3 비율 
set.seed(415)
idx = sample(nrow(wdbc_df), 0.7*nrow(wdbc_df))
wdbc_train = wdbc_df[idx, ] # 훈련 데이터 
wdbc_test = wdbc_df[-idx, ] # 검정 데이터 
dim(wdbc_train) # 398  31
dim(wdbc_test) # 171  31

# 4. rpart 분류모델 생성 
model_wdbc <- rpart(diagnosis ~ ., data = wdbc_train)
model_wdbc # - 중요한 x변수 확인 
# 1) root 398 136 B (0.65829146 0.34170854) 
#    398:전체관측치, 136:다른범주관측치, B:높은비율범주 (분류비율)
#   2) points_mean< 0.04923 255  10 B (0.96078431 0.03921569) *
#   3) points_mean>=0.04923 143  17 M (0.11888112 0.88111888)
#  번호) <분류조건> 분류된관측치 다른범주관측치 높은비율범주 (분류비율) 단노드여부     

# tree 시각화 
rpart.plot(model_wdbc)
# 노드-상 : 높은비율 범주(label) 
# 노드-중 : 더미변수=1 분류비율 
# 노드-하 : 분류된 전체비율 

398 - 136 # 262
table(wdbc_train$diagnosis)
# B   M 
#262 136

left_node <- subset(wdbc_train, points_mean< 0.04923)
dim(left_node) # 255  31
table(left_node$diagnosis)
#  B   M 
#245  10 

right_node <- subset(wdbc_train, points_mean >= 0.04923)
dim(right_node) # 143  31
table(right_node$diagnosis)
# B   M 
#17 126

# [정리] points_mean 평균 크기가 클 수록 악성 비율이 높다.

# 노드번호 제공 
fancyRpartPlot(model_wdbc)


# 5. 모델 평가 
y_pred <- predict(object=model_wdbc, newdata=wdbc_test, type='prob') # y 확률 예측 
y_pred2 <- predict(object=model_wdbc, newdata=wdbc_test, type='class') # y class 예측 
# type = “prob” : 확률 예측(기본)
# type = “class” : 범주 예측 
range(y_pred) # 0.01626016 ~ 0.98373984
dim(y_pred) # 171   2
str(y_pred) # num [1:171, 1:2]

# [시험] 각 관측치별 악성의 확률 -> DF
y_pred_df <- data.frame(m_prop = y_pred[,2])
str(y_pred_df)
head(y_pred_df)


# 확률예측 : 악성(M) 기준 
y_pred <- ifelse(y_pred[,2] > 0.5, 1, 0) 
y_pred

y_true <- wdbc_test$diagnosis # y 관측치 

# 혼동행렬 
tab <- table(y_true, y_pred)
tab
#      y_pred
#y_true  B(0)  M(1)
#     B(0) 92(TN)  3(FP)
#     M(1) 13(FN) 63(TP)


acc <- sum(diag(tab)) / sum(tab)
cat('분류정확도 =',acc) # 분류정확도 = 0.9064327

table(wdbc_test$diagnosis)
# B  M 
#95 76

# 정확률 : 예측치 기준 
p <- tab[2,2] / sum(tab[,2]) # TP / FP+TP

# 재현률 : 관측치 기준 
r <- tab[2,2] / sum(tab[2,]) # TP / FN+TP

f1 <- 2*((p*r) / (p+r)) # 조화평균 
f1 # 0.8873239



################################
# iris 데이터셋 : 다항분류 
################################
# 단계1. 실습데이터 생성 
data(iris)
set.seed(415)
idx = sample(nrow(iris), 0.7*nrow(iris))
train = iris[idx, ]
test = iris[-idx, ]
dim(train) # 105 5
dim(test) # 45  5

table(train$Species)

# 단계2. 분류모델 생성 
# rpart(y변수 ~ x변수, data)
?rpart # cp = 0.01(가지치기 기본값)

# 가지치기(cp)
# - 과적합 해결방안 
# - cp 범위 : 0 ~ 1, default=0.01
# 0에 가까울 수록 tree depth 깊어짐(과적합(분산) 증가, 오류율 감소) 
# 1에 가까울 수록 tree depth 낮아짐(과적합(분산) 감소, 오류율 증가)  


model = rpart(Species~., data=train) # iris의 꽃의 종류(Species) 분류 
model

names(model)

# 현재 cp 정보 제공 
model$cptable
#   CP nsplit        rel error     xerror       xstd
#1 0.5147059      0 1.00000000 1.16176471 0.06504241
#2 0.4558824      1 0.48529412 0.60294118 0.07351540 : tree depth=1
#3 0.0100000      2 0.02941176 0.04411765 0.02510482 : 현재 cp(tree depth=2) 

# 분류모델 시각화 - rpart.plot 패키지 제공 
prp(model) # 간단한 시각화   
rpart.plot(model) # rpart 모델 tree 출력
fancyRpartPlot(model) # node 번호 출력(rattle 패키지 제공)


# 단계3. 분류모델 평가  
pred <- predict(model, test, type="prob") # 비율 예측(기본) 
pred <- predict(model, test, type="class") # 분류 예측 

# 1) 분류모델로 분류된 y변수 보기 
table(pred)

# 2) 분류모델 성능 평가 
tab <- table(pred, test$Species)

sum(diag(tab)) / sum(tab) # 0.9111111


###################################
### 회귀트리 
###################################
# 독립변수(설명변수) : 범주형 또는 연속형 변수  
# 종속변수(반응변수) : 연속형 변수
# 회귀트리 평가방법 : MSE, R2 score 


library(ggplot2) # dataset 사용 
data(mpg)
str(mpg) 

# 단계1 : 학습데이터와 검정데이터 샘플링
idx <- sample(nrow(mpg), nrow(mpg)*0.7)
train <- mpg[idx, ]
test <- mpg[-idx, ]
dim(train) # 163 11


# 단계2 : 학습데이터 이용 분류모델 생성 
model <- rpart(cty ~ displ + cyl + year, data = train)
model # 중요변수 : displ
range(train$cty) # 9 ~ 35 : 연속형 


# 단계3 : 검정데이터 이용 예측치 생성 및 평가 
y_pred <- predict(model, test) # type='class' : 범주형에서 사용 
y_pred
y_true <- test$cty

# 결정계수(R2 score) : 1 수렴정도(Y변수 표준화 X) 
R2_score <- cor(y_true, y_pred)^2 # R^2
R2_score# 0.8382528

# 평균제곱오차(MSE) : 0 수렴정도(Y변수 표준화 0) 
MSE <- mean((y_true-y_pred)^2)
MSE # 3.247795


# 단계4 : 분류분석 결과 시각화
rpart.plot(model) 

range(train$cty)# 9 ~ 35
# [해설] 엔진 크기가 클 수록 cty 연비는 낮고, 작을 수록 cty 연비는 좋다. 

displ_2.6_g <- subset(train, displ >= 2.6)
displ_2.6_l <- subset(train, displ < 2.6)

mean(displ_2.6_g$cty) # 14.64815
mean(displ_2.6_l$cty) # 21.09091


################################
###  k겹 교차검증 
################################
install.packages('cvTools')
library(cvTools)


# 단계1 : # K겹 교차검정을 위한 샘플링 
cross <- cvFolds(nrow(iris), K=5) # 5등분 균등 분할 
cross # Fold(dataset 구분)   Index(행번호)
str(cross) # List of 5

# 5개 데이터셋 구성 : which 이용 색인으로 사용   
dataset1 <- cross$subsets[cross$which==1, 1]
dataset2 <- cross$subsets[cross$which==2, 1]
dataset3 <- cross$subsets[cross$which==3, 1]
dataset4 <- cross$subsets[cross$which==4, 1]
dataset5 <- cross$subsets[cross$which==5, 1]

dataset1
length(dataset1) # 30
length(dataset5) # 30

# 단계2 : K겹 교차검정 
K = 1:5 # 5겹(1 2 3 4 5) 
ACC <- c() # 분류정확도 누적  

for(i in K){ # 5회 반복 
  # 1. train/test split 
  idx <- cross$subsets[cross$which==i, 1] # [i=1, 1]행번호 추출  
  test <- iris[idx, ] # 평가셋(30) 
  train <- iris[-idx, ] # 훈련셋(120) 
  # 2. 모델 생성 
  model <- rpart(Species ~ ., data = train)
  # 3. 예측치 
  y_pred <- predict(model, test, type='class')  
  # 4. 혼동행렬 
  t <- table(test$Species, y_pred)
  # 5. 분류정확도 
  ACC <- c(ACC,  sum(diag(t)) / sum(t) )
}


ACC # 0.8666667 0.9666667 1.0000000 0.9333333 0.9333333

# 단계3 : 교차검정 평가 
cat('분류정확도 산술평균 =', mean(ACC))


#####################################
### Entropy vs GINI : 불확실성 척도 
#####################################

# 사례1) p1 : 앞면, p2 : 뒷면 - 불확실성이 높은 경우 
p1 = 0.5; p2 = 0.5 # 확률변수 

# Entropy = -sum(p * log(p))
e1 <- -(p1 * log2(p1)) + -(p2 * log2(p2))
e1 <- -sum(c(p1 * log2(p1), p2 * log2(p2)))
e1 # 1 -> 엔트로피 가장 큰 경우 

# GINI = p1 * (1-p1) + p2 * (1-p2)
gini1 <- p1 * (1 - p1) + p2 * (1 - p2)
gini1 # 0.5 -> gini 계수 가장 큰 경우 


# 사례2) p1 : 앞면, p2 : 뒷면 - 불확실성이 낮은 경우 
p1 = 0.9; p2 = 0.1

e2 <- -(p1 * log2(p1)) + -(p2 * log2(p2))
e2 # 0.4689956

gini2 <- p1 * (1 - p1) + p2 * (1 - p2)
gini2 # 0.18


###################################
### 분류분석 적용 예(ppt.18)
###################################
# - 먹구름(x1)과 돌풍(x2)에 대한 비 유무 예측  

#먹구름(x1) 돌풍(x2)  비(y) 
#   1        1         yes
#   1        1         yes
#   1        0         no
#   0        1         no
#   0        1         no


# 단계1 : 전체 엔트로피(Entropy Root)
hit_prop = 2/5 # yes 확률 
no_hit_prop = 3/5 # no 확률 
# Entropy 수식 
entropyRoot = - (hit_prop*log2(hit_prop) + no_hit_prop*log2(no_hit_prop))
entropyRoot # 0.9709506


# 단계2 : 각 x변수 엔트로피 
# 1) x1 변수(YYN/NN) 엔트로피 
YYN = -(2/3*log2(2/3)) + -(1/3*log2(1/3)) # x1=1 entropy  
NN = -(0) + -(2/2*log2(2/2)) # x1=0  entropy 
x1_Entropy = (3/5)*YYN + (2/5)*NN # 가중치 * entropy 적용 
x1_Entropy # 0.5509775


# 2) x2 변수(YYNN/N) 엔트로피 
YYNN = -(2/4*log2(2/4)) + -(2/4*log2(2/4)) # x2=1 entropy   
N = -(0) + -(1/1*log2(1/1))# x2=0 entropy
x2_Entropy = (4/5)*YYNN + (1/5)*N 
x2_Entropy # 0.8


# 단게3 : 정보이득(Information gain) 
x1_infoGain <- entropyRoot - x1_Entropy  # 0.4199731
x2_infoGain <- entropyRoot - x2_Entropy  # 0.1709506

# 해설 : 정보이득이 높은 x1 변수를 중요변수로 선정