데이터마이너를 꿈꾸며

독일 신용평가 데이터를 활용한 신경망 모형. 목표변수 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 그래프가 무엇을 의미하는지 더 자세히 해석할 능력이 필요함.

-> 교수님 답변: 신경망 모형의 해석은 그림을 보고 해석하기가 상당히 힘듭니다. 워낙 복잡한 함수의 결합이기 때문입니다. 다만, 화살표 위의 수치(절대값)를 보고 연결강도가 강한지 여부를 판단할 수 있습니다. 신경망 모형의 태생적인 한계점인 것 같습니다.

참고로 교재에서 사용하였던 neuralnet 패키지는 plot을 제공합니다. 과거 강의에서 사용했던 다른 패키지는 직접적으로 plot을 산출할 수는 없었습니다.

 

3. 위 신경망 모형을 통해서 오분류율은 40.4%로 나왔는데 너무 높은건 아닌지 생각됨.

 -> 교수님 답변: 오분율은 상대적인 것이긴 하지만 높은 수준으로 보입니다. 교재에서는 예측 판별 기준으로 0.5라는 값을 사용했는데 실제 실무에서는 이 값을 적절하게 변경시키면서 오분류율을 낮추어주는 것이 좋을 것입니다. 예를 들어 원 데이터가 1이 0에 비해 많이 분포되어 있다면, 0.5보다 큰 값을 기준치로 삼는 편이 좋습니다.

 

4. 은닉층의 개수는 어떻게 설정하면 좋을지?

 

5. 가장 궁금한 것은 이 모델을 어떻게 실무에 적용할지 잘 모르겠음. 위 독일 신용평가 데이터로 신경망 모형을 만들고 오분류율도 체크하고, 실제값과 예측값도 비교하였는데, 이것이 의미하는 것들. 가령, 변수의 중요성, 그래서 신용도가 좋은 경우는 어느 경우이고, 또 다른 고객 데이터 셋이 있는데 이 새로운 셋에서 어떤 고객, 변수, 값들을 추출해야 우리가 원하는 우수한 고객을 알아낼 수 있는지, 실제 비즈니스 적용 포인트가 가장 궁금함.

 -> 교수님 답변: 실무에서 신경망 모형을 잘 해석하여 활용하기는 어려운 것이 사실입니다. 다만, 상대적으로 예측력이 높다는 점을 살려 새로운 데이터의 모든 변수들을 활용하여 1이나 0값을 예측하거나 목표변수의 값 자체를 산출하는 목적으로는 유용성이 높다고 봅니다. 변수의 해석이나 변수의 선택보다는 예측이나 분류 자체의 목적으로 사용하기에 적합하다고 보시면 좋습니다.

 

6. 전반적으로 의사결정나무, 앙상블 모형에 이어 데이터 마이닝 분야에 깊은 관심을 가지게 되는 좋은 계기.

 

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

보스턴 하우징을 이용한 신경망 모형. 목표변수를 주택가격의 중간값인 medv(연속형 변수)로 하고, 나머지 변수를 입력변수로 하는 신경망 모형

1) 데이터 입력

> set.seed(100)

> library(MASS)

> library(neuralnet)

> bdat = Boston

> str(bdat)

'data.frame':     506 obs. of  15 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 ...

$ medv.hat: num  23.7 23.7 32.4 32.4 32.4 ...

 

주의: chas, rad를 수치형 변수로 바꾸고, 앙상블 모형에서 사용하였던 medv.hat 변수를 삭제해야 한다.

신경망에서는 범주형 데이터를 수치화하여 적용한다. 대체로 순서가 의미 있는 범부형 데이터는 수치 전환을 한 뒤 표준화하여 사용하게 된다.

 

2) 데이터 및 타입 변경

> bdat<-bdat[,-15]

> bdat$chas = as.numeric(bdat$chas)

> bdat$rad = as.numeric(bdat$rad)

> class(bdat$chas);class(bdat$rad)

[1] "numeric"

[1] "numeric"

 

3) 50% 랜덤 추출

> i = sample(1:nrow(bdat), round(0.5*nrow(bdat)))

> max1 = apply(bdat, 2, max)

> min1 = apply(bdat, 2, min)

결과 해석 : 50% 랜덤 추출하여 훈련데이터(train)로 저장하고 나머지는 검증 데이터(test)로 저장

 

 

4) 변수의 표준화 과정

> sdat = scale(bdat,center=min1,scale = max1 - min1) 

> sdat = as.data.frame(sdat)

> head(sdat,3)

crim          zn               indus  chas                nox                   rm                 age                   dis   rad

1 0.0000000000000 0.18 0.06781524927    0 0.3148148148 0.5775052692 0.6416065911 0.2692031391 0.000

2 0.0002359225392 0.00 0.24230205279    0 0.1728395062 0.5479977007 0.7826982492 0.3489619802 0.125

3 0.0002356977440 0.00 0.24230205279    0 0.1728395062 0.6943858977 0.5993820803 0.3489619802 0.125

tax               ptratio              black                     lstat         medv

1 0.2080152672 0.2872340426 1.0000000000 0.08967991170 0.4222222222

2 0.1049618321 0.5531914894 1.0000000000 0.20447019868 0.3688888889

3 0.1049618321 0.5531914894 0.9897372535 0.06346578366 0.6600000000

 

함수 설명:

> sdat = scale(bdat,center=min1,scale = max1 - min1) 신경망에서 사용하는 수치형 변수를 0과 1사이의 값으로 바꾸어 주기 위한 단계. 위의 [0-1] 변환은 (x-min(x) /(max(x)-min(x)) 의 수식으로도 계산 가능. 

 

 

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

> train = sdat[i,] #학습, 훈련샘플 training

> test = sdat[-i,] #테스트샘플 test

> n = names(train) ;n

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

[11] "ptratio" "black"   "lstat"   "medv"  

> form = as.formula(paste('medv~',paste(n[!n %in% 'medv'],collapse = '+'))) 

> form

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat

> nn1 = neuralnet(form,data=train,hidden = c(5,3),linear.output = T)

> summary(nn1) #작성된 신경망모형의 오브젝트의 구조를 정리

                       Length    Class     Mode   

call                          5   -none-     call   

response              253   -none-     numeric

covariate            3289   -none-     numeric

model.list                 2   -none-     list   

err.fct                      1   -none-     function

act.fct                      1   -none-     function

linear.output             1   -none-     logical

data                       14   data.frame list   

net.result                 1   -none-     list   

weights                    1   -none-     list   

startweights              1   -none-     list   

generalized.weights   1   -none-     list   

result.matrix            95   -none-     numeric

함수 설명

> form = as.formula(paste('medv~',paste(n[!n %in% 'medv'],collapse = '+'))) #종속변수는 mdev, 나머지는 독립변수이다 틸다 ~. 와같은 개념. 종속변수를 제외한 n 안에 있는 모든 변수명을 다 집어넣고 더하여라. 변수명을 직접 입력해도 된다.

실제 form라고 입력하면 medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat

 

라는 공식을 볼 수 있다.

 

> nn1 = neuralnet(form,data=train,hidden = c(5,3),linear.output = T) #은닉층은 2개이고 처음은 5개 나머지는 3개, 회귀의 문제이므로 linear.output = T

 

 

질문: 은닉층의 수는 어떻게 결정할까? 통상 은닉층의 마디수는 입력층의 마디수 두배를 넘지 않도록 해야한다. 위에서는 c(5,3)으로 결정하였는데 이 수치가 최적의 수치인지는 어떻게 알까? 더 공부가 필요하다.

 

<참고> 딥러닝(deep learning)

기계학습 기법 중 하나로 신경망모형으로부터 비롯된 딥러닝(deep learning)은 기본적으로 은닉층이 많이 쌓여 가면서 복잡하고 깊은 구조로 발전하면서 deep 이라는 이름이 붙여졌다. 입력변수와 출력변수 간 복잡한 관계를 가중치를 통해 조정할 수 있는 구조와 매커니즘을 갖고 있다.

 

 

> plot(nn1)

 

 

①로 표시된 것이 상수항에 해당하고 가중치는 각각의 화살 표 위에 출력된다. 해석 공부도 더 필요하다.

처음 은닉층은 5개, 두번째 은닉층은 3개를 가지고 있는 신경망 모형이 완성된다. 

6) 모형 추정: 위의 그래프는 50%의 training data만을 사용하였으니 나머지 50% test data를 쓰자.

> pred.nn0 = compute(nn1,test[,1:13])

> summary(pred.nn0)
               Length Class  Mode  
    neurons      3    -none- list  
    net.result 253    -none- numeric

> names(pred.nn0)

[1] "neurons"    "net.result"

> pred0 = pred.nn0$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv)

> head(pred0)
                    [,1]    
    1 26.55282080
    2 25.15912380
    3 34.48852794
    4 31.56317828
    5 32.47594868
    6 25.32302689

 

함수 설명:

> pred.nn0 = compute(nn1,test[,1:13])

14번째 변수가 medv 목표변수.compute는 작성된 신경망모형을 이용하여 새로운 예에 적용하여 결과를 도출. nn1는 새롭게 예측에 적용할 자료, test[,1:13]는 신경망모형으로 적합한 결과 오브젝트

> pred0 = pred.nn0$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv)

목표변수를 조정전 변수값으로 환원하는 과정. 원래 가지고 있는 값으로 환원. 역변환 과정

 

 

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

아래 함수 head(cbind(bdat[i,14],pred0),50)에서 cbind(bdat[i,14],는 실제값(training data)과 예측값(neural network)으로 구성된 행렬을 구성하는 함수

> head(cbind(bdat[i,14],pred0),50)

[,1]        [,2]

1  15.6 28.12482177

2  19.2 21.36405177

3  29.1 34.24035993

4  18.4 35.88436410

6  24.0 25.77273199

9  22.2 14.86779840

12 11.9 20.82595913

14 26.4 20.03825755

18 24.4 17.19077183

19 23.9 18.69934412

21 20.3 14.85391889

23  9.6 15.77921679

24 13.3 14.95074253

25 33.3 16.19467498

26 15.0 15.53094329

27 19.3 16.50978025

28 27.5 15.65043905

30 22.6 19.71434020

33 29.4 13.39051466

34 19.5 15.45310394

36 33.8 22.75202201

41 31.2 37.96410157

42 21.2 36.32787047

43 19.9 27.29495234

44 42.3 26.74797523

46 24.8 20.55343697

47 50.0 19.83522653

48 20.8 16.66204238

50 48.8 16.90901569

51 23.0 19.54979162

54 41.7 22.43724314

55 17.1 16.69040814

57 25.0 22.91956778

58 15.2 29.22726384

62  8.1 18.82030494

63  8.8 24.13500281

66 19.7 25.15181243

67 28.6 20.16650282

68 20.2 20.27218139

69 18.7 17.48036500

70 17.0 19.87776814

71 12.1 25.21432307

73 25.0 23.16314867

75 12.3 23.82067236

77 23.9 20.42068536

79 37.6 20.09872166

80 22.7 20.89416546

81  5.0 28.30078783

84 28.4 22.10685814

87 14.0 21.30514814

결과 해석: 26번째 행의 값처럼 유사한 결과값도 있지만 81번 결과 값처럼 서로 차이가 나는 경우도 있다.

 

 

8) 예측 정확도 평가: 모형 추정에 사용되지 않은 test data를 가지고 예측을 해보자. 

> pred.nn1 = compute(nn1,test[,1:13]) #목표 변수 제외, 새로운 변수로 간주하고 nn1에 적용

> pred1 = pred.nn1$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv) #목표변수를 조정전 변수값으로 환원. 역변환

> head(cbind(bdat[-i,14],pred1),50) #좌측은 test 변수의 실제값, 우측은 예측값(적합값), 14번째 값은 목표변수, pred1은 예측값(적합값)

[,1]        [,2]

1  24.0 28.12482177

2  21.6 21.36405177

3  34.7 34.24035993

4  33.4 35.88436410

6  28.7 25.77273199

9  16.5 14.86779840

12 18.9 20.82595913

14 20.4 20.03825755

18 17.5 17.19077183

19 20.2 18.69934412

21 13.6 14.85391889

23 15.2 15.77921679

24 14.5 14.95074253

25 15.6 16.19467498

26 13.9 15.53094329

27 16.6 16.50978025

28 14.8 15.65043905

30 21.0 19.71434020

33 13.2 13.39051466

34 13.1 15.45310394

36 18.9 22.75202201

41 34.9 37.96410157

42 26.6 36.32787047

43 25.3 27.29495234

44 24.7 26.74797523

46 19.3 20.55343697

47 20.0 19.83522653

48 16.6 16.66204238

50 19.4 16.90901569

51 19.7 19.54979162

54 23.4 22.43724314

55 18.9 16.69040814

57 24.7 22.91956778

58 31.6 29.22726384

62 16.0 18.82030494

63 22.2 24.13500281

66 23.5 25.15181243

67 19.4 20.16650282

68 22.0 20.27218139

69 17.4 17.48036500

70 20.9 19.87776814

71 24.2 25.21432307

73 22.8 23.16314867

75 24.1 23.82067236

77 20.0 20.42068536

79 21.2 20.09872166

80 20.3 20.89416546

81 28.0 28.30078783

84 22.9 22.10685814

87 22.5 21.30514814

 

결과해석: 전반적으로 실제값과 예측값이 서로 유사함을 알 수 있다.

 

> head(cbind(bdat[-i,14],pred1)[,1],10)
   1    2    3    4    5    6    7    8    9   10
24.0 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9
> head(cbind(bdat[-i,14],pred1),10)
   [,1]        [,2]
1  24.0 26.55282080
2  21.6 25.15912380
3  34.7 34.48852794
4  33.4 31.56317828
5  36.2 32.47594868
6  28.7 25.32302689
7  22.9 19.92609580
8  27.1 20.01203029
9  16.5 19.68915018
10 18.9 19.33552651
> obs <-cbind(bdat[-i,14],pred1)[,1]
> pdt <-cbind(bdat[-i,14],pred1)[,2]

> out = cbind(obs,pdt)
> head(out)
   obs         pdt
1 24.0 26.55282080
2 21.6 25.15912380
3 34.7 34.48852794
4 33.4 31.56317828
5 36.2 32.47594868
6 28.7 25.32302689
> which.min(abs(out$obs - out$pdt))
Error in out$obs : $ operator is invalid for atomic vectors
> out = as.data.frame(out)
> which.min(abs(out$obs - out$pdt))
[1] 43
> out[43,]
    obs         pdt
43 25.3 25.31048382
> which.max(abs(out$obs - out$pdt))
[1] 197
> out[197,]
    obs         pdt
372  50 18.13253008

 

 

9) PMSE: 예측된 값이 얼마나 잘 예측되었을까? 

> PMSE = sum((bdat[-i,14] - pred1)^2) / nrow(test)

> PMSE

[1] 18.70948041

함수 설명:

> PMSE = sum((bdat[-i,14] - pred1)^2) / nrow(test) #예측 평균 제곱 오차, 

 

 결과 해석: 예측평균제곱오차는 18.7

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

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

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

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

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

변수명	속성	변수 설명
crim	수치형(numeric)	per capita crime rate by town 타운별 1인당 범죄율
zn	수치형(numeric)	proportion of residential land zoned for lots over 25,000 sq.ft. 25,000평방피트를 초과하는 거주지역 비율
indus	수치형(numeric)	proportion of non-retail business acres per town. 비소매 사업지역의 토지 비율
chas	범주형(integer)	Charles River dummy variable (= 1 if tract bounds river; 0 otherwise). 찰스강 더비 변수 (강의 경계에 위치 = 1, 아니면 = 0)
nox	수치형(numeric)	nitrogen oxides concentration (parts per 10 million). 10ppm당 농축 일산화질소
rm	수치형(numeric)	average number of rooms per dwelling. 주택 1가구등 방의 평균 개수
age	수치형(numeric)	proportion of owner-occupied units built prior to 1940. 1940년 이전에 건축된 소유자 주택 비율
dis	수치형(numeric)	weighted mean of distances to five Boston employment centres. 5개의 보스턴 고용센터까지의 접근성 지수
rad	범주형(integer)	index of accessibility to radial highways. 방사형 도로까지의 접근성 지수
tax	수치형(numeric)	full-value property-tax rate per \$10,000. 10,000달러당 재산세율
ptratio	수치형(numeric)	pupil-teacher ratio by town. 타운별 학생/교사 비율
black	수치형(numeric)	1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town. 타운별 흑인의 비율
lstat	수치형(numeric)	lower status of the population (percent). 모집단의 하위계층의 비율
medv (목표변수)	수치형(numeric)	median value of owner-occupied homes in \$1000s. 본인 소유의 주택가격(중앙값)

 

 

 

 

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로 아주 작은 것도 모형이 적합하다는 것을 지지하다.  

 

 

 

[참고] 단계적선택법(stepwise selection)의 AIC는 1291.93이다. 후진소거법과 전친선택법은??

후진소거법(backward elimination)의 AIC는 1291.93

> fit.step.back = step(fit1,direction='backward')

Start:  AIC=1295.03

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat

 

          Df Sum of Sq    RSS    AIC

- chas     1      0.63 6321.5 1293.1

- indus    1     11.30 6332.2 1293.9

<none>                 6320.9 1295.0

- age      1     82.67 6403.5 1299.4

- tax      1    126.28 6447.1 1302.7

- nox      1    189.88 6510.7 1307.5

- black    1    192.42 6513.3 1307.7

- zn       1    203.44 6524.3 1308.5

- crim     1    228.82 6549.7 1310.5

- rad      8    721.85 7042.7 1332.0

- ptratio  1    706.41 7027.3 1344.9

- dis      1    860.51 7181.4 1355.6

- lstat    1    965.26 7286.1 1362.7

- rm       1   1330.92 7651.8 1386.7

 

Step:  AIC=1293.08

medv ~ crim + zn + indus + nox + rm + age + dis + rad + tax +

    ptratio + black + lstat

 

          Df Sum of Sq    RSS    AIC

- indus    1     11.00 6332.5 1291.9

<none>                 6321.5 1293.1

- age      1     82.48 6404.0 1297.4

- tax      1    130.45 6451.9 1301.1

- nox      1    189.27 6510.8 1305.5

- black    1    193.59 6515.1 1305.9

- zn       1    203.76 6525.2 1306.6

- crim     1    230.58 6552.1 1308.6

- rad      8    738.26 7059.8 1331.2

- ptratio  1    719.40 7040.9 1343.9

- dis      1    861.64 7183.1 1353.7

- lstat    1    965.11 7286.6 1360.7

- rm       1   1333.37 7654.9 1384.9

 

Step:  AIC=1291.93

medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio +

    black + lstat

 

          Df Sum of Sq    RSS    AIC

<none>                 6332.5 1291.9

- age      1     81.09 6413.6 1296.2

- tax      1    192.78 6525.3 1304.6

- black    1    196.55 6529.0 1304.9

- zn       1    220.63 6553.1 1306.7

- crim     1    225.50 6558.0 1307.1

- nox      1    239.09 6571.6 1308.1

- rad      8    791.09 7123.6 1333.6

- ptratio  1    732.81 7065.3 1343.6

- dis      1    857.27 7189.8 1352.1

- lstat    1    987.73 7320.2 1361.0

- rm       1   1380.21 7712.7 1386.5

 

> summary(fit.step.back )

 

Call:

lm(formula = medv ~ crim + zn + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat, data = boston)

 

Residuals:

    Min      1Q  Median      3Q     Max

-9.5200 -2.2850 -0.4688  1.7535 15.3972

 

Coefficients:

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

(Intercept)  30.252522   4.329907   6.987 9.64e-12 ***

crim         -0.104568   0.025533  -4.095 4.96e-05 ***

zn            0.045510   0.011235   4.051 5.97e-05 ***

nox         -12.366882   2.932651  -4.217 2.97e-05 ***

rm            3.583130   0.353644  10.132  < 2e-16 ***

age          -0.025822   0.010514  -2.456 0.014412 * 

dis          -1.253903   0.157029  -7.985 1.08e-14 ***

rad2          2.387130   1.160735   2.057 0.040278 * 

rad3          4.644091   1.062157   4.372 1.51e-05 ***

rad4          2.608777   0.944668   2.762 0.005977 **

rad5          3.116933   0.960550   3.245 0.001258 **

rad6          1.422890   1.150280   1.237 0.216705   

rad7          4.868388   1.240114   3.926 9.94e-05 ***

rad8          5.872144   1.180865   4.973 9.26e-07 ***

rad9          6.420553   1.393304   4.608 5.24e-06 ***

tax          -0.010571   0.002792  -3.787 0.000172 ***

ptratio      -0.837356   0.113420  -7.383 7.08e-13 ***

black         0.007949   0.002079   3.823 0.000149 ***

lstat        -0.357576   0.041718  -8.571  < 2e-16 ***

---

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

 

Residual standard error: 3.667 on 471 degrees of freedom

Multiple R-squared:  0.7907,    Adjusted R-squared:  0.7827

F-statistic: 98.83 on 18 and 471 DF,  p-value: < 2.2e-16

 

 

[참고] 전진선택법(forward selection)의 AIC는 1295.03

> fit.step.forward = step(fit1,direction='forward')

Start:  AIC=1295.03

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

tax + ptratio + black + lstat

 

> summary(fit.step.forward)

 

Call:

lm(formula = medv ~ crim + zn + indus + chas + nox + rm + age +

    dis + rad + tax + ptratio + black + lstat, data = boston)

 

Residuals:

    Min      1Q  Median      3Q     Max

-9.5220 -2.2592 -0.4275  1.6778 15.2894

 

Coefficients:

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

(Intercept)  30.120918   4.338656   6.942 1.29e-11 ***

crim         -0.105648   0.025640  -4.120 4.47e-05 ***

zn            0.044104   0.011352   3.885 0.000117 ***

indus        -0.046743   0.051044  -0.916 0.360274   

chas2         0.158802   0.736742   0.216 0.829435   

nox         -11.576589   3.084187  -3.754 0.000196 ***

rm            3.543733   0.356605   9.937  < 2e-16 ***

age          -0.026082   0.010531  -2.477 0.013613 * 

dis          -1.282095   0.160452  -7.991 1.05e-14 ***

rad2          2.548109   1.175012   2.169 0.030616 * 

rad3          4.605849   1.064492   4.327 1.85e-05 ***

rad4          2.663393   0.950747   2.801 0.005299 **

rad5          3.077800   0.962725   3.197 0.001483 **

rad6          1.314892   1.157689   1.136 0.256624   

rad7          4.864208   1.241760   3.917 0.000103 ***

rad8          5.772296   1.194221   4.834 1.82e-06 ***

rad9          6.195415   1.417826   4.370 1.53e-05 ***

tax          -0.009396   0.003070  -3.061 0.002333 **

ptratio      -0.828498   0.114436  -7.240 1.85e-12 ***

black         0.007875   0.002084   3.779 0.000178 ***

lstat        -0.354606   0.041901  -8.463 3.36e-16 ***

---

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

 

Residual standard error: 3.671 on 469 degrees of freedom

Multiple R-squared:  0.7911,    Adjusted R-squared:  0.7821

F-statistic: 88.78 on 20 and 469 DF,  p-value: < 2.2e-16

 

 

 

 

 

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

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

 

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

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

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

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

 

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