데이터마이너를 꿈꾸며

회귀앙상블 - 랜덤포레스트(링크) 

분류앙상블 – 배깅, 부스팅, 랜덤 포레스트

 1) 데이터 입력

> setwd('c:/Rwork')

> german=read.table('germandata.txt',header=T)

> german$numcredits = factor(german$numcredits)

> german$residence = factor(german$residence)

> german$residpeople = factor(german$residpeople)

2) 랜덤 포레스트 방법의 실행

> library(randomForest)

> rf.german<-randomForest(y~.,data=german,ntree=100,mytr=5,importance=T,na.action = na.omit)

> summary(rf.german)

Length Class  Mode    

call               7   -none- call    

type               1   -none- character

predicted       1000   factor numeric 

err.rate         300   -none- numeric 

confusion          6   -none- numeric 

votes           2000   matrix numeric 

oob.times       1000   -none- numeric 

classes            2   -none- character

importance        80   -none- numeric 

importanceSD      60   -none- numeric 

localImportance    0   -none- NULL    

proximity          0   -none- NULL    

ntree              1   -none- numeric 

mtry               1   -none- numeric 

forest            14   -none- list    

y               1000   factor numeric 

test               0   -none- NULL    

inbag              0   -none- NULL    

terms              3   terms  call

함수 설명:

randomForest(y~.,data=german,ntree=100,mytr=5,importance=T,na.action = na.omit ntree는 랜덤포레스트 방법에서 생성하게 될 분류기의 개수 B이다. 100개의 분류기는 분류나무를 이용하여 생성되도록 하였다. mtry는 중간노드마다 랜덤하게 선택되는 변수의 개수를 의미한다. 5개의 입력변수를 중간노드마다 모두 선택하고 이 중에서 분할점이 탐색되도록 하였다. importance는 변수의 중요도를 계산하는 옵션이고, na.action은 결측치를 처리하는 방법에 대한 사항이다. 변수의 중요도를 계산하게 하고, 결측치는 필요한 경우에만 삭제되도락 하였다.

2) 랜덤 포레스트 방법의 실행 – 변수 중요도 보기

> names(rf.german)

[1] "call"            "type"            "predicted"       "err.rate"        "confusion"     

[6] "votes"           "oob.times"       "classes"         "importance"      "importanceSD"  

[11] "localImportance" "proximity"       "ntree"           "mtry"            "forest"        

[16] "y"               "test"            "inbag"           "terms"         

> head(rf.german$predicted,10) #데이터의 예측집단을 출력

1    2    3    4    5    6    7    8    9   10

good  bad good  bad  bad good good good good  bad

Levels: bad good

> importance(rf.german)

bad        good MeanDecreaseAccuracy MeanDecreaseGini

check       13.14268724  9.63208341         15.1106884       44.838547

duration     3.33563217  8.10482760          8.6281030        41.273512

history      3.87863720  4.50449203            5.8685013        25.731461

purpose      2.67025503  3.09871408            4.1078354        35.651165

credit       2.44312087  4.16670498            4.8800356        52.641745

savings      7.48182326  2.84190645            6.1879918        21.705909

employment   1.91049595  2.70568977            3.2416484        23.910509

installment  0.02100147  3.49542966            3.0522029        15.714499

personal     2.10802207  1.83819513            3.0602966        15.365475

debtors     -0.17277289  4.39384503            3.8963219         6.718969

residence    0.74571096  0.90284661            1.1155901        17.932937

property     2.16016716  3.85454658            4.7010080        18.803008

age          2.69637542  4.35170316            5.4753585        41.451103

others       0.31569112  3.60499256            3.3679530        11.399935

housing      2.70314243  2.06074416            3.0737900         9.596539

numcredits  -0.24996827  0.95259106            0.6502566         7.944861

job          1.53070574  1.18486660            2.1420488        12.054190

residpeople  0.88657814 -0.43166449            0.1370845         4.234831

telephone    1.46824003 -0.24200291            0.6110284         6.427048

foreign      1.26297478 -0.05431566            0.5733125         1.519832

> importance(rf.german,type=1)

MeanDecreaseAccuracy

check                 15.1106884

duration               8.6281030

history                5.8685013

purpose                4.1078354

credit                 4.8800356

savings                6.1879918

employment             3.2416484

installment            3.0522029

personal               3.0602966

debtors                3.8963219

residence              1.1155901

property               4.7010080

age                    5.4753585

others                 3.3679530

housing                3.0737900

numcredits             0.6502566

job                    2.1420488

residpeople            0.1370845

telephone              0.6110284

foreign                0.5733125

결과 해석: check가 가장 중요하고 duration이 두번째로 중요

> order(importance(rf.german)[,'MeanDecreaseAccuracy'],decreasing = T)

[1]  1  2  6  3 13  5 12  4 10 14  7 15  9  8 17 11 16 19 20 18

> which.max(importance(rf.german)[,'MeanDecreaseAccuracy'])

check

1

> importance(rf.german)[which.max(importance(rf.german)[,'MeanDecreaseAccuracy']),]

bad                 good     MeanDecreaseAccuracy     MeanDecreaseGini

13.142687             9.632083                15.110688            44.838547

> rf.german$importance

bad          good MeanDecreaseAccuracy MeanDecreaseGini

check        5.728619e-02  2.488337e-02         3.461467e-02        44.838547

duration     1.262784e-02  2.059686e-02         1.820397e-02        41.273512

history      1.165482e-02  8.234252e-03         9.100115e-03        25.731461

purpose      1.001589e-02  5.877403e-03         7.125436e-03        35.651165

credit       8.957835e-03  9.469900e-03         9.342178e-03        52.641745

savings      1.873739e-02  4.572813e-03         8.795598e-03        21.705909

employment   5.377379e-03  4.391391e-03         4.647218e-03        23.910509

installment  5.083834e-05  4.492258e-03         3.109720e-03        15.714499

personal     4.903039e-03  2.481005e-03         3.188303e-03        15.365475

debtors     -1.380717e-04  3.198633e-03         2.202015e-03         6.718969

residence    1.833243e-03  1.342976e-03         1.400876e-03        17.932937

property     6.357328e-03  5.871030e-03         6.069060e-03        18.803008

age          1.026633e-02  8.969565e-03         9.239623e-03        41.451103

others       6.321219e-04  4.940365e-03         3.562374e-03        11.399935

housing      5.086080e-03  2.740150e-03         3.358222e-03         9.596539

numcredits  -3.793953e-04  7.653351e-04         4.256327e-04         7.944861

job          2.803065e-03  1.422938e-03         1.912685e-03        12.054190

residpeople  9.171770e-04 -2.801992e-04         7.286163e-05         4.234831

telephone    2.173474e-03 -2.108124e-04         4.743466e-04         6.427048

foreign      5.935867e-04 -1.538455e-05         1.547421e-04         1.519832

추가:

rf.german$importance #이건 importance(rf.german) 의 결과와 어떻게 다른지 좀 더 공부하자

3) 분류예측치를 구하고 정확도 평가

> pred.rf.german<-predict(rf.german,newdata=german)

> head(pred.rf.german,10)

  1    2    3    4    5    6    7    8    9   10

  good  bad good good  bad good good good good  bad

  Levels: bad good

> tab=table(german$y,pred.rf.german,dnn=c("Actual","Predicted"))

> tab

        Predicted

  Actual  bad   good

    bad    300    0

 good   0      700

> addmargins(tab)

  Predicted

  Actual  bad  good  Sum

  bad   300   0      300

  good  0      700  700

  Sum   300  700  1000

> 1-sum(diag(tab)/sum(tab))

  [1] 0

결과 해석: 기존 cart 분류나무 모형 오분류율은 19.6%, 배깅은 3.7%, 부스팅은 22.6%인 반면

랜덤 포레스트는 오분류율이 0%로 기존 그 어떠한 모형보다 대단히 우수한 결과를 보임을 알 수 있다

(분류나무 cart 모형 바로가기) / (분류앙상블 모형 – 배깅 바로가기) / (분류앙상블 모형 – 부스팅 바로가기)

4) 몇개의 분류기가 적당할까?

> plot(rf.german,'l')

결과 해석: x축은 분류기, y축은 OOB 오분류율을 의미 out of bag의 약자로 부트스트랩에 포함되지 못한 데이터를 의미. 결국 oob데이터는 검증데이터의 역할과 비슷. 가장 아래선은 목표변수 중 good에 대한 오분율, 가장 위쪽 선은 bad에 대한 오분율. 가운데는 전체 오분류율. 그림에 따르면 분류기 개수가 80개이상이면 안정적인 형태로 보인다

5) 훈련 데이터와 검증데이터로 분할하여 랜덤포레스트 방법을 평가해 보자.

> set.seed(1234)

> i=sample(1:nrow(german),round(nrow(german)*0.7)) #70% for training 훈련 data, 30% test 검증 데이

> german.train = german[i,]

> german.test = german[-i,]

> rf.train.german<-randomForest(y~.,data=german.train,ntree=100,mtry=5,importance=T,na.action = na.omit)

> pred.rf.train.german<-predict(rf.train.german,newdata=german.test)

> tab.train = table(german.test$y,pred.rf.train.german,dnn=c('Actual','Predicted'))

> tab.train

Predicted

Actual bad good

bad   39   60

good  15  186

> addmargins(tab.train)

Predicted

Actual   bad  good Sum

bad    39   60     99

good  15   186    201

Sum   54   246    300

> 1-sum(diag(tab.train) / sum(tab.train)) #오분류율이 25%로 어느정도 향상

[1] 0.25

결과 해석: 검증데이터에 대한 오분류율은 기존 cart 분류 나무는 26.3%, 배깅 역시 26.3%, 부스팅은 25.3%. 랜덤포레스트는 25%로, 조금 향상 되었음.

(분류나무 cart 모형 바로가기) / (분류앙상블 모형 – 배깅 바로가기) / (분류앙상블 모형 – 부스팅 바로가기)

출처: 데이터마이닝(장영재, 김현중, 조형준 공저) 

목표변수가 집단을 의미하는 범주형 의사결정나무 -> 분류나무모형

목표변수가 연속형 변수인 의사결정나무 -> 회귀나무모형

 목표 변수가 연속형 변수인 경우, 의사결정나무는 회귀나무라고 한다. 통계학에서 목표변수가 연속형일 때 흔히 수행하는 분석이 선형 회귀분석인 것을 연상해 본다면 명칭이 이해가 될 것이다. 일반적으로 범주형 변수들은 가변수화하여 0,1의 형태로 변형하여 회귀분석에 활용하지만, 범주형 변수의 수가 매우 많고, 각 범주형 변수의 개수도 많은 경우 해석은 어려워지고 분석도 복잡해 진다. 이러한 경우 회귀나무 모형을 사용하게 되면, 가변수를 생성할 필요 없이 범주형 입력 변수와 연속형 입력변수를 그대로 활용할 수 있게 되어 분석 및 그 해석이 용이해질 수 있다.

예제) 목표변수가 숫자인 보스턴하우징 데이터를 이용하여 cart 방법을 이용한 회귀나무 모형 구축

1) 데이터 읽기

> library(MASS)

> Boston$chas=factor(Boston$chas)

> Boston$rad=factor(Boston$rad)

> summary(Boston)

> Boston<-Boston[,-15]
> str(Boston)
'data.frame':   506 obs. of  14 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : Factor w/ 9 levels "1","2","3","4",..: 1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ black  : num  397 397 393 395 397 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ... 

2) cart 회귀나무모형 실행

> library(rpart)

> my.control<-rpart.control(xval=10,cp=0,minsplit=nrow(Boston)*0.05)

> fit.Boston=rpart(medv~.,data=Boston,method='anova',control=my.control)

> fit.Boston

n= 506

node), split, n, deviance, yval

* denotes terminal node

1) root 506 42716.30000 22.532810

   … 중간 생략…

  30) rad=2,7,8 10   140.04100 43.670000 *

     31) rad=1,3,4,5,24 17   135.72120 48.158820 *

함수 설명:

rpart.control(xval=10,cp=0,minsplit=nrow(Boston)*0.05) : 10-fold교차타당성 오분류율을 계산하고, 가지치기 방법에서 a값은 0이될때까지 순차적인 나무구조를 저장하도록 한다. 그리고 노드를 나누는 최소 자료이 수인 minsplit은 전체 데이터 수의 5% 수준에서 정의한다

3) 가지치기 수행 : ①최소 cp 값 찾기 

> printcp(fit.Boston)

Regression tree:
rpart(formula = medv ~ ., data = Boston, method = "anova", control = my.control)

Variables actually used in tree construction:
[1] age     crim    indus   lstat   nox     ptratio rad     rm    

Root node error: 42716/506 = 84.42

n= 506

           CP nsplit rel error  xerror     xstd
1  0.45274420      0   1.00000 1.00298 0.083114
2  0.17117244      1   0.54726 0.62762 0.058631
3  0.07165784      2   0.37608 0.43188 0.047879
4  0.03428819      3   0.30443 0.35161 0.043496
5  0.02661300      4   0.27014 0.32766 0.042850
6  0.01802372      5   0.24352 0.29534 0.041781
7  0.01348721      6   0.22550 0.27799 0.039217
8  0.01285085      7   0.21201 0.28299 0.039277
9  0.00844925      8   0.19916 0.26415 0.032335
10 0.00833821      9   0.19071 0.25667 0.031357
11 0.00726539     10   0.18238 0.25624 0.031407
12 0.00612633     11   0.17511 0.25343 0.031377
13 0.00480532     12   0.16898 0.24000 0.028645
14 0.00410785     13   0.16418 0.23897 0.027439
15 0.00394102     14   0.16007 0.23177 0.027282
16 0.00385379     15   0.15613 0.23177 0.027282
17 0.00223540     16   0.15228 0.22527 0.027268
18 0.00171691  17   0.15004 0.22441 0.027196
19 0.00153485     18   0.14832 0.22453 0.026980
20 0.00140981     19   0.14679 0.22585 0.026952
21 0.00135401     20   0.14538 0.22646 0.026946
22 0.00113725     22   0.14267 0.22674 0.026936
23 0.00098921     23   0.14153 0.22490 0.026948
24 0.00081924     24   0.14054 0.22674 0.027098
25 0.00074570     25   0.13972 0.22624 0.027097
26 0.00074096     27   0.13823 0.22671 0.027092
27 0.00040763     28   0.13749 0.22753 0.027088
28 0.00031464     29   0.13708 0.22749 0.027081
29 0.00019358     30   0.13677 0.22798 0.027098
30 0.00000000     31   0.13658 0.22816 0.027099

> names(fit.Boston)

[1] "frame"               "where"               "call"                "terms"              

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"            

> which.min(fit.Boston$cp[,4])
18
18
> which.min(fit.Boston$cptable[,'xerror'])
18
18

> fit.Boston$cp[18,]
          CP       nsplit    rel error       xerror         xstd
 0.001716908 17.000000000  0.150039975  0.224411843  0.027196436
> fit.Boston$cp[18]
[1] 0.001716908

> fit.Boston$cptable[which.min(fit.Boston$cptable[,'xerror']),]
          CP       nsplit    rel error       xerror         xstd
 0.001716908 17.000000000  0.150039975  0.224411843  0.027196436
> fit.Boston$cptable[which.min(fit.Boston$cptable[,'xerror'])]
[1] 0.001716908

> 0.22441 +0.027196
[1] 0.251606

결과 해석: 가장 작은 cp값은 18번째 0.001716908이나 가장 작은 xerror의 xstd 범위 내 즉, 0.22441 + 0.027196 = 0.251606 내 작은 cp 값 적용도 가능하다.

3) 가지치기 수행 : ②가장 적은 cp값을 찾았으니 prune 함수를 이용하여 가지기치를 수행하자. 

> fit.prune.Boston<-prune(fit.Boston,cp=0.00304027)

> print(fit.prune.Boston)

> names(fit.prune.Boston)

[1] "frame"               "where"               "call"                "terms"             

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"           

> fit.prune.Boston$variable.importance
              rm             lstat           indus             nox             age         ptratio             dis
23124.34862 16998.05183  5871.44959  5492.56562  5042.88012  4757.36213  4742.25575
            tax            crim              zn             rad           black            chas
 4155.77974  2258.47799  1760.50488  1390.73661   772.28149    11.40364

node), split, n, deviance, yval 는 각각 노드번호, 분할규칙, 해등 노드의 총 관칙치 수, 해당 노드의 분산, 해당 노드의 목표변수 예측치(회귀나무일 때는 평균값)의 순서로 정보를 제공하고 있고, 줄 마지막에 있는 * denotes terminal node는 * 표시가 최종노드임을 알려주고 있다.

fit.prune.Boston$variable.importance로 중요 변수도 알 수 있다. rm이 가장 중요하고 상대적으로 chas가 덜 중요한 변수이다.

4) cart 회귀나무 모형 그리기

> plot(fit.prune.Boston)

목표변수가 집단을 의미하는 범주형 의사결정나무 -> 분류나무모형

목표변수가 연속형 변수인 의사결정나무 -> 회귀나무모형

예제) 목표변수가 범주형 good, bad인 독일 신용평가 데이터를 이용하여 cart 방법을 이용한 의사결정나무 구축

1) 데이터 불러오기

> setwd('c:/Rwork')

> german<-read.table('germandata.txt',header=T)

> str(german)

'data.frame':        1000 obs. of  21 variables:

$ check      : Factor w/ 4 levels "A11","A12","A13",..: 1 2 4 1 1 4 4 2 4 2 ...

$ duration   : int  6 48 12 42 24 36 24 36 12 30 ...

$ history    : Factor w/ 5 levels "A30","A31","A32",..: 5 3 5 3 4 3 3 3 3 5 ...

$ purpose    : Factor w/ 10 levels "A40","A41","A410",..: 5 5 8 4 1 8 4 2 5 1 ...

$ credit     : int  1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...

$ savings    : Factor w/ 5 levels "A61","A62","A63",..: 5 1 1 1 1 5 3 1 4 1 ...

$ employment : Factor w/ 5 levels "A71","A72","A73",..: 5 3 4 4 3 3 5 3 4 1 ...

$ installment: int  4 2 2 2 3 2 3 2 2 4 ...

$ personal   : Factor w/ 4 levels "A91","A92","A93",..: 3 2 3 3 3 3 3 3 1 4 ...

$ debtors    : Factor w/ 3 levels "A101","A102",..: 1 1 1 3 1 1 1 1 1 1 ...

$ residence  : int  4 2 3 4 4 4 4 2 4 2 ...

$ property   : Factor w/ 4 levels "A121","A122",..: 1 1 1 2 4 4 2 3 1 3 ...

$ age        : int  67 22 49 45 53 35 53 35 61 28 ...

$ others     : Factor w/ 3 levels "A141","A142",..: 3 3 3 3 3 3 3 3 3 3 ...

$ housing    : Factor w/ 3 levels "A151","A152",..: 2 2 2 3 3 3 2 1 2 2 ...

$ numcredits : int  2 1 1 1 2 1 1 1 1 2 ...

$ job        : Factor w/ 4 levels "A171","A172",..: 3 3 2 3 3 2 3 4 2 4 ...

$ residpeople: int  1 1 2 2 2 2 1 1 1 1 ...

$ telephone  : Factor w/ 2 levels "A191","A192": 2 1 1 1 1 2 1 2 1 1 ...

$ foreign    : Factor w/ 2 levels "A201","A202": 1 1 1 1 1 1 1 1 1 1 ...

$ y          : Factor w/ 2 levels "bad","good": 2 1 2 2 1 2 2 2 2 1 ...

> german$numcredits<-factor(german$numcredits)

> german$residence<-factor(german$residence)

> german$residpeople<-factor(german$residpeople)

> class(german$numcredits);class(german$residence);class(german$residpeople)

[1] "factor"

[1] "factor"

[1] "factor"

2) cart 방법 적용

> library(rpart)

> my.control<-rpart.control(xval=10,cp=0,minsplit=5)

> fit.german<-rpart(y~.,data=german,method='class',control=my.control)

> fit.german #최초의 나무. 가지치기를 하지 않은 최대 크기의 나무 보기

n= 1000

node), split, n, loss, yval, (yprob)

* denotes terminal node

1) root 1000 300 good (0.300000000 0.700000000) 

2) check=A11,A12 543 240 good (0.441988950 0.558011050) 

.

. (너무 커서 중략)

.

    253) credit>=1273 10   0 good (0.000000000 1.000000000) *

  127) check=A14 122   1 good (0.008196721 0.991803279) *

함수 설명

 1. rpart.control:

- xval=10: 교타 타당성의 fold 개수, 디폴트는 10

- cp=0: 오분류율이 cp값 이상으로 향상되지 않으면 더 이상 분할하지 않고 나무구조 생성을 멈춘다. cp값이 0이면 오분류값이 최소, 디폴트는 0.01

- minsplit=5: 한 노드를 분할하기 위해 필요한 데이터의 개수. 이 값보다 적은 수의 관측치가 있는 노드는 분할하지 않는다. 디폴트는 20

2. r.part

 - method=class: 나무 모형을 지정한다. anova는 회귀나무, poisson 포아송 회귀나무, class는 분류나무 exp는 생존나무. 디폴트는 class

 - na.action=na.rpart: 목표변수가 결측치이면 전체 관측치를 삭제. 입력변수가 결측치인 경우에는 삭제하지 않는다.

결과 해석

중간노드를 분할하는 최소 자료의 수를 5개로 지정하였고, cp값은 0으로 하여 나무모형의 오분류값이 최소가 될 때 까지 분할을 진행하였다. 또한 10-fold 교차타당성을 수행하여 최적의 cp값을 찾도록 하였다. 나무가 너무나 큰 관계로 중간 부분을 생략하였고, 용이한 모형 분석을 위해 가지치기를 해보자.

3) 큰 나무를 줄이기 위한 가지치기 작업

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

Root node error: 300/1000 = 0.3

n= 1000

               CP nsplit    rel error  xerror       xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556      8   0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000    13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

결과 해석

10-fold 교차타당성 방법에 의한 오분율(xerror)이 최소가 되는 값은 0.85667이며 이때의 cp값은 0.01임을 알 수 있다. 이 때 분리의 횟수가 13회(nsplit=13)인 나무를 의미한다.

또는 아래와 같은 방법으로도 최소 오분류값(xerror)를 찾을 수 있다.

> names(fit.german)

[1] "frame"               "where"               "call"                "terms"             

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"           

> fit.german$cptable[,'xerror']

1         2         3         4         5         6         7         8         9

1.0000000 0.9466667 0.8633333 0.8700000 0.8866667 0.8800000 0.8566667 0.8700000 0.8733333

10        11        12        13        14        15        16        17

0.9200000 0.9100000 0.9200000 0.9200000 0.9700000 0.9766667 0.9766667 1.0133333

> which.min(fit.german$cptable[,'xerror'])

7

7

> fit.german$cptable[7,]

CP       nsplit     rel error       xerror        xstd

0.01000000 13.00000000  0.65333333  0.85666667  0.04606167

> fit.german$cptable[7]

[1] 0.01

> fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

[1] 0.01

> min.cp<-fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

> min.cp

[1] 0.01

> fit.prune.german<-prune(fit.german,cp=min.cp)

4) 오분류율이 최소인 cp값(=0.011)을 찾았으니 이 값을 기준으로 가지치기를 시행하자.

> fit.prune.german<-prune(fit.german,cp=0.01)

> fit.prune.german

결과 해석

node), split, n, loss, yval, (yprob) 기준으로 첫번째 결과를 분석하면 다음과 같다.

노드, 분할점, 개수, …공부 필요

16) duration>=47.5 36   5 bad (0.8611111 0.1388889) *

duration 변수 중 47.5보다 큰 경우, 전체 36(n)개를 bad(yval)로 분류하였고 그 중 5개(loss)가 good이다. 그리하여 bad로 분류되는 것은 31/36 = 0.8611111로 표기하게 되고, 5개의 loss는 5/36 = 1388889 로 그 확률을 볼 수 있다. 아래 plot에서는 bad 31/5로 표기

376) property=A123,A124 20   3 bad (0.8500000 0.1500000) *

377) property=A121,A122 45  14 good (0.3111111 0.6888889) *

property가 a123(car), a124(unknown / no property)의 경우 전체 20개를 bad로 분류하였고 3개의 loss 즉 good (3/20 = 0.15)로 분류하였다. 아래 plot에서는 bad 17/3로 표기 

property가 a121(real estate), a122(building society savings agreement)인 경우에는 전체 45개를 good으로 분류하였고 14개의 loss 즉 bad로 분류 (14/45=0.3111111), 아래 plot에서는 good 14/31로 표기

<< 17.6.18(일)>> 해석 부분 내용 추가

duration > = 22.5인 경우, 전체 고객은 237명이고, 이 중 신용도가 나쁜 사람의 비율은 56.5%이고 좋은 사람의 비율은43.5%로 103명이다. 따라서 duration > 22.5 그룹은 bad로 분류된다.

가지치기를 한 모형을 그림으로 나타내는 함수는 아래와 같다.

> plot(fit.prune.german,uniform = T,compress=T,margin=0.1)

> text(fit.prune.german,use.n=T,col='blue',cex=0.7)

왼쪽 가지의 가장 아랫부분의 분할점인 ‘purpose=acdeghj’는 purpose 변수의 범주값 중에서 알파벳 순서로, 1(=a), 3(=c), 4(=d), 5(=e), 7(=g), 8(=h), 10(=j)번째 범주값을 의미하며, fit.prune.german에서 각각 A40,A410,A42,A43,A45,A46,A49 임을 알 수 있다.

34) purpose=A40,A410,A42,A43,A45,A46,A49 137  52 bad (0.6204380 0.3795620) *

<< 17.6.18(일)>> 해석 부분 내용 추가

가장 우측의 duration > = 11.5가 아닌 경우, 신용다가 나쁜 / 좋은 사람의 비율은 9명 / . 4명이고, 신용도가 좋은 good으로 분류된다.

5) 나무수를 더 줄여보자.

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

Root node error: 300/1000 = 0.3

n= 1000

                CP nsplit   rel error  xerror      xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556      8   0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000     13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

5번째 단계이며 분리의 횟수가 8회(nsplit=8)인 나무는 교차타당성 오분류율이 0.88667로 최소는 아니지만 7번째 단계의 분리의 횟수 13회 나무 가지의 최소 오분류율 0.85667과는 크게 차이가 나지 않는다. 그리고 최소 오분류율 표준편차의 1배 범위(0.88667  <  0.85667 + 0.046062)에 있다. 이런 경우에는 5번째 단계이며 분리의 횟수가 8인 나무를 선택하는 경우도 있다.

5번째 단계이며 분리 횟수가 8인 cp값 0.0155556의 반올림 값 0.016 적용하여 다시 가지치기

> fit.prune.german<-prune(fit.german,cp=0.016)

> fit.prune.german

> plot(fit.prune.german,uniform=T,compress=T,margin=0.1)

> text(fit.prune.german,use.n=T,col='blue',cex=0.7)

6) 목표변수의 분류예측치를 구하고 그 정확도에 대해서 평가해 보자

 > fit.prune.german<-prune(fit.german,cp=0.01) 

 > pred.german=predict(fit.prune.german,newdata=german,type='class')

 > tab=table(german$y,pred.german,dnn=c('Actual','Predicted'))

 > tab

           Predicted

  Actual  bad  good

    bad  180   120

   good   76   624

함수 설명

predict(fit.prune.german,newdata=german,type='class'), type = class는 분류나무의 집단값 예측결과, 회귀나무라면 type = vector라고 해야 한다.

결과 해석

실제 good인데 good으로 예측한 것이 624개, 실제 bad인데 bad로 예측한 것이 180

따라서 오분류율은 {1000 – (624+180)} / 1000 = 19.6%

R코드를 이용하면 1-sum(diag(tab)) / sum(tab)

7) 마지막으로 독일신용평가데이터를 훈련데이터와 검증 데이터로 분할하여 분류나무를 평가해보자.

> set.seed(1234)

> i=sample(1:nrow(german),round(nrow(german)*0.7)) #70% for training훈련 data, 30% for test검증

> german.train=german[i,]

> german.test=german[-i,]

> fit.german<-rpart(y~.,data=german.train,method='class',control=my.control)

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german.train, method = "class",

        control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings     telephone 

Root node error: 201/700 = 0.28714

n= 700

CP  nsplit   rel error  xerror     xstd

1  0.05721393      0   1.00000 1.00000 0.059553

2  0.03482587      2   0.88557 1.00498 0.059641

3  0.02985075      5   0.78109 1.00000 0.059553

4  0.01990050    6   0.75124 0.95025 0.058631

5  0.01741294      8   0.71144 0.96020 0.058822

6  0.01492537     10   0.67662 1.00000 0.059553

7  0.01243781     14   0.61692 1.00000 0.059553

8  0.00995025     17   0.57711 1.00995 0.059728

9  0.00746269     35   0.39303 1.03980 0.060238

10 0.00621891     46   0.30846 1.06965 0.060722

11 0.00497512     50   0.28358 1.04975 0.060402

12 0.00331675     58   0.24378 1.09950 0.061181

13 0.00248756     61   0.23383 1.11940 0.061474

14 0.00124378     69   0.21393 1.14925 0.061894

15 0.00099502     73   0.20896 1.14925 0.061894

16 0.00000000     78   0.20398 1.14925 0.061894

> fit.prune.german<-prune(fit.german,cp=0.02)

> fit.prune.german

> p.german.test=predict(fit.prune.german,newdata=german.test,type='class')

> tab=table(german.test$y,p.german.test,dnn=c('Actual','Predicted'))

> tab

            Predicted

Actual  bad good

      bad   34   65

      good  14  187

> 1-sum(diag(tab))/sum(tab) #오분류율

[1] 0.2633333 

출처: 데이터마이닝(장영재, 김현중, 조형준 공저,knou press)

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0

보스턴 하우징 데이터 Housing Values in Suburbs of Boston

(출처: http://127.0.0.1:31866/library/MASS/html/Boston.html)

1. 데이터 불러오기

2. 선형 회귀 모형 만들기

-> 결과 해석: 최초 만든 회귀함수 fit1 = lm(medv~.,data=boston)에서, 가장 적절한 모형 선택 위한 변수 선택을 위해 step 함수를 사용한다. fit.step = step(fit1,direction='both') 이 함수에서 both는 단계적 선택법 적용을 의미한다. 결과적으로 lm(formula = medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio + black + lstat, data = boston)라는 모형이 만들어졌고, 최초 만든 모형 대비 indus, chas1 변수가 사라졌음을 알 수 있다. 또한 최종 모형에서 범주 rad2,3,4 등은 범주형 범수 중 특정 변수를 의미한다.

    입력변수 crim의 회귀계수 추정치는 음수이므로 crim이 증가함에 따라 목표변수medv는 감소한다. nox 변수의 회귀곗수는 -12인데, nox 변수가 올라갈 때 마다 medv 값은 내려간다. nox 변수는 10ppm당 농축 일산화질소를 뜻한다.

      rad변수는 9개 범주로 구성되어 있기 때문에 8개의 가변수가 생성되었다. 각 입력 변수의 t값의 절대값으 커서 대응하는 p-값은 0.05보다 작아서 유의하다고 할 수 있다. 단, rad6는 유의하지 않지만 다른 가변수가 유의하므로 제거되지 않고 여전히 모형에 포함된다. 

     R²은 79.07%로 적합한 선형 회귀모형으로 데이터를 설명할 수 있는 부분이 약 80%로 높고, F-검정의 p-value도 2.2e-16로 아주 작은 것도 모형이 적합하다는 것을 지지하다.  

3. 어떤 변수들이 제거 되었을까?

->  해석: fit.step$anova라는 함수를 통해 최종 모형에서 제거된 변수를 알 수 있다. 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택하게 되는데, 여기서는 chas와 indus가 제거되었음을 일목요연하게 알 수 있다.

4. 목표 예측값을 알아보자

         ->   함수 predict 는 다양한 모형 적합결과로부터 예측값을 계산할 때 사용하고, 이중 type 

           옵션은 예측 형태를 입력하는 것으로, 목표값을 예측할 때 ‘response’를 사용한다.

->   목표변수가 연속형인 경우에 모형의 예측력 측도로서 MSE(mean squared error)를 주로 사용한다. 관측치(yi)와 예측치 Ŷi의 차이가 적을수록 그 모형의 예측력은 높다고 할 수 있다. 이를 시각적으로 확인하기 위해서는 이들을 가로축 및 세로축에 놓고 그린 산점도가 45도 대각선을 중심으로 모여 있으면 예측력이 좋다고 할 수 있다. 

출처: 데이터마이닝(장영재, 김현중, 조형준 공저)