제6.2장 비모수적 방법 - 1. 생명표 방법

KNOU/2 보건 정보 데이터 분석 2017. 1. 26. 09:36
반응형

비모수적 방법 non-parametric method

6.2.1생명표 방법 life table method

 

생명표 방법 개념 추가 설명

 

예] 협심증(angina pectoris)이 있는 2,418명의 남성들에 대한 생존자료이다. 이 생명표에서 경과기관에 따른 생존확률을 구해보자

 

진단 후 경과 기관(단위: 년)

사망자 수

중도절단 수

( 0 – 1 ]

456

0

( 1 – 2 ]

226

39

( 2 – 3 ]

152

22

( 3 – 4 ]

171

23

( 4 – 5 ]

135

24

( 5 – 6 ]

125

107

( 6 – 7 ]

83

133

( 7 – 8 ]

74

102

( 8 – 9 ]

51

68

( 9 – 10 ]

42

64

( 10 – 11 ]

43

45

( 11 – 12 ]

34

53

( 12 – 13 ]

18

33

( 13 – 14 ]

9

27

( 14 – 15 ]

6

23

( 15 –

0

30

 

 

1. 데이터 입력

> library(KMsurv)

> setwd('C:/Rwork')

> 협심증환자자료 = read.csv('협심증환자자료.csv')

> head(협심증환자자료)

  time censor freq

1  0.5      1  456

2  0.5      0    0

3  1.5      1  226

4  1.5      0   39

5  2.5      1  152

6  2.5      0   22

> attach(협심증환자자료)

 

2. lifetab 함수 연구

> lifetab

function (tis, ninit, nlost, nevent)

(이하 생략)

-> tis: 시간, 길이는 17, ninit = 전체 데이터, 2418명, nlost = 중도절단, nevent 사건발생 

 

3. nevent 자료 만들기

> which(censor==1) #사망 데이터

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31

> 집단.사망=협심증환자자료 [which(censor==1),]

> head(집단.사망)

   time censor freq

1   0.5      1  456

3   1.5      1  226

5   2.5      1  152

7   3.5      1  171

9   4.5      1  135

11  5.5      1  125

 

 4. nlost 자료 만들기

> which(censor==0) #절단 데이터

 [1]  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32

> 집단.절단=협심증환자자료[which(censor==0),]

> head(집단.절단)

   time censor freq

2   0.5      0    0

4   1.5      0   39

6   2.5      0   22

8   3.5      0   23

10  4.5      0   24

12  5.5      0  107

 

5. ninit전체데이터 자료 만들기

> 사망자수=집단.사망[,3]

> 사망자수

 [1] 456 226 152 171 135 125  83  74  51  42  43  34  18   9   6   0

> 절단자수=집단.절단[,3]

> 절단자수

 [1]   0  39  22  23  24 107 133 102  68  64  45  53  33  27  23  30

> 합=사망자수+절단자수

> 합

 [1] 456 265 174 194 159 232 216 176 119 106  88  87  51  36  29  30

> sum(합) #2418

[1] 2418

 

6. tis 시간 변수 자료 만들기

> 년 = floor(집단.사망$time) #floor 내림값

> 년

 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

> lt=length(년)

> lt

[1] 16

> length(년)

[1] 17

> 년[lt+1] = NA

> 년[lt+1]

[1] NA

 

7. 생명표 만들기 

> 생명표=lifetab(년,sum(합),절단자수,사망자수)

> 생명표

      nsubs nlost  nrisk nevent      surv        pdf     hazard     se.surv

0-1    2418     0 2418.0    456 1.0000000 0.18858561 0.20821918 0.000000000

1-2    1962    39 1942.5    226 0.8114144 0.09440394 0.12353102 0.007955134

2-3    1697    22 1686.0    152 0.7170105 0.06464151 0.09440994 0.009179397

3-4    1523    23 1511.5    171 0.6523689 0.07380423 0.11991585 0.009734736

4-5    1329    24 1317.0    135 0.5785647 0.05930618 0.10804322 0.010138361

5-6    1170   107 1116.5    125 0.5192585 0.05813463 0.11859583 0.010304216

6-7     938   133  871.5     83 0.4611239 0.04391656 0.10000000 0.010379949

7-8     722   102  671.0     74 0.4172073 0.04601094 0.11671924 0.010450930

8-9     546    68  512.0     51 0.3711964 0.03697464 0.10483042 0.010578887

9-10    427    64  395.0     42 0.3342218 0.03553750 0.11229947 0.010717477

10-11   321    45  298.5     43 0.2986843 0.04302654 0.15523466 0.010890741

11-12   233    53  206.5     34 0.2556577 0.04209376 0.17941953 0.011124244

12-13   146    33  129.5     18 0.2135639 0.02968456 0.14937759 0.011396799

13-14    95    27   81.5      9 0.1838794 0.02030570 0.11688312 0.011765989

14-15    59    23   47.5      6 0.1635737 0.02066194 0.13483146 0.012259921

15-NA    30    30   15.0      0 0.1429117         NA         NA 0.013300258

           se.pdf   se.hazard

0-1   0.007955134 0.009697769

1-2   0.005975178 0.008201472

2-3   0.005069200 0.007649121

3-4   0.005428013 0.009153696

4-5   0.004945997 0.009285301

5-6   0.005033980 0.010588867

6-7   0.004690538 0.010962697

7-8   0.005175094 0.013545211

8-9   0.005024599 0.014659017

9-10  0.005307615 0.017300846

10-11 0.006269963 0.023601647

11-12 0.006847514 0.030646128

12-13 0.006682743 0.035110295

13-14 0.006514794 0.038894448

14-15 0.008035120 0.054919485

15-NA          NA          NA


nsubs  nlost   nrisk  nevent       surv        
생존자수 중도절단수 유효인원수 사망자수 생존함수
pdf     hazard      se.surv se.pdf se.hazard
확률밀도 함수 위험함수 생존함수의 표준오차 확률밀도 함수의 표준오차 위험함수의 표준오차


-> 5번째 열 surv는 생존 함수, 7번째 열 hazard는 위험 함수. 협심증 환자의 약 19%가 1년 이내, 약 28%가 2년 이내에 사망하는 것으로 추정된다. 따라서 협심증 환자가 2년 이상 생존할 확률은 약 72%가 된다는 것을 알 수 있다. 또한 5년 이상 생존할 확률은 약 52%이다.

 

 

#전체 환자 중 생존확율이 70% 되는 지점

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> which.min(abs(생명표$surv-0.7))

[1] 3

> 생명표[3,]

    nsubs nlost nrisk nevent      surv        pdf     hazard     se.surv

2-3  1697    22  1686    152 0.7170105 0.06464151 0.09440994 0.009179397

       se.pdf   se.hazard

2-3 0.0050692 0.007649121

-> abs:절대값, which.min 최소값, 또는 생명표에서 3번째 행을 보면 survival이 70%에 가장 가까이 접근했음을 알 수 있다. 2년 이상 생존할 확률은 72%

 

#전체 환자 중 생존확율이 50% 되는 지점

> which.min(abs(생명표$surv-0.5))

[1] 6

> 생명표[6,]

    nsubs nlost  nrisk nevent      surv        pdf    hazard    se.surv

5-6  1170   107 1116.5    125 0.5192585 0.05813463 0.1185958 0.01030422

        se.pdf  se.hazard

5-6 0.00503398 0.01058887

-> 5년이 지나면 surv 즉 생존할 확률이 약 52%가 됨을 알 수 있다.

 

 

8. 데이터 시각화 – 생존함수 추정치 그래프

 > plot(년[1:lt], 생명표[,5],type='s',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)


> plot(년[1:lt], 생명표[,5],type='o',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)         


 

 

9. 데이터 시각화 – 위험함수 추정치 그래프

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> mean(생명표$hazard)

[1] NA

> 생명표$hazard

 [1] 0.20821918 0.12353102 0.09440994 0.11991585 0.10804322 0.11859583

 [7] 0.10000000 0.11671924 0.10483042 0.11229947 0.15523466 0.17941953

[13] 0.14937759 0.11688312 0.13483146         NA

> mean(생명표$hazard,na.rm=T) #NA 값 제외한 평균

[1] 0.1294874

> plot(년[1:lt], 생명표[,7],type='s',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)

> abline(h=0.1294874)

 

> plot(년[1:lt], 생명표[,7],type='o',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)>

> abline(h=0.1294874)

 

 

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

반응형

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.3장 비모수적 방법을 이용한 생존함수의 비교  (0) 2017.01.26
제6.2장 비모수적 방법 - 2. 누적한계추정법  (0) 2017.01.26
제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우  (2) 2017.01.13
제5.1장 공분산 분석 - 1. 공변량이 하나인 경우  (0) 2017.01.06
제4장 범주형 자료의 분석 - 4.3 로짓분석  (0) 2017.01.04
Posted by 마르띤
,

제3장 나무모형 - 회귀나무모형

KNOU/2 데이터마이닝 2016. 10. 26. 14:08
반응형

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

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


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


> plot(fit.prune.Boston,uniform=T)

 


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

 

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

 

  

 

 


5) 목표변수의 적합값을 구하고 평가하자. 먼저 MSE 구하기

> Boston$medv.hat = predict(fit.prune.Boston,newdata=Boston,type='vector')

> mean((Boston$medv - Boston$medv.hat)^2)

 [1] 10.8643

> plot(Boston$medv,Boston$medv.hat,xlab='Observed values',ylab='Fitted values',xlim=c(0,50),ylim=c(0,50))
> abline(0,1)

 

결과 해석: MSE(Mean Squared Error, 평균오차제곱합)은 10.8643 이고, plot 그래프를 통해 cart 회귀나무모형의 적합값과 실제값의 일치도를 보이고 있음을 알 수 있다. 회귀나무 모형에서는 최종노드에는 적합값으로 관찰치들의 평균값을 사용하기 때문에 동일한 값을 가진 적합값이 많음을 확인할 수 있다.


 

6) 보스턴하우징 데이터를 훈련 데이터와 검증 데이터로 분할하여 회귀나무를 평가해 보자.

> set.seed(1234)

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

 

> 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 검증 데이터

 

#Fit a CART model by training data

 

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

> printcp(fit.train.Boston)

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

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

Root node error: 28893/354 = 81.618

n= 354

           CP nsplit rel error  xerror     xstd
1  0.45225287      0   1.00000 1.00649 0.099965
2  0.17404032      1   0.54775 0.64452 0.071140
3  0.06762461      2   0.37371 0.46448 0.058349
4  0.04776357      3   0.30608 0.41314 0.057818
5  0.02547541      4   0.25832 0.34693 0.051952
6  0.02433056      5   0.23284 0.32946 0.051911
7  0.01063840      6   0.20851 0.31645 0.051778
8  0.00718287      7   0.19787 0.31064 0.051954
9  0.00584084      8   0.19069 0.30989 0.052698
10 0.00558293     10   0.17901 0.30937 0.052600
11 0.00327437     11   0.17343 0.30188 0.052057
12 0.00323971     12   0.17015 0.29575 0.052015
13 0.00266789     13   0.16691 0.29202 0.051794
14 0.00242096     14   0.16424 0.29298 0.051791
15 0.00217440     15   0.16182 0.29230 0.051816
16 0.00164045   17   0.15748 0.28260 0.049406
17 0.00090517     18   0.15583 0.28542 0.049919
18 0.00081553     19   0.15493 0.28338 0.049884
19 0.00079750     20   0.15411 0.28347 0.049883
20 0.00067797     21   0.15332 0.28384 0.049878
21 0.00000000     22   0.15264 0.28332 0.049290

> which.min(fit.train.Boston$cptable[,'xerror'])
16
16

> fit.train.Boston$cptable[16]
[1] 0.001640451
> fit.train.Boston$cptable[16,]
          CP       nsplit    rel error       xerror         xstd
 0.001640451 17.000000000  0.157475033  0.282599587  0.049406453

>  0.28260 + 0.049406
[1] 0.332006

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

> fit.prune.train.Boston


 

 

 

> summary(fit.prune.train.Boston)
Call:
rpart(formula = medv ~ ., data = Boston.train, method = "anova",
    control = my.control)
  n= 354

            CP nsplit rel error    xerror       xstd
1  0.452252869      0 1.0000000 1.0064946 0.09996494
2  0.174040318      1 0.5477471 0.6445168 0.07114014
3  0.067624608      2 0.3737068 0.4644789 0.05834934
4  0.047763574      3 0.3060822 0.4131392 0.05781821
5  0.025475412      4 0.2583186 0.3469320 0.05195198
6  0.024330558      5 0.2328432 0.3294558 0.05191119
7  0.010638405      6 0.2085127 0.3164491 0.05177769
8  0.007182873      7 0.1978743 0.3106449 0.05195410
9  0.005840842      8 0.1906914 0.3098942 0.05269778
10 0.005582930     10 0.1790097 0.3093730 0.05259987

Variable importance
  lstat      rm   indus    crim     nox     age     dis     tax     rad ptratio
     26      19      14      11      11      11       3       2       1       1

Node number 1: 354 observations,    complexity param=0.4522529
  mean=22.15763, MSE=81.61758
  left son=2 (215 obs) right son=3 (139 obs)
  Primary splits: #improve는 개선도를 뜻한다.
      lstat   < 9.63     to the right, improve=0.4522529, (0 missing)
      rm      < 6.9715   to the left,  improve=0.4283989, (0 missing)
      indus   < 6.66     to the right, improve=0.2809072, (0 missing)
      ptratio < 19.9     to the right, improve=0.2685387, (0 missing)
      nox     < 0.6635   to the right, improve=0.2328252, (0 missing)
  Surrogate splits:
      indus < 7.625    to the right, agree=0.825, adj=0.554, (0 split)
      nox   < 0.519    to the right, agree=0.802, adj=0.496, (0 split)
      rm    < 6.478    to the left,  agree=0.794, adj=0.475, (0 split)
      crim  < 0.08547  to the right, agree=0.785, adj=0.453, (0 split)
      age   < 64.8     to the right, agree=0.782, adj=0.446, (0 split)

 

> plot(fit.prune.train.Boston,uniform=T,margin=0.1)

> text(fit.prune.train.Boston,use.n=T,col='blue',cex=0.8)

 

 

 


7) 예측치: 

이제 이 회귀나무를 활용하여 30% 검증 데이터에 적용시켜서 예측치를 구해보자. 그리고 그 예측치의 오차를 계산하기 위해 예측평균오차제곱합인 PMSE(Predicted Mean Square Error)를 구해보자

> medv.hat.test<-predict(fit.prune.train.Boston,newdata=Boston.test,type='vector')

> mean((Boston.test$medv-medv.hat.test)^2) #PSME,

[1] 13.95258

결과 해석: 예측평균오차제곱합은 13.95258로 계산되었다.

 

PMSE 계산시, Boston.train$medv가 아닌 Boston.test$medv임을 유의하자.

위에서 구한 관찰값과 적합값의 MSE(Mean Squared Error, 평균오차제곱합) 10.8643 과 비교하면 다소 높음을 알 수 있다. 



복습과 무한반복 연습이 필요해 보인다


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

반응형

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제4장 앙상블 모형 - 회귀앙상블모형 - 랜덤포레스트  (1) 2016.11.02
제4장 앙상블 모형  (0) 2016.11.02
제3장 나무모형 - 분류나무모형  (4) 2016.10.18
제2장 회귀모형 - 로지스틱 회귀모형 연습  (0) 2016.09.14
제2장 회귀모형 - 선형회귀 연습  (0) 2016.09.14
Posted by 마르띤
,

제2장 회귀모형 - 선형회귀 연습

KNOU/2 데이터마이닝 2016. 9. 14. 10:10
반응형

목표변수가 연속형인 경우 -> 선형 회귀모델, 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. 데이터 불러오기

> library(MASS)

> range(Boston$medv)

[1]  5 50

> stem(Boston$medv)

The decimal point is at the | 

4 | 006

6 | 30022245

8 | 1334455788567

10 | 2224455899035778899

12 | 013567778011112333444455668888899

14 | 0111233445556689990001222344666667

16 | 01112234556677880111222344455567888889

18 | 01222334445555667778899990011112233333444444555566666778889999

20 | 0000011111223333444455566666677888990001122222444445566777777788999

22 | 00000001222223344555666667788889999000011111112222333344566777788889

24 | 001112333444455566777888800000000123

26 | 24456667011555599

28 | 01244567770011466889

30 | 111357801255667

32 | 0024579011223448

34 | 679991244

36 | 01224502369

38 | 78

40 | 37

42 | 38158

44 | 084

46 | 07

48 | 358

50 | 0000000000000000

 

> i=which(Boston$medv==50)#본인 소유의 주택가격(중앙값)

> Boston[i,]

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

162 1.46336  0 19.58    0 0.6050 7.489  90.8 1.9709   5 403    14.7 374.43  1.73   50

163 1.83377  0 19.58    1 0.6050 7.802  98.2 2.0407   5 403    14.7 389.61  1.92   50

164 1.51902  0 19.58    1 0.6050 8.375  93.9 2.1620   5 403    14.7 388.45  3.32   50

167 2.01019  0 19.58    0 0.6050 7.929  96.2 2.0459   5 403    14.7 369.30  3.70   50

187 0.05602  0  2.46    0 0.4880 7.831  53.6 3.1992   3 193    17.8 392.63  4.45   50

196 0.01381 80  0.46    0 0.4220 7.875  32.0 5.6484   4 255    14.4 394.23  2.97   50

205 0.02009 95  2.68    0 0.4161 8.034  31.9 5.1180   4 224    14.7 390.55  2.88   50

226 0.52693  0  6.20    0 0.5040 8.725  83.0 2.8944   8 307    17.4 382.00  4.63   50

258 0.61154 20  3.97    0 0.6470 8.704  86.9 1.8010   5 264    13.0 389.70  5.12   50

268 0.57834 20  3.97    0 0.5750 8.297  67.0 2.4216   5 264    13.0 384.54  7.44   50

284 0.01501 90  1.21    1 0.4010 7.923  24.8 5.8850   1 198    13.6 395.52  3.16   50

369 4.89822  0 18.10    0 0.6310 4.970 100.0 1.3325  24 666    20.2 375.52  3.26   50

370 5.66998  0 18.10    1 0.6310 6.683  96.8 1.3567  24 666    20.2 375.33  3.73   50

371 6.53876  0 18.10    1 0.6310 7.016  97.5 1.2024  24 666    20.2 392.05  2.96   50

372 9.23230  0 18.10    0 0.6310 6.216 100.0 1.1691  24 666    20.2 366.15  9.53   50

373 8.26725  0 18.10    1 0.6680 5.875  89.6 1.1296  24 666    20.2 347.88  8.88   50

 

> boston=Boston[-i,] #최대값 50인 관측치 16개를 찾아 제거

> boston$chas = factor(boston$chas) #범주형으로 변경

> boston$rad = factor(boston$rad) #범주형으로 변경

> table(boston$rad)

1   2   3   4   5   6   7   8  24

19  24  37 108 109  26  17  23 127

 

> boston$chas <- as.factor(boston$chas)

> boston$rad <- as.factor(boston$rad)

> class(boston$rad);class(boston$chas)
[1] "factor"
[1] "factor"

 

[참고] 아래와 같은 방법으로 이용하면 모든 변수를 수치로 변경할 수 있다.

> for(i in 1:ncol(boston))if(!is.numeric(boston[,i])) boston[,i]=as.numeric(boston[,i])
> str(boston)
'data.frame':   490 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   : num  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    : num  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. 선형 회귀 모형 만들기

#선형회귀모형 만들기

> fit1 = lm(medv~.,data=boston) #목표변수 = medv, 선형회귀모형 함수, ~.는 목표 변수를 제외한 모든 변수를 입력변수로 사용

> summary(fit1)

 

  Call:

    lm(formula = medv ~ ., 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   

    chas1         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 ***

    rad24         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

 

또는 아래와 같은 방법도 가능하다

> names(boston)
 [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"   
 [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"  
> bn <- names(boston)

> f <- as.formula(paste('medv~',paste(bn[!bn %in% 'medv'],collapse='+')))
> f
medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +
    tax + ptratio + black + lstat
> fit2 <- lm(f,data=boston)
> summary(fit2)

Call:
lm(formula = f, 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 

  

 #가장 적절한 모형 선택 위한 변수 선택

> fit.step = step(fit1,direction='both') #both는 단계적 선택법 적용

  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

  + chas     1      0.63 6320.9 1295.0

  - 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

  + indus    1     11.00 6321.5 1293.1

  + chas     1      0.32 6332.2 1293.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) #최종모형, rad는 범주형 변수를 가변수로 변환한 것.#AIC가 가장 작은 변수가 단계적 선택법에 의해 변수들이 정의 됨

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

    rad24         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

-> 결과 해석: 최초 만든 회귀함수 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는 유의하지 않지만 다른 가변수가 유의하므로 제거되지 않고 여전히 모형에 포함된다. 

     

     R2은 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.all = lm(medv~.,data=boston)

> fit.step = step(fit.all, direction = "both")

  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

  + chas     1      0.63 6320.9 1295.0

  - 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

  + indus    1     11.00 6321.5 1293.1

  + chas     1      0.32 6332.2 1293.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

 

> names(fit.step)
 [1] "coefficients"  "residuals"     "effects"       "rank"        
 [5] "fitted.values" "assign"        "qr"            "df.residual" 
 [9] "contrasts"     "xlevels"       "call"          "terms"       
[13] "model"         "anova" 

 

 > fit.step$anova #최종모형에서 제거된 변수를 알 수 있다.

  Step Df   Deviance Resid. Df Resid. Dev      AIC

  1         NA         NA       469   6320.865 1295.031

  2  - chas  1  0.6261633       470   6321.491 1293.079

3 - indus  1 10.9964825       471   6332.487 1291.931

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


 

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

> yhat=predict(fit.step,newdata=boston,type='response') #목표값 예측 시, type='response'

> head(yhat) #예측된 값 산출

  1        2        3        4        5        6

  26.59831 24.00195 28.99396 29.60018 29.07676 26.41636

> plot(fit.step$fitted,boston$medv, xlim=c(0,50),ylim=c(0,50),xlab="Fitted",ylab="Observed")#실제값과 가까운지 평가

> abline(a=0,b=1) # or abline(0,1)

> mean((boston$medv-yhat)^2) #MSE

  [1] 12.92344

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

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


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

 

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

반응형

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제3장 나무모형 - 회귀나무모형  (0) 2016.10.26
제3장 나무모형 - 분류나무모형  (4) 2016.10.18
제2장 회귀모형 - 로지스틱 회귀모형 연습  (0) 2016.09.14
제2장 회귀모형 - 선형회귀, 로지스틱회귀  (0) 2016.09.14
1장 데이터과학과 데이터마이닝  (0) 2016.08.26
Posted by 마르띤
,

티스토리툴바