데이터마이너를 꿈꾸며

독일 신용평가 데이터를 활용한 신경망 모형. 목표변수 y는 good / bad의 범주형 데이터로 모든 변수를 수치화 한 후 신경망 모형을 

1) 데이터 입력

> set.seed(1000)

> library(neuralnet)

> library(dummy)

> setwd('c:/Rwork')

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

추가 공부: dummy화란?

2) 데이터 및 타입 변경

> dvar=c(4,9,10,15,17) #명목변수 지정 purpose(a43,a40..), personal,  debtors, housing, job

> german2 = dummy(x=german[,dvar]) #명목변수를 더미변수화

> head(german2,1)

purpose_A40 purpose_A41 purpose_A410 purpose_A42 purpose_A43 purpose_A44 purpose_A45

1           0           0            0           0           1           0           0

purpose_A46 purpose_A48 purpose_A49 personal_A91 personal_A92 personal_A93 personal_A94

1           0           0           0            0            0            1            0

debtors_A101 debtors_A102 debtors_A103 housing_A151 housing_A152 housing_A153 job_A171

1            1            0            0            0            1            0        0

job_A172 job_A173 job_A174

1        0        1        0

> german2 = german2[,-c(10,14,17,20,24)] #더미변수생성

> head(german,1)

check  duration  history  purpose  credit  savings  employment  installment  personal  debtors

1   A11        6     A34      A43    1169      A65        A75           4      A93    A101

Residence  property  age  others  housing  numcredits  job  residpeople  telephone  foreign    y

1        4     A121  67   A143    A152          2 A173           1      A192    A201   good

> german2 = cbind(german[,-dvar],german2) #변수 결함

> str(german2)
'data.frame':   1000 obs. of  40 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 ...
 $ 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 ...
 $ 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 ...
 $ numcredits  : int  2 1 1 1 2 1 1 1 1 2 ...
 $ 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 ...
 $ purpose_A40 : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 1 2 ...
 $ purpose_A41 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
 $ purpose_A410: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ purpose_A42 : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 2 1 1 1 ...
 $ purpose_A43 : Factor w/ 2 levels "0","1": 2 2 1 1 1 1 1 1 2 1 ...
 $ purpose_A44 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ purpose_A45 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ purpose_A46 : Factor w/ 2 levels "0","1": 1 1 2 1 1 2 1 1 1 1 ...
 $ purpose_A48 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ purpose_A49 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ personal_A91: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 2 1 ...
 $ personal_A92: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 1 1 ...
 $ personal_A93: Factor w/ 2 levels "0","1": 2 1 2 2 2 2 2 2 1 1 ...
 $ personal_A94: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 2 ...
 $ debtors_A101: Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 2 2 2 ...
 $ debtors_A102: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ debtors_A103: Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 1 1 ...
 $ housing_A151: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
 $ housing_A152: Factor w/ 2 levels "0","1": 2 2 2 1 1 1 2 1 2 2 ...
 $ housing_A153: Factor w/ 2 levels "0","1": 1 1 1 2 2 2 1 1 1 1 ...
 $ job_A171    : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ job_A172    : Factor w/ 2 levels "0","1": 1 1 2 1 1 2 1 1 2 1 ...
 $ job_A173    : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 2 1 1 1 ...
 $ job_A174    : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 2 ...

> nrow(german2);ncol(german2)
[1] 1000
[1] 40

> for(i in 1:ncol(german2)) if(!is.numeric(german2[,i])) german2[,i] = as.numeric(german2[,i])#여타 순서가 있는 범주형 변수의 수치형 변수화

> german2$y = ifelse(german$y == 'good',1,0) #목표변수 변환

> head(german$y)

[1] good bad  good good bad  good

Levels: bad good

> head(german2$y)

[1] 1 0 1 1 0 1

 ## 중요 : 신경망에서는 범주형 데이터를 수치화하여 적용한다. ##

3) 75% 랜덤 추출

> i = sample(1:nrow(german2),round(0.75*nrow(german2))) #75%랜덤 추출

> length(i)
[1] 750

4) 변수의 표준화 과정

> max2 = apply(german2, 2, max)

> min2 = apply(german2, 2, min)

> gdat = scale(german2, center = min2, scale = max2 - min2) # 변수조정(0,1,dummy variable은 변화 없음)

> gdat = as.data.frame(gdat) #데이터 프레임 형태로 변경, 데이터 준비 끝!

> str(gdat)

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

$ check       : num  0 0.333 1 0 0 ...

$ duration    : num  0.0294 0.6471 0.1176 0.5588 0.2941 ...

$ history     : num  1 0.5 1 0.5 0.75 0.5 0.5 0.5 0.5 1 ...

$ credit      : num  0.0506 0.3137 0.1016 0.4199 0.2542 ...

$ savings     : num  1 0 0 0 0 1 0.5 0 0.75 0 ...

$ employment  : num  1 0.5 0.75 0.75 0.5 0.5 1 0.5 0.75 0 ...

$ installment : num  1 0.333 0.333 0.333 0.667 ...

$ residence   : num  1 0.333 0.667 1 1 ...

$ property    : num  0 0 0 0.333 1 ...

$ age         : num  0.8571 0.0536 0.5357 0.4643 0.6071 ...

$ others      : num  1 1 1 1 1 1 1 1 1 1 ...

$ numcredits  : num  0.333 0 0 0 0.333 ...

$ residpeople : num  0 0 1 1 1 1 0 0 0 0 ...

$ telephone   : num  1 0 0 0 0 1 0 1 0 0 ...

$ foreign     : num  0 0 0 0 0 0 0 0 0 0 ...

$ y           : num  1 0 1 1 0 1 1 1 1 0 ...

$ purpose_A40 : num  0 0 0 0 1 0 0 0 0 1 ...

$ purpose_A41 : num  0 0 0 0 0 0 0 1 0 0 ...

$ purpose_A410: num  0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A42 : num  0 0 0 1 0 0 1 0 0 0 ...

$ purpose_A43 : num  1 1 0 0 0 0 0 0 1 0 ...

$ purpose_A44 : num  0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A45 : num  0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A46 : num  0 0 1 0 0 1 0 0 0 0 ...

$ purpose_A48 : num  0 0 0 0 0 0 0 0 0 0 ...

$ personal_A91: num  0 0 0 0 0 0 0 0 1 0 ...

$ personal_A92: num  0 1 0 0 0 0 0 0 0 0 ...

$ personal_A93: num  1 0 1 1 1 1 1 1 0 0 ...

$ debtors_A101: num  1 1 1 0 1 1 1 1 1 1 ...

$ debtors_A102: num  0 0 0 0 0 0 0 0 0 0 ...

$ housing_A151: num  0 0 0 0 0 0 0 1 0 0 ...

$ housing_A152: num  1 1 1 0 0 0 1 0 1 1 ...

$ job_A171    : num  0 0 0 0 0 0 0 0 0 0 ...

$ job_A172    : num  0 0 1 0 0 1 0 0 1 0 ...

$ job_A173    : num  1 1 0 1 1 0 1 0 0 0 ...

5) 신경망 모델 구축 및 신경망 그래프 그리기

> train = gdat[i,] #학습샘플과 테스트 샘플 추출

> test = gdat[-i,]

> gn = names(german2)

> gn

[1] "check"        "duration"     "history"      "credit"       "savings"      "employment" 

[7] "installment"  "residence"    "property"     "age"          "others"       "numcredits" 

[13] "residpeople"  "telephone"    "foreign"      "y"            "purpose_A40"  "purpose_A41"

[19] "purpose_A410" "purpose_A42"  "purpose_A43"  "purpose_A44"  "purpose_A45"  "purpose_A46"

[25] "purpose_A48"  "personal_A91" "personal_A92" "personal_A93" "debtors_A101" "debtors_A102"

[31] "housing_A151" "housing_A152" "job_A171"     "job_A172"     "job_A173"   

> f = as.formula(paste('y~',paste(gn[!gn %in% 'y'],collapse = '+')))

> f

y ~ check + duration + history + credit + savings + employment +

 installment + residence + property + age + others + numcredits +

   residpeople + telephone + foreign + purpose_A40 + purpose_A41 +

   purpose_A410 + purpose_A42 + purpose_A43 + purpose_A44 +

 purpose_A45 + purpose_A46 + purpose_A48 + personal_A91 +

   personal_A92 + personal_A93 + debtors_A101 + debtors_A102 +

   housing_A151 + housing_A152 + job_A171 + job_A172 + job_A173

> nn1 = neuralnet(f,data=train,hidden=c(3,2),linear.output=F) #은닉층은 2개, 첫번째 노드는 3개, 두번째 노드는 2개. 분류의 경우 linear.output = F

> summary(nn1)

Length Class      Mode   

call                   5  -none-     call   

response            750  -none-     numeric

covariate          25500  -none-     numeric

model.list              2  -none-     list   

err.fct                 1  -none-     function

act.fct                 1  -none-     function

linear.output           1  -none-     logical

data                  35  data.frame list   

net.result              1  -none-     list   

weights                1  -none-     list   

startweights            1  -none-     list   

generalized.weights     1  -none-     list   

result.matrix         119  -none-     numeric

> plot(nn1)

은닉층이 2개인 신경망 모형의 그래프가 완성된다. 이 그래프에 대한 해석을 좀 더 공부하고 싶은데 아직은 잘 모르겠음.

6) 모형 추정:

> dim(german2)[1]
[1] 1000

> dim(german2)[2]
[1] 35

> colnames(test)[16]
[1] "y"

> pred.nn0 = compute(nn1,train[,c(1:15,17:dim(german2)[2])]) #학습데이터의 실제값과 예측값 비교, 16번째 열의 값은 y

함수 설명:

> pred.nn0 = compute(nn1,train[,c(1:15,17:dim(german2)[2])]) 16번째 변수가 y 목표변수. compute는 작성된 신경망모형을 이용하여 새로운 예에 적용하여 결과를 도출. nn1는 새롭게 예측에 적용할 자료, train[,c(1:15,17:dim(german2)[2])]는 신경망모형으로적합한 결과 오브젝트

7) 학습샘플의 실제값과 예측값을 비교해보자.

> head(cbind(german2[1,16],round(pred.nn0$net.result,10)))

[,1]         [,2]

931    1 0.7786470581

546    1 0.0000005387

56     1 0.8208161458

883    1 0.9999999722

11     1 0.0000004232

354    1 0.0000046419

#왼쪽이 실제값, 오른쪽이 예측값. 분류의 문제이므로 값은 0과 1사이. 0.5를 cut off값을 둘 수 있다. 또는 0.3미만 폐기, 0.7이하 보류, 0.7 초과만 사용한느 cutoff도 가능. 왼쪽이 실제 값, 오른쪽이 학습된 데이터. 4번째 행은 실제 1의 값을 1에 가깝게 예측하였고, 5번째 행은 실제 1이지만 0에 가깝게 예측한 사례. german2[,16]은 16번째 열 즉 y값임을 알겠는데 german2[1,16]은 뭘까… 궁금

8) 예측 정확도 평가

> pred.nn1 = compute(nn1,test[,c(1:15,17:dim(german2)[2])]) #test data를 바탕으로 판단해보자

> pred.nn2 = ifelse(pred.nn1$net.result>0.5,1,0) #0.5를 경계로 1과 0 분류

> head(cbind(german2[-i,16],pred.nn2)) #테스트 샘플의 실제값과 예측값

[,1]  [,2]

3     1    1

12    1    0

13    1    1

14    1    0

15    0    1

26    1    1

> sum(german2[-i,16]!=pred.nn2) / length(german2[-i,16])

[1] 0.404

> sum(german2[-i,16]!=pred.nn2)

[1] 101

> length(german2[-i,16])

[1] 250

#테스트 샘플의 오분류율, 16번째 값은 목표변수, sum(german2[-i,16]!=pred.nn2)는 pred.nn2와 같지 않은  값을 전체길이 length(german2[-i,16])로 나눔. i를 빼고 16번째 컬럼을 사용

> library(nnet)

> nnet1 =nnet(f,data=train,size = 3, linout = F)

# weights:  109

initial  value 181.570914

iter  10 value 113.679457

iter  20 value 96.943318

iter  30 value 82.803121

iter  40 value 73.239858

iter  50 value 70.807278

iter  60 value 69.865795

iter  70 value 69.476434

iter  80 value 69.158350

iter  90 value 69.039026

iter 100 value 68.929899

final  value 68.929899

stopped after 100 iterations

> pred.nnet1 = predict(nnet1,test[,c(1:15,17:dim(german2)[2])])

> pred.nnet2 = ifelse(pred.nnet1>0.5,1,0)

> head(cbind(german2[-i,16],pred.nnet2)) #테스트 샘플의 실제값과 예측값

[,1] [,2]

3     1    1

12    1    0

13    1    1

14    1    1

15    0    0

26    1    1

> sum(german2[-i,16]!=pred.nnet2) / length(german2[-i,16]) #테스트 샘플 예측의 오분류율

[1] 0.408

이 부분은 교재에 별도 설명이 없어서 추가 공부가 필요함.

소감: 알파고 딥마이닝으로 인해 관심을 가지게 된 신경망 모형. 이론 공부도 해보고 R도 따라해보니 약 50%정도 이해된 상태. 궁금한 점은 아래와 같음.

1. 명복 변수 중 더미화 하지 않은 것들도 있음.

-> 교수님 답변: 해당 신경망 모형에서는 숫자로 입력되어 있는 범주형 변수들은 수치변수로 그대로 사용하고 표준화만 하였습니다. 순서가 있는 범주형 변수라고 판단한 변수였습니다. 해당 변수들을 제외하고 나머지 변수들은 변환이 필요하여 dummy 함수를 사용하였습니다.

2. 위 신경망 plot 그래프가 무엇을 의미하는지 더 자세히 해석할 능력이 필요함.

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

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

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

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 분류나무 모형 바로가기)

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

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

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

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

독일신용평가 데이터 셋

독일신용평가 데이터(German Credit Data)는 머신러닝 저장소에 탑재되어 있는 데이터로 분류의 예제로 많이 활용된다.

1. 데이터 불러오기

2. 로지스틱 회귀 분석 시작

->  <!--[endif]-->해석: fit.step = step(fit.all, direction='both'를 통해 AIC가 가장 작은 모형을 찾는다.

check는 4개의 범주(checkA11 계좌 없음 / A12 잔액 없음 / A13 잔액 200 이하 / A14 잔액 200 이상)를 가지므로 3개의 가변 수 생성. 추정된 회귀계수는 모두 양수이므로, A12~A14 즉 계좌가 있는 경우 계좌 없음(A11)대비 신용이 좋을 확률(Y=1)이 더 높다. 대출기간인 duration은 마이너스의 값을 지니므로 대출 기간이 오래 될 수록 신용도는 낮아진다. 모델의 AIC는 980.91로, AIC가 클 경우 그 모형은 적합하지 않기 때문에, 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택한다.

3. 예측함수 및 정오분류표 작성

->  해석: 1로 예측할 확률이 임계치(threshold) 0.5보다 클 경우에는 1로, 0.5이하일 경우에는 0으로 예측. 실제로는 0인데 0으로 예측한 경우가 158개, 1인데 1로 분류한 경우가 618개이다.반면에 0인데 1로 오분류한 경우가 142개, 1인데 0으로 오분류한 경우가 82개이다.

4. 예측력 측도

-> 해석: 민감도는 실제 양성(Y=1)일 때 양성으로 예측할 확률, 특이도는 실제 음성(Y=0)일 때 음성으로 예측할 확률이다. 예측정확도(prediction accuracy)는 실제 양서일 때 양성으로, 음성일 때 음성으로 제대로 예측할 확률로 민감도와 특이도의 가중평균이다. 오분류율(misclassification rate)는 양성일 때 음성으로, 음성일 때 양성으로 잘못 예측할 확률이다.

5. ROC 곡선 및 AUC 생성

->  민감도와 특이도는 임계치에 다라 달라지고 임계치는 상황에 따라 다르게 결정할 수 이다. 여러 가능한 임계치에 대해 ‘1-특이도(Specificity)’를 가로축에, 민감도를 세로축에 놓고 그린 그래프를 ROC(Receiver operating characteristic) 곡선이라 한다. 민감도와 특이도가 높을수록 예측력이 좋다고 할 수 있기 때문에 ROC 곡선이 좌상단에 가까울수록 ROC 곡선 아래 면적인 AUC(area under the ROC curve)가 커지고, 예측력이 좋다고 할 수 있다.이 독일신용평가 데이터에 적합한 로지스틱 회귀모형에 대한 예측력의 측도인 AUC는 최대치 1보다 다소 작은 0.831로 상당히 높음을 알 수 있다.