데이터마이너를 꿈꾸며

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

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

부스팅 방법 실행

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(rpart)

> library(adabag)

> my.control<-rpart.control(xval=0,cp=0,maxdepth=1)

> boo.german<-boosting(y~.,data=german,boos=T,mfianl=100,control=my.control)

> names(boo.german)

[1] "formula"    "trees"      "weights"    "votes"      "prob"       "class"      "importance"

[8] "terms"      "call"     

> summary(boo.german)

Length Class   Mode     

formula       3   formula call    

trees       100   -none-  list    

weights     100   -none-  numeric 

votes      2000   -none-  numeric 

prob       2000   -none-  numeric 

class      1000   -none-  character

importance   20   -none-  numeric 

terms         3   terms   call    

call          6   -none-  call

함수 설명: boos=T, 표본추출에 의한 분류기 생성방식, boos=F 가중치를 반영한 분류기 생성방식 사용, mfinal은 부스팅 방법에서 생성하게 될 분류기의 개수 B, 100개의 분류기는 rpart 분류나무를 이용하여 생성되고, 각 분류나무는 최대분할의 깊이는 1이라는 옵션, 즉 100개의 스텀프나무를 생성

2) 부스팅 방법의 실행 - 중요도 보기 

> names(boo.german)

[1] "formula"    "trees"      "weights"    "votes"      "prob"       "class"      "importance"

[8] "terms"      "call"

> head(boo.german$trees) #B=100개의 분류나무를 출력할 수 있다

[[1]]

n= 1000

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

* denotes terminal node

1) root 1000 294 good (0.2940000 0.7060000) 

2) check=A11 275 119 bad (0.5672727 0.4327273) *

3) check=A12,A13,A14 725 138 good (0.1903448 0.8096552) *

… 중간 생략 …

> boo.german$importance

age        check       credit      debtors     duration   employment      foreign

5.757338101 40.631489260  5.445313887  0.008302238 15.889422076  2.918286176  0.133462457

history      housing  installment          job   numcredits       others     personal

12.835057694  0.312497120  1.276787836  0.000000000  0.000000000  1.941629429  1.282855587

property      purpose    residence  residpeople      savings    telephone

1.402513542  6.311856846  0.537511661  0.000000000  3.315676089  0.000000000

> order(boo.german$importance,decreasing=T) #bagging 방법과 중요도 변수가 조금은 다르네.

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

> boo.german$importance[2]

check

40.63149

> which.max(boo.german$importance)

check

2

> boo.german$importance[which.max(boo.german$importance)]

check

40.63149

> boo.german$importance[5]

duration

15.88942

> boo.german$importance[8]

history

12.83506

> importanceplot(boo.german)

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

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

> names(pred.boo.german)

  [1] "formula"   "votes"     "prob"      "class"     "confusion" "error"   

> head(pred.boo.german$prob,10) 

  [,1]      [,2]

  [1,] 0.28278825 0.7172117

  [2,] 0.56170615 0.4382938

  [3,] 0.17385176 0.8261482

  [4,] 0.50498723 0.4950128

  [5,] 0.52209125 0.4779088

  [6,] 0.35035420 0.6496458

  [7,] 0.18447524 0.8155248

  [8,] 0.44024939 0.5597506

  [9,] 0.08690876 0.9130912

  [10,] 0.48660566 0.5133943

> pred.boo.german$confusion

  Observed Class

  Predicted Class bad good

           bad  135   61

           good 165  639

  > addmargins(pred.boo.german$confusion)

                 Observed Class

  Predicted Class  bad good  Sum

           bad   135   61  196

           good  165  639  804

  Sum  300  700 1000

  > 1-sum(diag(pred.boo.german$confusion)) / sum(pred.boo.german$confusion)

  [1] 0.226

함수 설명: > head(pred.boo.german$prob,10) #prob벡터는 부스팅 방법에서 집단별 투표비율을 출력

결과 해석: 기존 cart 분류나무 모형 오분류율은 19.6%, 배깅은 3.7%, 

부스팅의 오분류율은 22.6%로 기존 cart 분류나무 모형 대비 다소 결과가 떨어짐을 알 수 있다.

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

4) 분류 앙상블에서 몇 개의 분류기가 적당한 것인지를 알 아 볼 수 있다.

> evol.german <- errorevol(boo.german,newdata=german) 

> plot.errorevol(evol.german)

결과 해석: x축은 분류기의 개수, y축은 오분류율을 의미, 그림에 따르면 분류기의 개수가 30개가 넘어가면 비교적 안정적인 결과를 보인다. 하지만 90개 이상인 경우 약간 오분류율이 감소하는 경향도 있는 듯하다.따라서 독일신용평가데이터에는 배깅 앙상블의 크기를 최소 90개 이상으로 정하면 된다. 

5) 데이터의 70%를 랜덤하게 훈련 데이터로 할당하고, 나머지 30%는 검증 데이터로 할당

#데이터를 입력

> set.seed(1234)

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

> german.train = german[i,]

> german.test = german[-i,]

#부스팅 앙상블 분류 모형 만들기

> my.control<-rpart.control(xval=0,cp=0,maxdepth=1)

> boo.train.german<-boosting(y~.,data=german.train,boos=T,mfianl=100,control=my.control)

#중요한 변수 분석

> names(boo.train.german)

[1] "formula"    "trees"      "weights"    "votes"      "prob"       "class"      "importance"

[8] "terms"      "call"     

> boo.train.german$importance

age       check      credit     debtors    duration  employment     foreign     history

3.9363991  32.5300152   4.3733165   0.0000000  16.0942441   6.2892310   0.0000000  13.8030374

housing installment         job  numcredits      others    personal    property     purpose

0.3568900   2.0703206   0.0000000   0.6247335   1.1360575   0.7039699   0.0000000  12.6654775

residence residpeople     savings   telephone

0.0000000   0.0000000   5.4163077   0.0000000

> order(boo.train.german$importance,decreasing = T)

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

> boo.train.german$importance[2]

check

32.53002

> boo.train.german$importance[5]

duration

16.09424

> boo.train.german$importance[8]

history

13.80304

> importanceplot(boo.train.german)

#예측하기

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

> pred.boo.train.german$confusion

Observed Class

Predicted Class bad good

bad   42   19

good  57  182

> addmargins(pred.boo.train.german$confusion)

Observed Class

Predicted Class bad good Sum

bad   42   19  61

good  57  182 239

Sum   99  201 300

#오분류율 계산

> 1-sum(diag(pred.boo.train.german$confusion)) / sum(pred.boo.train.german$confusion)

[1] 0.2533333

결과 해석: 검증데이터에 대한 오분류율은 기존 cart 분류 나무는 26.33%, 배깅 역시 26.3%, 

              반면 부스팅은 25.3%으로 계산되어, 조금 향상되었음.

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

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

4.2 범주형 자료의 검정

#예 4.2) 비타민 c가 감기치료에 효과가 있는지 점검. 대조군(control) 그룹 140명에게는 플라시보를 주고 처리군(treat) 그룹 139명에게는 매일 비타민 c를 투여하였다. 아래 분할표를 가지고 비타민 C가 감기에 효과가 있는지 점검

H0 복용군과 비복용군의 감기 이환율 같다
H1 복용군과 비복용군의 감기 이환율 다르다

1) 자료 입력

> vitamin = matrix(c(31,109,17,122),nrow=2,byrow=T)

> dimnames(vitamin) = list(vitamin=c('ctr','trt'),flu=c('y','n'))

> vitamin

          flu

vitamin  y   n

    ctr 31 109

    trt 17 122

> round(vitamin/sum(vitamin),2)

            flu

vitamin    y    n

    ctr 0.11 0.39

    trt 0.06 0.44

> addmargins(vitamin)

          flu

vitamin  y   n Sum

    ctr 31 109 140

    trt 17 122 139

    Sum 48 231 279

2) 데이터 시각화

> par(mfrow=c(1,2))

> dotchart(vitamin)

> dotchart(t(vitamin))

> par(mfrow=c(1,1))

> mosaicplot(vitamin)

결과 해석: 비타민 복용군과 비복용군의 감기 이환율이 동일하지 않음을 알 수 있다.

3)카이제곱 검정 실행

#카이검정

> chisq.test(vitamin)

Pearson's Chi-squared test with Yates' continuity correction

data:  vitamin

X-squared = 4.1407, df = 1, p-value = 0.04186

결과해석
대립가설 H0 : 복용군과 비복용군의 감기 이환율 같다
귀무가설 H1 : 복용군과 비복용군의 감기 이환율 다르다

p-value : 0.04186

결정: p-value값이 0.05보다 작으므로 H0를 기각, 비타민 복용군과 비복용군 간 이환율은 다르다.

관찰도수, 기대도수, 잔차를 보는 법 

#관찰도수

> names(chisq.test(vitamin))

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed"  "expected"

[8] "residuals" "stdres"  

> chisq.test(vitamin)$observed

flu

vitamin  y   n

ctr 31 109

trt 17 122

#기대도수

> chisq.test(vitamin)$expected

flu

vitamin        y       n

ctr 24.08602 115.914

trt 23.91398 115.086

#피어슨잔차

> chisq.test(vitamin)$residual

flu

vitamin         y          n

ctr  1.408787 -0.6421849

trt -1.413846  0.6444908

출처: 보건정보데이터 분석(이태림 저)

제4장 범주형 자료의 분석

4.1 범주형 자료와 분할표

분할표 table 만들기 연습

> medi = read.table('c:/Rwork/medication.txt',header=T)

> head(medi,3)

癤퓆o medication surv

1     1      treat    y

2     2      treat    n

3     3      treat    y

열 이름이 깨져 보임.

> colnames(medi)<-c('no','medication','surv')

> head(medi,3)

no medication surv

1  1      treat    y

2  2      treat    n

3  3      treat    y

분할표 작성, treat = 처리군 , control = 대조군

> attach(medi)

> tab <- table(medication, surv)

> colnames(tab) = c('die','survival')

> rownames(tab) = c('trt','ctr')

> tab

surv

medication die survival

trt   1        4

ctr   3        2

> addmargins(tab)

surv

medication die survival Sum

trt   1        4   5

ctr   3        2   5

Sum   4        6  10

> tab/sum(tab)

surv

medication die survival

trt 0.1      0.4

ctr 0.3      0.2

> addmargins(tab/sum(tab))

surv

medication die survival Sum

trt 0.1      0.4 0.5

ctr 0.3      0.2 0.5

Sum 0.4      0.6 1.0

출처: 보건정보데이터 분석(이태림 저)

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

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

이번에는 배깅 방법을 이용하여 분류 앙상블 모형을 진행해보자

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)

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

[1] "factor"

[1] "factor"

[1] "factor"

2) 배깅 방법의 실행

> library(rpart)

> library(adabag)

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

> bag.german <- bagging(y~.,data = german, mfinal = 50, control = my.control)

> summary(bag.german)

Length Class   Mode    

formula        3  formula call    

trees         50  -none-  list    

votes       2000  -none-  numeric 

prob        2000  -none-  numeric 

class       1000  -none-  character

samples    50000  -none-  numeric 

importance    20  -none-  numeric 

terms          3  terms   call    

call           5  -none-  call

함수 설명

 rpart.control(xval = 0, cp = 0, minsplit = 5, maxdepth = 10) : mfinal은 배깅 방법에서 생성하게 될 분류기의 개수  B, 각 50개의 분류기는 rpart 분류나무를 이용하여 생성되고, 각각의 분류나무는 노드의 최소 데이터 수는 5이고 최대 분할의 깊이는 10이라는 옵션으로 생성된다.

3) 가장 중요한 변수가 뭔지 알아보자.

> names(bag.german)

[1] "formula"    "trees"      "votes"      "prob"       "class"      "samples"    "importance" "terms"    

[9] "call"     

> bag.german$importance

age        check       credit        debtors     duration      employment  foreign     history      housing

7.1242723 16.4701060  13.2520285  2.0382461  10.3813110   6.4053643   0.1100654   6.7854018   0.9387335

installment   job         numcredits   others      personal    property     purpose      residence    residpeople

1.7740456   1.8354714   0.9098175   2.6632713   3.0909168   3.8218936  11.4089057   3.5700495   0.3570470

savings      telephone

6.5166935   0.5463590

결과 해석: check 변수가 가장 중요한 입력 변수이고, credit이 두 번째로 중요한 역할을 하는 입력변수이다.

> order(bag.german$importance)

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

> order(bag.german$importance,decreasing=T)

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

> bag.german$importance[2]

check

16.47011

> bag.german$importance[3]

credit

13.25203

> which.max(bag.german$importance)

check

2

> which.min(bag.german$importance)

foreign

7

> bag.german$importance[which.max(bag.german$importance)]

check

16.47011

> bag.german$importance[which.min(bag.german$importance)]

foreign

0.1100654

> importanceplot(bag.german)

4) 목표변수의 분류예측치를 구하고 그 정확도에 대해 평가하는 방법에 대해 알아보자

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

> names(pred.bag.german)

  [1] "formula"   "votes"     "prob"      "class"     "confusion" "error"   

> head(pred.bag.german$prob,10) #각 집단 별 투표비율

      [,1] [,2]

  [1,] 0.04 0.96

  [2,] 0.92 0.08

  [3,] 0.02 0.98

  [4,] 0.30 0.70

  [5,] 0.72 0.28

  [6,] 0.16 0.84

  [7,] 0.00 1.00

  [8,] 0.16 0.84

  [9,] 0.00 1.00

  [10,] 0.86 0.14

  > head(pred.bag.german$class,10) #관측치마다 최종 예측 집단 출력

  [1] "good" "bad"  "good" "good" "bad"  "good" "good" "good" "good" "bad"

  > pred.bag.german$confusion #실제 목표변수의 값과 예측 목표변수의 값이 어느정도 유사한지 행렬의 형태로 보이고 있다.

                  Observed Class

  Predicted Class bad good

               bad  268    5

              good  32  695

> addmargins(pred.bag.german$confusion)

Observed Class

Predicted Class  bad good  Sum

bad   268    5  273

good   32  695  727

Sum   300  700 1000

> sum(pred.bag.german$confusion)

[1] 1000

> diag(pred.bag.german$confusion)

bad good

268  695

> 1-sum(diag(pred.bag.german$confusion)) /sum(pred.bag.german$confusion)

[1] 0.037

결과 해석: 오분류율이 3.7%로 기존 19.6% 보다 배깅은 대단히 우수한 결과를 보임을 알 수 있다

(분류나무 cart 모형 바로가기)

5) 분류 앙상블에서 몇 개의 분류기가 적당한 것인지를 알 아 볼 수 있다.

> evol.german <- errorevol(bag.german,newdata=german)

> plot.errorevol(evol.german)

결과 해석: x축은 분류기의 개수, y축은 오분류율을 의미, 그림에 따르면 분류기의 개수가 40개가 넘어가면 비교적 안정적인 결과를 보인다. 따라서 독일신용평가데이터에는 배깅 앙상블의 크기를 최소 40개 이상으로 정하면 된다.

6) 훈련데이터와 검증 데이터로 분할하여 배깅 방법으로 평가

> set.seed(1234)

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

> german.train = german[i,]

> german.test = german[-i,]

> bag.train.german <- bagging(y~., data = german.train, mfinal = 50, control = my.control)

> bag.train.german$importance

age        check        credit       debtors      duration    employment     foreign     history

5.9281898  14.8985632  12.2570321   2.1746067  13.6936647   7.2300214   0.0000000   6.3576993

Housing  installment        job    numcredits      others     personal    property     purpose

1.0506904   2.3547243   1.4585356   0.4620406   1.9891509   3.0157372   3.8454948  12.3522852

Residence  residpeople     savings   telephone

3.7829112   0.4905634   5.8943300   0.7637592

> order(bag.train.german$importance,decreasing=T)

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

> bag.train.german$importance[2]

check

14.89856

> bag.train.german$importance[5]

duration

13.69366

> which.max(bag.train.german$importance)

check

2

> which.min(bag.train.german$importance)

foreign

7

> importanceplot(bag.train.german)

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

> pred.bag.train.german$confusion

Observed Class

Predicted Class bad good

bad   44   24

good  55  177

> addmargins(pred.bag.train.german$confusion)

Observed Class

Predicted Class bad  good Sum

bad    44   24   68

good   55  177  232

Sum    99  201  300

> 1-sum(diag(pred.bag.train.german$confusion)) / sum(pred.bag.train.german$confusion)

[1] 0.2633333

결과 해석: 검증데이터에 대한 오분류율은 26.33%로 계산되어, 기존 분류 나무의 검증데이터 오분류율 26.33%과 동일하네… 왜 동일하지..조금이라도 향상되었어야 하는데.. 더 공부하자.

(기존 cart 분류나무 모형 바로가기)

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

보스턴 하우징 데이터 – 랜덤포레스트 방법의 회귀앙상블 모형

1) 데이터 읽기

> Boston$chas=factor(Boston$chas)

> Boston$rad=factor(Boston$rad)

> 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 ...

Factor 함수를 이용하여 범주형으로 변경하고, medv 변수를 목표변수로 다른 변수를 입력변수로 사용한다.

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

>library(randomForest)

> rf.boston<-randomForest(medv~.,data=Boston,ntree=100,mtry=5,importance=T,na.action=na.omit)

> rf.boston

Call:

  randomForest(formula = medv ~ ., data = Boston, ntree = 100,      mtry = 5, importance = T, na.action = na.omit)

Type of random forest: regression

Number of trees: 100

No. of variables tried at each split: 5

Mean of squared residuals: 9.743395

% Var explained: 88.46

> summary(rf.boston)

Length Class  Mode    

call              7    -none- call    

type              1    -none- character

predicted       506    -none- numeric 

mse             100    -none- numeric 

rsq             100    -none- numeric 

oob.times       506    -none- numeric 

importance       26    -none- numeric 

importanceSD     13    -none- numeric 

localImportance   0    -none- NULL    

proximity         0    -none- NULL    

ntree             1    -none- numeric 

mtry              1    -none- numeric 

forest           11    -none- list    

coefs             0    -none- NULL    

y               506    -none- numeric 

test              0    -none- NULL    

inbag             0    -none- NULL    

terms             3    terms  call

 함수 설명

ntree=100, 분류기 개수. 디폴트는 500개

mtry=5, 중간노드마다 랜덤하게 선택되는 변수들의 개수. 디폴트는 분류나무의 경우 sqrt(p), 회귀나무의 경우 p/3

importance=T, 변수의 중요도 계산, 디폴트는 F

na.action=na.omit, 결측치를 처리하는 방법, 변수의 중요도를 계산하게 하고, 결측치는 필요한 경우에만 삭제.

names()함수를 통해 rf.boston에 저장된 오브젝트의 리스트를 불러내어, $predicted를 이용하여 훈련 데이터의 예측 집단을 출력할 수 있다.

> names(rf.boston)

[1] "call"            "type"            "predicted"       "mse"             "rsq"           

[6] "oob.times"       "importance"      "importanceSD"    "localImportance" "proximity"     

[11] "ntree"           "mtry"            "forest"          "coefs"           "y"              

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

> head(rf.boston$predicted,30)

1        2        3        4        5        6        7        8        9       10

28.25382 22.55963 34.14192 35.45333 34.06798 26.81151 21.01950 16.78839 17.80599 19.23591

11       12       13       14       15       16       17       18       19       20

21.02440 21.23466 21.80889 20.05162 19.30557 20.21721 21.61349 18.46000 18.14724 19.96174

21       22       23       24       25       26       27       28       29       30

14.10136 18.55984 16.05801 15.04825 16.70996 15.70548 17.84748 14.82048 18.88633 20.64939

importance()함수를 통해 계산된 입력변수의 중요도를 알 수 있다.

> importance(rf.boston,type=1)

%IncMSE

crim     8.325232

zn       2.061869

indus    5.130483

chas     1.030915

nox      9.211906

rm      17.090802

age      5.229782

dis      8.322716

rad      5.342500

tax      4.604745

ptratio  7.102056

black    5.292651

lstat   14.652271

결과를 보면 rm변수과 lstat 변수의 중요도가 가장 높음을 알 수 있다.

함수 설명

 type=1,은 정분류율의 평균감소값, 2는 불순도의 평균감소값을 이용하여 계산

목표변수의 적합값을 구하고 평가하기 위해 평균오차제곱합(mse)를 계산.

> names(Boston)

[1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"     "dis"     "rad"   

[10] "tax"     "ptratio" "black"   "lstat"   "medv"  

> Boston$medv.hat = predict(rf.boston,newdata=Boston)

> mean((Boston$medv-Boston$medv.hat)^2) #mean square error(mse)

[1] 1.915207

기존 선형회귀 한 회귀 분류 나무 모형 결과의 평균오차제곱합 mean((Boston$medv-Boston$medv.hat)^2) = 10.8643  대비 랜덤포레스트의 평균오차제곱합이 1.915207로 설명력이 매우 증가되었음을 알 수 있다. 랜덤포레스트의 경우 부트스트랩을 이용하기 때문에 확률임의추출에 의한 변동성이 있을 수 있다. 따라서 모델링을 할 때 마다 결과가 다르기 때문에, 랜덤포레스트를 수차례 반복 시행하고 예측결과의 평균값을 취하는 경우도 있다.

(기존 보스턴 하우징 데이터 회귀나무모형 사례 분석 링크 바로가기)

랜덤 포레스트 회귀앙상블의 적합값과 실제값의 일치도를 보자. 예측일치도가 우수함을 알 수 있다.

> plot(Boston$medv,Boston$medv.hat,xlab='Observed Values',ylab='Fitted Values')

> abline(0,1)

기존 분류 회귀의 나무모형의 적합값과 실제값의 일치도와 비교해봐도 매우 우수함을 알 수 있다.

이 의사결정 나무를 활용하여 30% 검증 데이터에 적용시켜서 분류예측치를 구해보자. 그리고 그 예측치를 구해보자.

> set.seed(1234)

> nrow(Boston)

[1] 506

> i=sample(1:nrow(Boston),round(nrow(Boston)*0.7))

> Boston.train = Boston[i,] #70% for training data 훈련 데이터

> Boston.test = Boston[-i,] #30% for test data 검증 데이터

> rf.train.boston<-randomForest(medv~.,data=Boston.train,ntree=100,importance=T,na.action=na.omit)

#obtain the predicted values

> medv.hat.test<-predict(rf.train.boston,newdata=Boston.test)

> mean((Boston.test$medv-medv.hat.test)^2) #predicted mean square error

[1] 4.114596

검증 데이터에 대한 평균오차제곱합 mse는 4.11로 계산되었다. 기존 회귀 분류 나무모형 의 검증 데이터 오분류율 13.95258와 비교해보면 상당히 향상된 결과임을 알 수 있다.

(기존 보스턴 하우징 데이터 회귀나무모형 사례 분석 링크 바로가기)

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

앙상블(ensemble)모형이란 주어진 데이터를 이용하여 여러 개의 서로 다른 예측 모형을 생성한 후, 이러한 예측 모형의 예측 결과를 종합하여 하나의 최종 예측결과를 도출해 내는 방법을 말한다. 목표변수의 형태에 따라 분류분석에도 사용 가능하고, 회귀분석에도 사용 가능하다. 분류분석에 사용한다면 분류앙상블, 회귀분석에 사용한다면 회귀앙상블이라 부를 수 있다. 현실적으로 앙상블 모형은 대부분 분류모형에서 사용되고 있는 실정이다. 이유는 데이터마이닝의 영역에서 더 자주 필요로 하는 모형이 분류모형이기 때문이라고 추측된다.

데이터를 이용하여 생성해 낸 한 분류모형의 결과를 분류기(classifier)라 하자. 예측집단을 종합하는 방법으로는 주도 다수결 방식이 사용되고 있다. 다수결 방식에 따라 아래와 같이 구분할 수 있다.

① 단순 다수결 방식: 만약 예측치 중에서 6개의 분류기가 1이라고 예측하고, 5개의 분류기가 0이라고 예측했다면, 다수결 방식에 의해서 이 관찰치는 1이라고 최종 결론을 내린다. 배깅, 랜덤포레스트 방법이 단순 다수결 방식을 사용한다.

② 가중 다수결 방식: 각 분류기마다 가중치인 w_i를 고려해야 한다. w_i는 각 분류기 오류율의 역수 개념이다. 성능이 우수한 분류기에 가중치를 더 부여하는 것이다. 부스팅 방법이 가중다수결 방식을 사용한다.

앙상블 모형의 종류에 따른 구분은 다음과 같다.

① 배깅 방법: 배깅(bagging) 방법은 Breiman(1996)에 의해 개발된 분류 앙상블 방법이다. Bagging은 bootstrap aggregating의 약어로 훈련 데이터로부터 부트스트랩 데이터를 B번 생성하여 부트스트랩 데이터마다 분류기를 생성한 후 그 예측결과를 앙상블하는 방법이다. 배깅 방법은 불안정한 분류방법의 예측력을 획기적으로 향상시킨다고 알려져 있다. 분류나무 중에서 가지치기를 사용하지 않은 최대나무가 더 불안정한 불류방법이기 때문에 그 효과가 더욱 좋다.

② 부스팅 방법: 부스팅(boosting) 방법은 Freund and Schapire(1997)에 의해 개발된 분류 앙상블 방법이다. 부스팅에 사용되는 분류기는 오분류율을 랜덤하게 예측하는 것보다 조금이라도 좋은 예측모형이기만 하면 효과가 있다고 알려져 있다. 이는 예측력이 약한 분류 모형을 결합하여 강한 예측모형을 만드는 과정으로, 가장 많이 실행되는 알고리즘은 아다부스트(AdaBoost: adaptive boosing)방법이다. 그 방법으로는 두 가지 방식으로 수행할 수 있다.

 1) 가중치를 반영한 분류기 생성 방식

 2) 표본추출에 의한 분류기 생성방식

③ 랜덤포레스트 방법: randor forest 방법은 부트스트랩을 이용한 데이터의 변화 및 나무모형 분할방법에 랜덤성을 가미하여 나무 모형이 배깅과 부스팅보다 훨씬 더 다양해지도록 유도하는 아이디어를 가지고 있고, 더 정확한 예측력을 가지고 있다고 알려져 있다. 이 방법은 Breiman(2001)이 고안한 방법으로, 입력변수의 개수가 많을 때 그 효과가 극대화된다.

아래 나무모형과 앙상블 모형을 비교한 table, 각 방법을 클릭하면 해당 글로 이동

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

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

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

3.2 여러 집단의 비교
3.2.1 1개의 요인을 고려하는 경우
p.86, 자폐아, 정상아, 지진아에 대한 혈청 항원 농도에 대해 조사를 하였다. 
이들 사이에 면역 이상에 대한 차이가 있다고 할 수 있는가?
귀무가설 H0: u1 = u2 = u3, 대립가설 H1: not H0

1) 데이터 입력

> a<-c(755,365,820,900,170,300,325,385,380,215,400,343,415,345,410,460,225,440,400,360,435,450,360)

> b<-c(165,390,290,435,235,345,320,330,205,375,345,305,220,270,355,360,335,305,325,245,285,370,345,345,230,370,285,315,195,270,305,375,220)

> c<-c(380,510,315,565,715,380,390,245,155,335,295,200,105,105,245)

> boxplot(a,b,c,col='yellow',names=c('자폐아','정상아','지진아'))

> library(vioplot)

> vioplot(a,b,c,col='yellow',names=c('자폐아','정상아','지진아'))

2) 각 그룹의 평균과 분산

> sapply(list(a,b,c),mean)

[1] 419.9130 305.0000 329.3333

> sapply(list(a,b,c),var)

[1] 31693.356  4071.875 29224.524

3) 등분산 검정

> sera = c(a,b,c)

  > group = factor(rep(1:3,c(length(a),length(b),length(c))))

  > fligner.test(sera~group)   

  Fligner-Killeen test of homogeneity of variances

  data:  sera by group

  Fligner-Killeen:med chi-squared = 6.8506, df = 2, p-value = 0.03254

결과 해석: p-값이 0.03254이어서 등분산성에 조금은 문제가 있음을 알 수 있다.

4) one way Anova

> out = aov(sera~group)

> out

  Call:

    aov(formula = sera ~ group)

  Terms:

                       group Residuals

  Sum of Squares   185159.3 1236697.2

  Deg. of Freedom         2        68

  Residual standard error: 134.8582

  Estimated effects may be unbalanced

 > summary(out)

              Df   Sum Sq Mean Sq F value  Pr(>F)  

  group        2  185159   92580   5.091 0.00871 **

  Residuals   68 1236697   18187                  

  ---

    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석:

질문:이들 사이에 면역 이상에 대한 차이가 있다고 할 수 있는가?
귀무가설 H0: u1 = u2 = u3,

대립가설 H1: not H0

p – value: 0.00871

결정: 귀무가설을 기각, 세 그룹 사이 면역 이상에 대한 차이가 있다.

5) 모형 적합성 검토 = 오차검토

> par(mfrow=c(2,2))

> plot(out)

결과 해석: 정규성 분포에 약간 문제가 있지만 큰 문제는 아니다.

6) average profile plot 평균 반응 프로파일 그림 – 효과 크기를 알 수 있는 plot

> plot.design(sera~group)

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다.

7) 다중 비교: 어느 그룹 간 차이가 있는지 보자.

> tukey = TukeyHSD(out)

> tukey

  Tukey multiple comparisons of means

  95% family-wise confidence level

  Fit: aov(formula = sera ~ group)

  $group

             diff        lwr       upr     p adj

  2-1 -114.91304 -202.68435 -27.14174 0.0070326

  3-1  -90.57971 -197.82092  16.66150 0.1142305

  3-2   24.33333  -76.28971 124.95638 0.8315432

 > plot(tukey)

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다. 그룹1과 2의 차이만이 신뢰도구간을 0을 포함하지 안으므로 유의미하게 다르다고 결론을 내릴 수 있다.

8) LSD 최소유의차검정 test

> pairwise.t.test(sera,group)

Pairwise comparisons using t tests with pooled SD

data:  sera and group

1      2    

2 0.0076 -    

3 0.0938 0.5642

P value adjustment method: holm

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다.

3.2.2 2개의 요인을 고려하는 경우
(1) 반복이 없을 때

예제) 장비 사용에 대한 3가지 방법을 연령별로 다르게 교육. 숙지 시간이 연령, 방법에 따라 다른가?

귀무가설 h0: u1 = u2 = u3, 대립가설 h1: not h0

1. 데이터 읽기

> setwd('c:/Rwork')

> data=read.table('device.txt',header=T)

> head(data)

ages way hour

1 under20   A    7

2   20~29   A    8

3   30~39   A    9

4   40~49   A   10

5 above50   A   11

6 under20   B    9

> tail(data)

ages way hour

10 above50   B   12

11 under20   C   10

12   20~29   C   10

13   30~39   C   12

14   40~49   C   12

15 above50   C   14

2. two way ANOVA

> out = aov(hour~ages+way,data=data)

> out

Call:

  aov(formula = hour ~ ages + way, data = data)

Terms:

                      ages        way Residuals

Sum of Squares  24.933333 18.533333  3.466667

Deg. of Freedom         4         2         8

Residual standard error: 0.6582806

Estimated effects may be unbalanced

> summary(out)

Df Sum Sq Mean Sq F value   Pr(>F)   

ages         4 24.933   6.233   14.38 0.001002 **

way          2 18.533   9.267   21.39 0.000617 ***

Residuals    8  3.467   0.433                     

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: p value가 0.05보다 적으므로 H0를 기각, h1를 받아들인다. 숙지 시간은 연령, 방법에 따라 서로 유의하게 다르다.

3. 모형적합성 검토 = 오차검토

> par(mfrow = c(2,2))

> plot(out)

결과 해석: 오차의 등분산성 및 정규성에 문제가 없음을 알 수 있다.

4. 다중비교, 왜 서로 유의한 차이가 났을까?

4.1) 나이 별 보기

> attach(data)

> pairwise.t.test(hour,ages)

Pairwise comparisons using t tests with pooled SD

data:  hour and ages

20~29 30~39 40~49 above50

30~39     1.00     -     -     -     

40~49     1.00  1.00     -     -     

above50   0.18  0.66   0.91   -     

under20   1.00  1.00   1.00  0.13  

P value adjustment method: holm

결과해석: 50대 이상이 다른 나이 수준보다 높았다. 결국 나이 수준이 숙지시간에 차이를 보인 것은 50대 이상이 다른 나이의 수준과 차이가 났기 때문이다.

4.2) 교육 방법 별 보기

> pairwise.t.test(hour,way)

Pairwise comparisons using t tests with pooled SD

data:  hour and way

A     B   

B 0.549 -   

C 0.061 0.125

P value adjustment method: holm

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: 방법은 A와 C간 차이가 났다.

위에서 본 두 함수 pairwise.t.test(hour,ages) ,pairwise.t.test(hour,way) 에 대하여 그래프를 그리면 아래와 같다.

>par(mfrow=c(1,2))

> plot.design(hour~ages)

> plot.design(hour~way)

5. 다중 비교

> tukey = TukeyHSD(out)

> tukey

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = hour ~ ages + way, data = data)

$ages                  diff        lwr        upr     p adj

30~39-20~29      1.0000000 -0.8568723  2.8568723 0.4057524

40~49-20~29      1.3333333 -0.5235390  3.1902056 0.1877558

above50-20~29    3.3333333  1.4764610  5.1902056 0.0017351

under20-20~29   -0.3333333 -2.1902056  1.5235390 0.9676094

40~49-30~39      0.3333333 -1.5235390  2.1902056 0.9676094

above50-30~39    2.3333333  0.4764610  4.1902056 0.0154324

under20-30~39   -1.3333333 -3.1902056  0.5235390 0.1877558

above50-40~49    2.0000000  0.1431277  3.8568723 0.0348816

under20-40~49   -1.6666667 -3.5235390  0.1902056 0.0810838

under20-above50 -3.6666667 -5.5235390 -1.8097944 0.0009146

$way diff       lwr      upr      p adj

B-A  0.6 -0.5896489 1.789649 0.3666717

C-A  2.6  1.4103511 3.789649 0.0006358

C-B  2.0  0.8103511 3.189649 0.0034083

> plot(tukey) 

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: 그래프에서도 나이가 20대 미만과 50대 이상에서, 방법은 3번과 1번 그리고 3번과 2번에서 신뢰구간을 0을 포함하지 않으므로 유의한 차이가 있었음을 알 수 있다.

(2) 반복이 있을 때
예) 세 종류의 호르몬 처리와 성별에 따라 혈액 칼슘값에 차이가 있는지 알아보기 위해 남녀 각 15명씩을 선정하여 이들을 세 그룹으로 나누어 세 가지 호르몬 처리를 한 후 혈액 칼슘을 측정하였다.
성별에 따라 혈액 칼슘에 차이가 있는가? 처리와 성별에 대한 교호작용이 존재하는가?

H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

1. 데이터 입력

> data=read.csv('calcium.csv')

> head(data)

sex way   cal

1   M   A 16.87

2   M   A 16.18

3   M   A 17.12

4   M   A 16.83

5   M   A 17.19

6   F   A 15.86

> tail(data)

sex way   cal

25   M   C 24.46

26   F   C 30.54

27   F   C 32.41

28   F   C 28.97

29   F   C 28.46

30   F   C 29.65

2. two way anova

> out = aov(cal~sex*way,data=data)

> out

Call:

  aov(formula = cal ~ sex * way, data = data)

Terms:

  sex       way   sex:way Residuals

Sum of Squares     4.0627 1146.6420    3.8454   76.2924

Deg. of Freedom         1         2         2        24

Residual standard error: 1.782933

Estimated effects may be unbalanced

> summary(out

Df Sum Sq Mean Sq F value   Pr(>F)   

sex          1    4.1     4.1   1.278    0.269   

way          2 1146.6   573.3 180.355 3.47e-15 ***

sex:way      2    3.8     1.9   0.605    0.554   

Residuals   24   76.3     3.2                    

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석: 처리에 대한 p-value가 0.0001보다 적게 나와 처리수준 간 모평균 차이가 없다라는 귀무가설을 기각. 성별과 처리도 p value가 0.05를 넘어 교호작용은 없다.

3. 모형적합성 검토 = 오차검토

> par(mfrow=c(2,2))

> plot(out)

결과 해석: 모형적합성 검토, 잔차도를 그려본 결과 오차의 등분산성에 약간의 문제는 있으나 큰 문제는 없음

4. 교호작용 검토 

> par(mfrow=c(1,1))

> with(data,interaction.plot(sex,way,cal))

결과 해석: with(data,interaction.plot(sex,way,cal)) #두 개의 선이 비슷한 거리를 유지하면서 평행에 가까우므로 interaction  교호작용이 없음을 알 수 있다. interaction.plot은 두 그룹변수의 조합으로 y의 평균을 그래프에 넣어 두 그룹 변수가 서로 y의 평균에 영향을 주는지 보는 방법

5. 다중비교

> attach(data)

The following object is masked _by_ .GlobalEnv:

  sex

The following objects are masked from data (pos = 3):

  cal, sex, way

The following objects are masked from data (pos = 4):

  cal, sex, way

> pairwise.t.test(cal,sex)

Error in tapply(x, g, mean, na.rm = TRUE) :

  arguments must have same length

왜 오류나는지 모르겠다. 더 공부 필요.

 > pairwise.t.test(cal,way)

Pairwise comparisons using t tests with pooled SD

data:  cal and way

A       B     

B 0.052   -     

C 8.4e-16 1.2e-14

P value adjustment method: holm

결과 해석:C와 A, 그리고 C와 B간 방법이 유의하게 차이가 났다.

> tukey = TukeyHSD(out)

> tukey

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = cal ~ sex * way, data = data)

$sex

diff        lwr     upr    p adj

M-F 0.736 -0.6076702 2.07967 0.269434

$way

diff        lwr       upr     p adj

B-A  1.609 -0.3822165  3.600217 0.1295236

C-A 13.845 11.8537835 15.836217 0.0000000

C-B 12.236 10.2447835 14.227217 0.0000000

$`sex:way`

diff        lwr       upr     p adj

M:A-F:A  1.548 -1.9385413  5.034541 0.7421633

F:B-F:A  1.956 -1.5305413  5.442541 0.5236718

M:B-F:A  2.810 -0.6765413  6.296541 0.1661169

F:C-F:A 14.716 11.2294587 18.202541 0.0000000

M:C-F:A 14.522 11.0354587 18.008541 0.0000000

F:B-M:A  0.408 -3.0785413  3.894541 0.9990770

M:B-M:A  1.262 -2.2245413  4.748541 0.8686490

F:C-M:A 13.168  9.6814587 16.654541 0.0000000

M:C-M:A 12.974  9.4874587 16.460541 0.0000000

M:B-F:B  0.854 -2.6325413  4.340541 0.9720701

F:C-F:B 12.760  9.2734587 16.246541 0.0000000

M:C-F:B 12.566  9.0794587 16.052541 0.0000000

F:C-M:B 11.906  8.4194587 15.392541 0.0000000

M:C-M:B 11.712  8.2254587 15.198541 0.0000000

M:C-F:C -0.194 -3.6805413  3.292541 0.9999760

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 성별간에는 유의한 차이는 없지만 방법에는 유의한 차이가 났다. C와 A, 그리고 C와 B간 방법이 유의하게 차이가 났다.

> par(mfrow=c(2,2))

> plot(tukey)

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 처리 C와 A, C와 B간 유의하게 서로 달랐다.

<또 다른 방법> 

위의 R 코드를 다른 방법으로 해보면 아래와 같다.(출처: R을 이용한 통계 분석, 안재형 지음)

> boxplot(cal~way+sex,col='red',data=data)

교호작용이 있는지 본 후

> with(data,interaction.plot(sex,way,cal))

결과 해석: 두개의 선이 서로 만나지 않으므로 교호작용이 존재하지 않는다는 것을 알 수 있다. (교호작용: 두 그룹 변수가 서로 y의 평균에 영향을 주는지 보는 방법)

분산분석표를 구한다. 교호작용이 존재하면 곱하기(sex*way), 존재하지 않으면 더하기(sex+way)

> out2=lm(cal~sex+way,data=data)

> anova(out2)

Analysis of Variance Table

Response: cal

Df  Sum Sq Mean Sq  F value    Pr(>F)   

sex        1    4.06    4.06   1.3181    0.2614   

way        2 1146.64  573.32 186.0089 3.944e-16 ***

Residuals 26   80.14    3.08                      

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(out2)

Call:

  lm(formula = cal ~ sex + way, data = data)

Residuals:

  Min      1Q  Median      3Q     Max

-5.8170 -0.5815 -0.0335  0.6623  4.3730

Coefficients:

 Estimate Std. Error t value Pr(>|t|)   

(Intercept)  15.6960     0.6411  24.484  < 2e-16 ***

sexM          0.7360     0.6411   1.148   0.2614   

wayB          1.6090     0.7851   2.049   0.0506 . 

wayC         13.8450     0.7851  17.634 5.53e-16 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.756 on 26 degrees of freedom

Multiple R-squared:  0.9349,   Adjusted R-squared:  0.9274

F-statistic: 124.4 on 3 and 26 DF,  p-value: 1.532e-15

결과 해석:

-   sexM(M-F, p-value는 0.2614)의 추정치가 0.7360으로 유의하지는 않다.

-   wayC(C-A, p-value는 5.533-16)의 추청치가 13.8450으로 유의하고 평균은 A보다 높다.

다중비교

> library(multcomp)

> tukey2=glht(out2,linfct=mcp(way='Tukey'))

> tukey2

General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Linear Hypotheses:

  Estimate

B - A == 0    1.609

C - A == 0   13.845

C - B == 0   12.236

> summary(tukey2)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = cal ~ sex + way, data = data)

Linear Hypotheses:

  Estimate Std. Error t value Pr(>|t|)   

B - A == 0   1.6090     0.7851   2.049    0.121   

C - A == 0  13.8450     0.7851  17.634   <0.001 ***

C - B == 0  12.2360     0.7851  15.584   <0.001 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(tukey2)

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 방법 C와 A, C와 B는 신뢰구간을 0을 포함하지 않으므로 유의한 차이가 있다는 결론을 내린다 (p-value < 0.001)

출처: 보건정보데이터 분석(이태림 저), R을 이용한 누구나 하는 통계분석(안재형 저)

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

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

 목표 변수가 연속형 변수인 경우, 의사결정나무는 회귀나무라고 한다. 통계학에서 목표변수가 연속형일 때 흔히 수행하는 분석이 선형 회귀분석인 것을 연상해 본다면 명칭이 이해가 될 것이다. 일반적으로 범주형 변수들은 가변수화하여 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)

제3장 연속형 자료의 분석

3.1 두 집단의 평균비교

3.1.1 독립표본의 평균비교 two sample test

예제) 흡연자 집단과 비흡연자 집단 간 폐 파괴지수를 측정하였다. 높은 수치는 폐의 손상이 크다는 것을 뜻한다. 흡연자와 비흡연자의 폐 파괴지수의 평균이 같다고 할 수 있는가? (각 그룹에서의 관측치들은 정규분포를 따르는 모집단으로부터 독립적으로 얻어진 것이며 두 그룹에서의 모분산은 같다고 가정하자. )

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

#1. 자료 입력

> smoke=c(16.6,13.9,11.3,26.5,17.4,15.3,15.8,12.3,18.6,12,24.1,16.5,21.8,16.3,23.4,18.8)

> nonsmoke=c(18.1,6,10.8,11,7.7,17.9,8.5,13,18.9)

> sapply(list(smoke,nonsmoke),mean)

[1] 17.53750 12.43333

> sapply(list(smoke,nonsmoke),sd)

[1] 4.475247 4.849227

#2. 정규성 검정

> qqnorm(smoke,main='smoke')

> qqline(smoke,col='blue')

> shapiro.test(smoke)

Shapiro-Wilk normality test

data:  smoke

W = 0.94511, p-value = 0.4163

결과 해석: shapiro.test의 결과에 따라 p value = 0.4163 > 0.05 이므로 귀무가설 기각 못한다, 즉 정규분포를 따른다

> qqnorm(nonsmoke,main = 'nonsmoke')

> qqline(nonsmoke,col='red')

> shapiro.test(nonsmoke)

Shapiro-Wilk normality test

data:  nonsmoke

W = 0.90366, p-value = 0.274

#boxplot과 vioplot

> boxplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

> library(vioplot)

> vioplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

결과 해석:  두 집단에 차이가 있음을 알 수 있다.

#3. 두 모분산 비교 (양측검정)

#대립 가설의 형태: alternative = c('two.sided','less','greater')

> var.test(smoke,nonsmoke)

F test to compare two variances

data:  smoke and nonsmoke

F = 0.8517, num df = 15, denom df = 8, p-value = 0.7498

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

  0.2076714 2.7243799

sample estimates:

  ratio of variances

0.8517046

결과 해석: p value 가 0.7498로 분산이 같다는 귀무가설을 기각할 수 없다. 즉, 등분산 가정

#4.두 모분산 비교 (양측검정) - 등분산 가정

> t.test(smoke,nonsmoke,var.equal = T)

Two Sample t-test

data:  smoke and nonsmoke

t = 2.658, df = 23, p-value = 0.01405

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  1.131680 9.076653

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석:

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

결정: p value는 p value는 0.01405로 두 모 평균이 같다는 귀무가설을 기각한다. 즉, 두 모평균이 서로 다르다.

#5.두 모분산 비교 (양측검정) - 이분산 가정

> t.test(smoke,nonsmoke)

Welch Two Sample t-test

data:  smoke and nonsmoke

t = 2.5964, df = 15.593, p-value = 0.01978

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  0.9279143 9.2804190

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석: 등분산 가정과 큰 차이는 없다.

3.1.2 짝지은 표본의 평균비교 paired sample test

예) 환자 15명에게 혈압강하제를 12주 투입 후 혈압을 비교하였다. 새로운 약은 효과적인가?

귀무가설 h0: u1-u2 = 0

대립가설 h1: u1 > u2

1) 데이터 입력

> before=c(90,56,49,64,65,88,62,91,74,93,55,71,54,64,54)

> after=c(72,55,56,57,62,79,55,72,73,74,58,59,58,71,61)

> diff = before - after

2) 정규성 차이: shapiro – wilk test

> qqnorm(diff)

> qqline(diff,col='red')

> shapiro.test(diff)

Shapiro-Wilk normality test

data:  diff

W = 0.90982, p-value = 0.1345

결과 해석: shapiro test 결과 p value는 0.1345로써 정규분포를 이루고 있다고 할 수 있다.

3) Paired sample test

> mean(diff) ; sd(diff)

[1] 4.533333

[1] 9.425396

> t.test(before, after, paired = T, alternative = 'greater') #μ 복용 전 > μ 복용 후

Paired t-test

data:  before and after

t = 1.8628, df = 14, p-value = 0.0418

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

  0.2469617       Inf

sample estimates:

  mean of the differences

4.533333

귀무가설 H0: 약 복용전과 복용 후의 혈압 수치는 같다. μ 복용 전 = μ 복용 후

대립가설 H1: 약 복용전 대비 복용 후의 혈압 수치기 다 낮다. μ 복용 전 > μ 복용 후

결론: 단측 검정에 대한 p value가 0.0418로서 유의수준 5%에서 두 그룹의 혈압 차이가 없다는 귀무가설을 기각할 만한 충분한 증거가 있으므로 새로운 약이 혈압을 내린다고 볼 수 있다.

출처: 보건 정보 데이터 분석(이태림 저자)

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

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

예제) 목표변수가 범주형 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)