데이터마이너를 꿈꾸며

비모수적 방법 non-parametric method

6.2.2 누적한계추정법 product-limit method

생명표 방법: 기간을 1년, 6개월 등 특정 단위로 나눠 구분

누적한계추정법: 매 사건이 발생할 때 마다 해당 시점을 표기

누적한계추정법(product-limit method)은 생존함수를 추정하는 대표적인 방법 중 하나로 연구자의 이름을 따서 Kaplan-Meier추정법이라고도 한다. 모든 환자의 생존시간 또는 중도절단 시간이 각각 관찰되었다고 하자. 모든 자료의 생존시간 도는 중도절단 x1,…,투을 순서대로 배열한 것을 t1<t2<…<tn이라 하고, δi = 0, 그렇지 않은 경우 δi =1로 정의한다. 즉 중도 절단 되지 않고 사건이 발생한 경우 status = 1, 중도절단된 경우 stats = 0 로 표기한다.

예] 신장이식수술을 받은 15명 환자들의 호전기간(remission duration)이 아래와 같다고 하자. 여기서 +로 표기된 것은 중도절단된 자료를 가르킨다. 각 시점에서 생존함수를 누적한계추정법을 이용하여 추정해보자

표. 신장이식 환자들의 호전기간 (단위 : 일)

1. 데이터 불러오기

> library(survival)

> setwd('c:/Rwork')

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

> kidney #status =1 사건, 0 = 절단, 절단 표시 매우 중요

time status

1   3.0      1

2   4.0      0

3   4.5      1

4   4.5      1

5   5.5      1

6   6.0      1

7   6.4      1

8   6.5      1

9   7.0      1

10  7.5      1

11  8.4      0

12 10.0      0

13 10.0      1

14 12.0      1

15 15.0      1

> attach(kidney)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

> summary(fit1)

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

time n.risk n.event survival std.err lower 95% CI upper 95% CI

3.0     15       1    0.933  0.0644       0.8153        1.000

4.5     13       2    0.790  0.1081       0.6039        1.000

5.5     11       1    0.718  0.1198       0.5177        0.996

6.0     10       1    0.646  0.1275       0.4389        0.951

6.4      9       1    0.574  0.1320       0.3660        0.901

6.5      8       1    0.503  0.1336       0.2984        0.846

7.0      7       1    0.431  0.1324       0.2358        0.787

7.5      6       1    0.359  0.1283       0.1781        0.723

10.0      4       1    0.269  0.1237       0.1094        0.663

12.0      2       1    0.135  0.1135       0.0258        0.703

15.0      1       1    0.000     NaN           NA           NA

> fit1

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

n  events  median 0.95LCL 0.95UCL

15      12       7       6      NA

Error: invalid multibyte character in parser at line 1

-> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

6.5일 에서의 생존확률은 50.3%, 7.0일 이 지나면 생존율이 43.1%로 50%이하가 된다.

3. 사망시점의 사분위수 추정치와 그의 신뢰구간, 가령 사망자가 50%되는 시점은 언제인가?

> quantile(fit1,probs=c(0.25,0.5,0.75),conf.int=T) #신뢰구간까지

$quantile

25   50   75

5.5  7.0 12.0

$lower

25  50  75

4.5 6.0 7.0

$upper

25  50  75

7.5  NA  NA

-> 생존확률이 75%인 경우, 즉 사망확률이 25%인 경우의 시점은 5.5일, 생존확률이 50%인 경우 7.0일, 생존확율이 25%인 경우의 시점은 12.0 일

4. 누적한계추정치의 95% 신뢰구간 그래프

> plot(fit1,xlab='time',ylab='Survival function',lwd=2)

> legend(0.2,0.3,c('KM estimate','95% CI'),lty=c(1,2))

-> 점선은 신뢰구간을 의미, 실선은 생존함수추정치.생존확율이 50%인 경우는 약 7.0일임을 알 수 있다.

출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)

로짓분석이란?

 생존여부나 교통사고 발성 여부와 같이 반응변수(Y)는 범주형(명목형(남자 1, 여자2), 순서형(초졸1, 중졸2, 고졸3, 대졸4))이고 이에 대한 설명변수(X)는 범주형과 이산형(방의 개수: 2개)dl 혼합된 경우 관련성 여부를 규명하기 위해서는 모형으로 로짓logit 모형을 적용할 수 있다.

•  오즈: 기본이 되는 변수로 오즈odds를 들 수 있는데 이것은 두 확률의 비를 의미하는 것이다. 

odds=p/(1-p) 

•   로짓: 흡연 산모가 미숙아를 출산할 확률이 비흡연 산모를 1로 했을 때 2.5배라고 발표하는 것이 그 예이다. 로짓은 오즈에 자연로그를 취한 형태를 의미한다. 

logit=ln(p/1-p) 

• 로짓 모형: 범주형 자료 분석에서 실험자들의 설명변수(X)에 대한 통제가 가능한 경우 반응변수(Y)에 대한 로짓모형을 적용할 수 있다.

예) B.J.T Morgan은 폴란드 바르샤바의 3918명 여성들을 대상으로 초사한 초경자료를 가지고 범주형 자료의 회귀분석(을 소개하였다). 연령과 월경 사이에 어떠한 관계가 있는지 알아보자.

> 초경연령자료 = read.table('c:/Rwork/바르샤바초경연령자료.csv',sep=',',header=T)

> head(초경연령자료)

  평균연령 초경경험자 군의총수

1     9.21          0      376

2    10.21          0      200

3    10.58          0       93

4    10.83          2      120

5    11.08          2       90

6    11.33          5       88

> attach(초경연령자료)

> 확률 = 초경경험자/군의총수

> plot(평균연령,확률)

#로짓분석

> logit <- glm(확률~평균연령,data=초경연령자료,family='binomial')

> logit

Call:  glm(formula = 확률 ~ 평균연령, family = "binomial", data = 초경연령자료)

Coefficients:

(Intercept)     평균연령 

    -20.907        1.608 

Degrees of Freedom: 23 Total (i.e. Null);  22 Residual

Null Deviance:      19.06

Residual Deviance: 0.2214       AIC: 11.38

> summary(logit)

Call:

glm(formula = 확률 ~ 평균연령, family = "binomial", data = 초경연령자료)

Deviance Residuals:

     Min        1Q    Median        3Q       Max 

-0.20043  -0.08458  -0.05277   0.06104   0.13296 

Coefficients:

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

(Intercept) -20.9067     8.1314  -2.571   0.0101 *

평균연령      1.6077     0.6244   2.575   0.0100 *

---

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

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 19.06046  on 23  degrees of freedom

Residual deviance:  0.22137  on 22  degrees of freedom

AIC: 11.383

Number of Fisher Scoring iterations: 6

결과: 로짓모형을 추정하면 아래와 같다. 즉, X가 1단위 증가할 때 변화하는 로그 오즈(odds)의 비율은 1.6077이다.

logit(p) = -20.9067 + 1.6077 X 

> exp(coef(logit))

 (Intercept)     평균연령

8.323842e-10 4.991105e+00

결과해석: 나이가 한 살 증가할 때 초경 경험자의 오즈 odds는 4.991 (exp(1.6077))배 더 많다고 할 수 있다.

그룹의 50%가 초경을 하는 연령을 구하기 위해 p=0.5일 때의 연령 X를 구하면 

이고 이것이 바로 유효 중앙값 ED₅₀(Effective Dose of 50%)이 되어 50%가 초경을 경험한 소녀들의 나이가 된다. 이를 계산하면 아래 수식을 통해 x=13, 즉 약 13세임을 알 수 있다.

0 = -20.9067 + 1.6077 X

이식은 양의 용량이 늘어날수록 반응률이 높아지는 반응곡선, 즉 용량-반응곡선(dose response curve)에도 응용된다. 예를 들면 100명의 고열 환자에게 해열제를 투여했을 때 환자 50%의 열이 떨어지는 효과를 보여주는 약의 용량인 ED₅₀을 구하는데 적용할 수 있다.

예) M.J.R Healy는 비타민 E의 용량에 따른 임신한 쥐의 숫자를 아래와 같이 발표하였다. 투입된 비타민 E의 용량과 임신에 대한 관계를 알아보자.

> 용량 = c(3.75,5,6.25,7.5,10,15)

> 수= c(5,10,10,10,11,11)

> 임신 = c(0,2,4,8,10,11)

> 쥐임신<-cbind(용량,수,임신)

> 쥐임신

      용량 수 임신

[1,]  3.75  5    0

[2,]  5.00 10    2

[3,]  6.25 10    4

[4,]  7.50 10    8

[5,] 10.00 11   10

[6,] 15.00 11   11

> 확률2 = 쥐임신[,3] / 쥐임신[,2]

> 확률2

[1] 0.0000000 0.2000000 0.4000000 0.8000000 0.9090909 1.0000000

#로짓분석

> logit2 <- glm(확률2~log10(쥐임신[,1]),family='binomial')

> logit2

Call:  glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")

Coefficients:

       (Intercept)  log10(쥐임신[, 1]) 

            -12.42               15.35 

Degrees of Freedom: 5 Total (i.e. Null);  4 Residual

Null Deviance:      4.297

Residual Deviance: 0.136        AIC: 6.315

> summary(logit2)

Call:

glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")

Deviance Residuals:

       1         2         3         4         5         6 

-0.23099   0.11888  -0.09852   0.15618  -0.16518   0.08464 

Coefficients:

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

(Intercept)          -12.42      10.39  -1.195    0.232

log10(쥐임신[, 1])    15.35      12.76   1.203    0.229

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4.29706  on 5  degrees of freedom

Residual deviance: 0.13603  on 4  degrees of freedom

AIC: 6.3152

Number of Fisher Scoring iterations: 6

결과: glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")를 보면 모형을 적합하는데 용량(X) 대신 log ₁₀(X)를 사용하였다. 위의 결과를 통해 다음 추정식을 회귀할 수 있다. 

따라서 비타민 E 용량이 10배 증가함에 따라 변화하는 임신 로그 오즈의 비율은15.35임을 알 수 있다.

> exp(coef(logit2))

       (Intercept) log10(쥐임신[, 1])

      4.033060e-06       4.631322e+06

비타민 용량 E가 10배 증가할 때 임신의 오즈는 4.631(=exp(15.35))배 더 많다고 할 수 있다. 

전체 그룹 중 50%가 임신하는 비타민 용량 E는 아래와 같이 구할 수 있다.

수식을 이용하면 ED₅₀(Effective Dose of 50%)은 0 = -12.42 + 15.35 log₁₀x이므로 ED₅₀(Effective Dose of 50%)의 추정량은 6.44mg이 된다.

> 10^(12.42/15.35)

[1] 6.443481

출처: 보건정보데이터분석(이태림, 이재원, 김주한, 장대흥 공저)

p.46 제2장 보건 정보 데이터의 기초 분석

2.3 자료의 기술 및 요약

> setwd('c:/Rwork/')

> 담즙과포화비율자료 = read.table('담즙과포화비율.txt',header=T)

> attach(담즙과포화비율자료)

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자               40

2 남자               88

> plot(담즙의과포화비율,type='p',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> par(new=T)

> plot(담즙의과포화비율,type='h',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> length(담즙의과포화비율)

[1] 60

> length(담즙의과포화비율) #길이

[1] 60

> sort(담즙의과포화비율,decreasing=T) #내림차순

[1] 146 142 137 128 127 123 123 120 118 116 112 111 110 110 107 106 106  98  91  90  89  88  88

[24]  88  87  87  86  86  86  84  84  82  80  80  80  79  78  77  76  76  75  74  73  73  69  67

[47]  66  66  65  65  58  58  57  56  55  52  52  47  40  35

> sum(담즙의과포화비율) #합

[1] 5185

> cumsum(담즙의과포화비율) #누적합

[1]   40  128  238  295  381  518  596  707  795  875  961 1041 1088 1194 1259 1333 1399 1478

[19] 1536 1659 1746 1834 1924 1980 2053 2165 2275 2393 2445 2551 2618 2683 2735 2819 2905 2940

[37] 3056 3132 3187 3260 3349 3476 3563 3705 3782 3858 3916 4007 4114 4212 4340 4424 4570 4645

[55] 4765 4845 4927 5050 5116 5185

> mean(담즙의과포화비율);median(담즙의과포화비율) #평균과 중앙값

[1] 86.41667

[1] 84

> mean(담즙의과포화비율,trim=1/10) #10% 절삭

[1] 85.4375

> var(담즙의과포화비율);sd(담즙의과포화비율) #분산과 표준편차

[1] 657.1624

[1] 25.63518

> fivenum(담즙의과포화비율) #다섯숫자요약

[1]  35.0  68.0  84.0 106.5 146.0

> quantile(담즙의과포화비율)

0%    25%    50%    75%   100%

35.00  68.50  84.00 106.25 146.00

> IQR(담즙의과포화비율) #3사분위수-1사분위수

[1] 37.75

> quantile(담즙의과포화비율)[4]-quantile(담즙의과포화비율)[2]

75%

37.75

> mad(담즙의과포화비율) #Median Absolut Deviation 각 데이터에서 중앙값을 뺀 후 절대값을 취한 값들의 중앙값

[1] 26.6868

> max(담즙의과포화비율);min(담즙의과포화비율)

[1] 146

[1] 35

> range(담즙의과포화비율)

[1]  35 146

> R=max(담즙의과포화비율)-min(담즙의과포화비율)

> R                          

[1] 111

2. 4 표와 그래프를 이용한 자료의 요약 

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자               40

2 남자               88

> 계급 = cut(담즙의과포화비율, breaks=c(20,40,60,80,100,120,140,160))

> head(계급)

[1] (20,40]   (80,100]  (100,120] (40,60]   (80,100]  (120,140]

Levels: (20,40] (40,60] (60,80] (80,100] (100,120] (120,140] (140,160]

> table(계급)

계급

(20,40]   (40,60]   (60,80]  (80,100] (100,120] (120,140] (140,160]

2         8        18        15        10         5         2

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

> hist(담즙의과포화비율,breaks=c(20,40,60,80,100,120,140,160),main='히스토그램')

> rug(담즙의과포화비율)

> stem(담즙의과포화비율) #나무 줄기 그림

The decimal point is 1 digit(s) to the right of the | 

2 | 5

4 | 072256788

6 | 556679334566789

8 | 000244666778889018

10 | 667001268

12 | 033787

14 | 26

> stem(담즙의과포화비율,scale=2) #줄기의 마디 2배로 늘리기

The decimal point is 1 digit(s) to the right of the | 

3 | 5

4 | 07

5 | 2256788

6 | 556679

7 | 334566789

8 | 000244666778889

9 | 018

10 | 667

11 | 001268

12 | 03378

13 | 7

14 | 26

> library(vioplot)

> vioplot(담즙의과포화비율,names='담즙의과포화비율',col='yellow')

> n=length(담즙의과포화비율)

> plot(sort(담즙의과포화비율),(1:n)/n,type='s',ylim=c(0,1),main='Ogive of bile supersaturation',ylab='ECDF',xlab='담즙과포화비율') #ogive 오자이브 그래프

> rug(담즙의과포화비율)

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저)

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

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

 1) 데이터 입력

> setwd('c:/Rwork')

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

> german$numcredits = factor(german$numcredits)

> german$residence = factor(german$residence)

> german$residpeople = factor(german$residpeople)

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

> library(randomForest)

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

> summary(rf.german)

Length Class  Mode    

call               7   -none- call    

type               1   -none- character

predicted       1000   factor numeric 

err.rate         300   -none- numeric 

confusion          6   -none- numeric 

votes           2000   matrix numeric 

oob.times       1000   -none- numeric 

classes            2   -none- character

importance        80   -none- numeric 

importanceSD      60   -none- numeric 

localImportance    0   -none- NULL    

proximity          0   -none- NULL    

ntree              1   -none- numeric 

mtry               1   -none- numeric 

forest            14   -none- list    

y               1000   factor numeric 

test               0   -none- NULL    

inbag              0   -none- NULL    

terms              3   terms  call

함수 설명:

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

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

> names(rf.german)

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

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

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

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

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

1    2    3    4    5    6    7    8    9   10

good  bad good  bad  bad good good good good  bad

Levels: bad good

> importance(rf.german)

bad        good MeanDecreaseAccuracy MeanDecreaseGini

check       13.14268724  9.63208341         15.1106884       44.838547

duration     3.33563217  8.10482760          8.6281030        41.273512

history      3.87863720  4.50449203            5.8685013        25.731461

purpose      2.67025503  3.09871408            4.1078354        35.651165

credit       2.44312087  4.16670498            4.8800356        52.641745

savings      7.48182326  2.84190645            6.1879918        21.705909

employment   1.91049595  2.70568977            3.2416484        23.910509

installment  0.02100147  3.49542966            3.0522029        15.714499

personal     2.10802207  1.83819513            3.0602966        15.365475

debtors     -0.17277289  4.39384503            3.8963219         6.718969

residence    0.74571096  0.90284661            1.1155901        17.932937

property     2.16016716  3.85454658            4.7010080        18.803008

age          2.69637542  4.35170316            5.4753585        41.451103

others       0.31569112  3.60499256            3.3679530        11.399935

housing      2.70314243  2.06074416            3.0737900         9.596539

numcredits  -0.24996827  0.95259106            0.6502566         7.944861

job          1.53070574  1.18486660            2.1420488        12.054190

residpeople  0.88657814 -0.43166449            0.1370845         4.234831

telephone    1.46824003 -0.24200291            0.6110284         6.427048

foreign      1.26297478 -0.05431566            0.5733125         1.519832

> importance(rf.german,type=1)

MeanDecreaseAccuracy

check                 15.1106884

duration               8.6281030

history                5.8685013

purpose                4.1078354

credit                 4.8800356

savings                6.1879918

employment             3.2416484

installment            3.0522029

personal               3.0602966

debtors                3.8963219

residence              1.1155901

property               4.7010080

age                    5.4753585

others                 3.3679530

housing                3.0737900

numcredits             0.6502566

job                    2.1420488

residpeople            0.1370845

telephone              0.6110284

foreign                0.5733125

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

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

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

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

check

1

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

bad                 good     MeanDecreaseAccuracy     MeanDecreaseGini

13.142687             9.632083                15.110688            44.838547

> rf.german$importance

bad          good MeanDecreaseAccuracy MeanDecreaseGini

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

추가:

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

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

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

> head(pred.rf.german,10)

  1    2    3    4    5    6    7    8    9   10

  good  bad good good  bad good good good good  bad

  Levels: bad good

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

> tab

        Predicted

  Actual  bad   good

    bad    300    0

 good   0      700

> addmargins(tab)

  Predicted

  Actual  bad  good  Sum

  bad   300   0      300

  good  0      700  700

  Sum   300  700  1000

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

  [1] 0

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

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

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

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

> plot(rf.german,'l')

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

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

> set.seed(1234)

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

> german.train = german[i,]

> german.test = german[-i,]

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

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

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

> tab.train

Predicted

Actual bad good

bad   39   60

good  15  186

> addmargins(tab.train)

Predicted

Actual   bad  good Sum

bad    39   60     99

good  15   186    201

Sum   54   246    300

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

[1] 0.25

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

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

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

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

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

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

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

1. 데이터 불러오기

2. 선형 회귀 모형 만들기

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

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

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

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

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

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

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

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

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

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

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