반응형

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

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


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



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


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

반응형
Posted by 마르띤
,
반응형

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

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

 

예제) 목표변수가 범주형 good, bad인 독일 신용평가 데이터를 이용하여 cart 방법을 이용한 의사결정나무 구축


1) 데이터 불러오기

> setwd('c:/Rwork')

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

> str(german)

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

$ check      : Factor w/ 4 levels "A11","A12","A13",..: 1 2 4 1 1 4 4 2 4 2 ...

$ duration   : int  6 48 12 42 24 36 24 36 12 30 ...

$ history    : Factor w/ 5 levels "A30","A31","A32",..: 5 3 5 3 4 3 3 3 3 5 ...

$ purpose    : Factor w/ 10 levels "A40","A41","A410",..: 5 5 8 4 1 8 4 2 5 1 ...

$ credit     : int  1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...

$ savings    : Factor w/ 5 levels "A61","A62","A63",..: 5 1 1 1 1 5 3 1 4 1 ...

$ employment : Factor w/ 5 levels "A71","A72","A73",..: 5 3 4 4 3 3 5 3 4 1 ...

$ installment: int  4 2 2 2 3 2 3 2 2 4 ...

$ personal   : Factor w/ 4 levels "A91","A92","A93",..: 3 2 3 3 3 3 3 3 1 4 ...

$ debtors    : Factor w/ 3 levels "A101","A102",..: 1 1 1 3 1 1 1 1 1 1 ...

$ residence  : int  4 2 3 4 4 4 4 2 4 2 ...

$ property   : Factor w/ 4 levels "A121","A122",..: 1 1 1 2 4 4 2 3 1 3 ...

$ age        : int  67 22 49 45 53 35 53 35 61 28 ...

$ others     : Factor w/ 3 levels "A141","A142",..: 3 3 3 3 3 3 3 3 3 3 ...

$ housing    : Factor w/ 3 levels "A151","A152",..: 2 2 2 3 3 3 2 1 2 2 ...

$ numcredits : int  2 1 1 1 2 1 1 1 1 2 ...

$ job        : Factor w/ 4 levels "A171","A172",..: 3 3 2 3 3 2 3 4 2 4 ...

$ residpeople: int  1 1 2 2 2 2 1 1 1 1 ...

$ telephone  : Factor w/ 2 levels "A191","A192": 2 1 1 1 1 2 1 2 1 1 ...

$ foreign    : Factor w/ 2 levels "A201","A202": 1 1 1 1 1 1 1 1 1 1 ...

$ y          : Factor w/ 2 levels "bad","good": 2 1 2 2 1 2 2 2 2 1 ...

> german$numcredits<-factor(german$numcredits)

> german$residence<-factor(german$residence)

> german$residpeople<-factor(german$residpeople)

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

[1] "factor"

[1] "factor"

[1] "factor"

 

2) cart 방법 적용

> library(rpart)

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

> fit.german<-rpart(y~.,data=german,method='class',control=my.control)

> fit.german #최초의 나무. 가지치기를 하지 않은 최대 크기의 나무 보기

n= 1000

 

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

* denotes terminal node

 

1) root 1000 300 good (0.300000000 0.700000000) 

2) check=A11,A12 543 240 good (0.441988950 0.558011050) 

.

. (너무 커서 중략)

.

    253) credit>=1273 10   0 good (0.000000000 1.000000000) *

  127) check=A14 122   1 good (0.008196721 0.991803279) *

함수 설명

 1. rpart.control:

- xval=10: 교타 타당성의 fold 개수, 디폴트는 10

- cp=0: 오분류율이 cp 이상으로 향상되지 않으면 이상 분할하지 않고 나무구조 생성을 멈춘다. cp값이 0이면 오분류값이 최소, 디폴트는 0.01

- minsplit=5: 노드를 분할하기 위해 필요한 데이터의 개수. 값보다 적은 수의 관측치가 있는 노드는 분할하지 않는다. 디폴트는 20

2. r.part

 - method=class: 나무 모형을 지정한다. anova 회귀나무, poisson 포아송 회귀나무, class 분류나무 exp 생존나무. 디폴트는 class

 - na.action=na.rpart: 목표변수가 결측치이면 전체 관측치를 삭제. 입력변수가 결측치인 경우에는 삭제하지 않는다.

 

결과 해석

중간노드를 분할하는 최소 자료의 수를 5개로 지정하였고, cp값은 0으로 하여 나무모형의 오분류값이 최소가 까지 분할을 진행하였다. 또한 10-fold 교차타당성을 수행하여 최적의 cp값을 찾도록 하였다. 나무가 너무나 관계로 중간 부분을 생략하였고, 용이한 모형 분석을 위해 가지치기를 해보자.

 

3) 나무를 줄이기 위한 가지치기 작업

> printcp(fit.german)

 

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

 

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

 

Root node error: 300/1000 = 0.3

 

n= 1000

 

               CP nsplit    rel error  xerror       xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556      8   0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000    13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

결과 해석

10-fold 교차타당성 방법에 의한 오분율(xerror)이 최소가 되는 값은 0.85667이며 이때의 cp값은 0.01임을 알 수 있다. 이 때 분리의 횟수가 13(nsplit=13)인 나무를 의미한다.

 

또는 아래와 같은 방법으로도 최소 오분류값(xerror)를 찾을 수 있다.

> names(fit.german)

[1] "frame"               "where"               "call"                "terms"             

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"           

> fit.german$cptable[,'xerror']

1         2         3         4         5         6         7         8         9

1.0000000 0.9466667 0.8633333 0.8700000 0.8866667 0.8800000 0.8566667 0.8700000 0.8733333

10        11        12        13        14        15        16        17

0.9200000 0.9100000 0.9200000 0.9200000 0.9700000 0.9766667 0.9766667 1.0133333

> which.min(fit.german$cptable[,'xerror'])

7

7

> fit.german$cptable[7,]

CP       nsplit     rel error       xerror        xstd

0.01000000 13.00000000  0.65333333  0.85666667  0.04606167

> fit.german$cptable[7]

[1] 0.01

> fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

[1] 0.01

> min.cp<-fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

> min.cp

[1] 0.01

> fit.prune.german<-prune(fit.german,cp=min.cp)


 


4) 오분류율이 최소인 cp(=0.011)을 찾았으니 이 값을 기준으로 가지치기를 시행하자.

> fit.prune.german<-prune(fit.german,cp=0.01)

> fit.prune.german


 

결과 해석

node), split, n, loss, yval, (yprob) 기준으로 첫번째 결과를 분석하면 다음과 같다.

노드, 분할점, 개수, …공부 필요

16) duration>=47.5 36   5 bad (0.8611111 0.1388889) *

 

duration 변수 47.5보다 경우, 전체 36(n)개를 bad(yval)로 분류하였고 그 중 5개(loss)가 good이다. 그리하여 bad로 분류되는 것은 31/36 = 0.8611111로 표기하게 되고, 5개의 loss는 5/36 = 1388889 로 그 확률을 볼 수 있다. 아래 plot에서는 bad 31/5 표기

 

376) property=A123,A124 20   3 bad (0.8500000 0.1500000) *

377) property=A121,A122 45  14 good (0.3111111 0.6888889) *

 

property a123(car), a124(unknown / no property) 경우 전체 20개를 bad로 분류하였고 3개의 loss 즉 good (3/20 = 0.15)로 분류하였다. 아래 plot에서는 bad 17/3 표기 

property a121(real estate), a122(building society savings agreement) 경우에는 전체 45개를 good으로 분류하였고 14개의 loss 즉 bad로 분류 (14/45=0.3111111), 아래 plot에서는 good 14/31 표기


<< 17.6.18(일)>> 해석 부분 내용 추가

duration > = 22.5인 경우, 전체 고객은 237명이고, 이 중 신용도가 나쁜 사람의 비율은 56.5%이고 좋은 사람의 비율은43.5%로 103명이다. 따라서 duration > 22.5 그룹은 bad로 분류된다.



가지치기를 한 모형을 그림으로 나타내는 함수는 아래와 같다.

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

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

 

왼쪽 가지의 가장 아랫부분의 분할점인 ‘purpose=acdeghj’purpose 변수의 범주값 중에서 알파벳 순서로, 1(=a), 3(=c), 4(=d), 5(=e), 7(=g), 8(=h), 10(=j)번째 범주값을 의미하며, fit.prune.german에서 각각 A40,A410,A42,A43,A45,A46,A49 임을 알 수 있다.

 

34) purpose=A40,A410,A42,A43,A45,A46,A49 137  52 bad (0.6204380 0.3795620) *



<< 17.6.18(일)>> 해석 부분
 내용 추가

가장 우측의 duration > = 11.5가 아닌 경우, 신용다가 나쁜 / 좋은 사람의 비율은 9명 / . 4명이고, 신용도가 좋은 good으로 분류된다.

 

 

5) 나무수를 더 줄여보자.

> printcp(fit.german)

 

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

 

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

 

Root node error: 300/1000 = 0.3

 

n= 1000

 

                CP nsplit   rel error  xerror      xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556       0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000     13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

5번째 단계이며 분리의 횟수가 8(nsplit=8)인 나무는 교차타당성 오분류율이 0.88667로 최소는 아니지만 7번째 단계의 분리의 횟수 13회 나무 가지의 최소 오분류율 0.85667과는 크게 차이가 나지 않는다. 그리고 최소 오분류율 표준편차의 1배 범위(0.88667  <  0.85667 + 0.046062) 있다. 이런 경우에는 5번째 단계이며 분리의 횟수가 8 나무를 선택하는 경우도 있다.

 

5번째 단계이며 분리 횟수가 8 cp 0.0155556의 반올림 값 0.016 적용하여 다시 가지치기

> fit.prune.german<-prune(fit.german,cp=0.016)

> fit.prune.german

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

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

 

 


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

 > fit.prune.german<-prune(fit.german,cp=0.01) 

 > pred.german=predict(fit.prune.german,newdata=german,type='class')

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

 > tab

           Predicted

  Actual  bad  good

    bad  180   120

   good   76   624

 

함수 설명

predict(fit.prune.german,newdata=german,type='class'), type = class 분류나무의 집단값 예측결과, 회귀나무라면 type = vector라고 해야 한다.


결과 해석

실제 good인데 good으로 예측한 것이 624, 실제 bad인데 bad로 예측한 것이 180

따라서 오분류율은 {1000 – (624+180)} / 1000 = 19.6%


R코드를 이용하면 1-sum(diag(tab)) / sum(tab)



7) 마지막으로 독일신용평가데이터를 훈련데이터와 검증 데이터로 분할하여 분류나무를 평가해보자.

> set.seed(1234)

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

> german.train=german[i,]

> german.test=german[-i,]

> fit.german<-rpart(y~.,data=german.train,method='class',control=my.control)

> printcp(fit.german)

 

Classification tree:

  rpart(formula = y ~ ., data = german.train, method = "class",

        control = my.control)

 

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings     telephone 

 

Root node error: 201/700 = 0.28714

 

n= 700

 

CP  nsplit   rel error  xerror     xstd

1  0.05721393      0   1.00000 1.00000 0.059553

2  0.03482587      2   0.88557 1.00498 0.059641

3  0.02985075      5   0.78109 1.00000 0.059553

4  0.01990050    6   0.75124 0.95025 0.058631

5  0.01741294      8   0.71144 0.96020 0.058822

6  0.01492537     10   0.67662 1.00000 0.059553

7  0.01243781     14   0.61692 1.00000 0.059553

8  0.00995025     17   0.57711 1.00995 0.059728

9  0.00746269     35   0.39303 1.03980 0.060238

10 0.00621891     46   0.30846 1.06965 0.060722

11 0.00497512     50   0.28358 1.04975 0.060402

12 0.00331675     58   0.24378 1.09950 0.061181

13 0.00248756     61   0.23383 1.11940 0.061474

14 0.00124378     69   0.21393 1.14925 0.061894

15 0.00099502     73   0.20896 1.14925 0.061894

16 0.00000000     78   0.20398 1.14925 0.061894

> fit.prune.german<-prune(fit.german,cp=0.02)

> fit.prune.german

 

> p.german.test=predict(fit.prune.german,newdata=german.test,type='class')

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

> tab

            Predicted

Actual  bad good

      bad   34   65

      good  14  187

> 1-sum(diag(tab))/sum(tab) #오분류율

[1] 0.2633333 


 

 

 


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

 

반응형
Posted by 마르띤
,
반응형

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

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


독일신용평가 데이터 셋

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


변수명

속성

변수 설명

check

범주형

자유예금형태
Status of existing checking account

duration

수치형

기간
Duration in month

history

범주형

과거신용정보
Credit history

purpose

범주형

목적
Purpose

credit

수치형

신용대출금액
Credit amount

savings

범주형

저축예금/채권
Savings account/bonds

employment

범주형

현직장 재직기간
Present employment since

installment

수치형

가처분소득 대비 적금비율
Installment rate in percentage of disposable income

personal

범주형

결혼상황 성별
Personal status and sex

debtors

범주형

여타 채무/채권
Other debtors / guarantors

residence

수치형

거주기간
Present residence since

property

범주형

재산
Property

age

수치형

나이
Age in years

others

범주형

여타적금
Other installment plans

housing

범주형

주거형태
Housing

numcredits

수치형

해당 은행 신용계좌
Number of existing credits at this bank

job

범주형

직업
Job

residpeople

수치형

부양가족수
Number of people being liable to provide maintenance for

telephone

범주형

전화소유
Telephone

foreign

범주형

외국인 노동자 여부
foreign worker

y

범주형

신용등급 양호 또는 불량
credit:Good or Bad

 

 



1. 데이터 불러오기

> setwd('c:/Rwork')

> german<-read.table('germandata.txt')

> head(german,2) # 값들의 변수명이 없음.


> names<-c("check","duration","history","purpose","credit","savings","employment","installment",         "personal",      "debtors",       "residence",     "property",       "age",   "others",         "housing",       "numcredits",    "job",   "residpeople",    "telephone",     "foreign"         ,"y")

> colnames(german)<-names

> head(german,2)


> german$y<-factor(german$y,levels=c(1,2),labels=c('good','bad'))

> head(german,2)


> summary(german)


#  residence,numcredits,residpeople 실제 범주형이지만 수치형으로 인식. 범주형으로 변환 필요 

> class(german$residence) #integer 수치형

[1] "integer"

> class(german$check) #factor 범주형

[1] "factor"

> german$residence = factor(german$residence)

> german$numcredits = factor(german$numcredits)

> german$residpeople = factor(german$residpeople)

> class(german$residence) #factor 변환

[1] "factor"

> class(german$numcredits) #factor 변환

[1] "factor"

> class(german$residpeople) #factor 변환

[1] "factor"

> table(german$residence)

1   2   3   4

130 308 149 413

> german$y<-ifelse(german$y=='good',1,0) #반응 good 1, bad 2 변환

 


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

> fit.all = glm(y~.,family = binomial,data=german) #로지스틱 회귀 분석

 

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

> gmn<-names(german)
> f<-as.formula(paste('y~',paste(gmn[!gmn%in%y],collapse='+')))
> fit.all.1<-glm(f,family = binomial, data=german)

> fit.step = step(fit.all, direction='both') #단계적 선택방법

Start:  AIC=993.44

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

  employment + installment + personal + debtors + residence +

  property + age + others + housing + numcredits + job + residpeople +

  telephone + foreign

 

Df Deviance     AIC

- job          3   888.00  988.00

- numcredits   3   890.25  990.25

- property     3   890.70  990.70

- residpeople  1   888.52  992.52

- age          1   889.37  993.37

- telephone    1   889.40  993.40

<none>             887.44  993.44

- employment   4   895.48  993.48

- housing      2   891.63  993.63

- residence    3   894.74  994.74

- debtors      2   894.80  996.80

- others       2   895.71  997.71

- personal     3   897.80  997.80

- foreign      1   894.16  998.16

- credit       1   895.07  999.07

- duration     1   896.25 1000.25

- installment  1   900.81 1004.81

- savings      4   908.55 1006.55

- history      4   911.01 1009.01

- purpose      9   922.07 1010.07

- check        3   957.33 1057.33

 

Step:  AIC=988

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

  employment + installment + personal + debtors + residence +

  property + age + others + housing + numcredits + residpeople +

  telephone + foreign

 

Df Deviance     AIC

- numcredits   3   890.85  984.85

- property     3   891.21  985.21

- residpeople  1   889.08  987.08

- employment   4   895.67  987.67

<none>             888.00  988.00

- housing      2   892.01  988.01

- age          1   890.05  988.05

- telephone    1   890.34  988.34

- residence    3   895.32  989.32

- debtors      2   895.25  991.25

- personal     3   898.31  992.31

- others       2   896.49  992.49

- foreign      1   894.77  992.77

+ job          3   887.44  993.44

- credit       1   895.72  993.72

- duration     1   897.14  995.14

- installment  1   901.56  999.56

- savings      4   909.71 1001.71

- history      4   911.44 1003.44

- purpose      9   922.89 1004.89

- check        3   957.60 1051.60

 

Step:  AIC=984.85

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

  employment + installment + personal + debtors + residence +

  property + age + others + housing + residpeople + telephone +

  foreign

 

Df Deviance     AIC

- property     3   894.03  982.03

- employment   4   898.02  984.02

- residpeople  1   892.07  984.07

- age          1   892.85  984.85

<none>             890.85  984.85

- housing      2   895.09  985.09

- telephone    1   893.29  985.29

- residence    3   898.52  986.52

+ numcredits   3   888.00  988.00

- debtors      2   898.27  988.27

- personal     3   901.17  989.17

- others       2   899.85  989.85

- foreign      1   898.00  990.00

+ job          3   890.25  990.25

- credit       1   898.64  990.64

- duration     1   899.76  991.76

- installment  1   904.66  996.66

- history      4   911.95  997.95

- savings      4   912.53  998.53

- purpose      9   926.15 1002.15

- check        3   959.38 1047.38

 

Step:  AIC=982.03

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

  employment + installment + personal + debtors + residence +

  age + others + housing + residpeople + telephone + foreign

 

Df Deviance     AIC

- residpeople  1   895.11  981.11

- employment   4   901.94  981.94

- telephone    1   895.95  981.95

<none>             894.03  982.03

- age          1   896.10  982.10

- housing      2   898.15  982.15

- residence    3   901.53  983.53

+ property     3   890.85  984.85

+ numcredits   3   891.21  985.21

- personal     3   903.97  985.97

- debtors      2   902.35  986.35

- foreign      1   901.07  987.07

+ job          3   893.45  987.45

- others       2   903.55  987.55

- credit       1   902.94  988.94

- duration     1   903.85  989.85

- installment  1   908.62  994.62

- savings      4   915.22  995.22

- history      4   915.59  995.59

- purpose      9   930.66 1000.66

- check        3   964.51 1046.51

 

Step:  AIC=981.11

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

  employment + installment + personal + debtors + residence +

  age + others + housing + telephone + foreign

 

Df Deviance     AIC

- employment   4   903.04  981.04

- age          1   897.04  981.04

<none>             895.11  981.11

- telephone    1   897.12  981.12

- housing      2   899.31  981.31

+ residpeople  1   894.03  982.03

- residence    3   902.80  982.80

- personal     3   904.04  984.04

+ property     3   892.07  984.07

+ numcredits   3   892.19  984.19

- debtors      2   903.15  985.15

- foreign      1   902.06  986.06

+ job          3   894.59  986.59

- others       2   904.70  986.70

- credit       1   903.73  987.73

- duration     1   904.80  988.80

- installment  1   909.03  993.03

- savings      4   916.06  994.06

- history      4   916.94  994.94

- purpose      9   932.01 1000.01

- check        3   965.87 1045.87

 

Step:  AIC=981.04

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

  installment + personal + debtors + residence + age + others +

  housing + telephone + foreign

 

Df Deviance     AIC

- age          1   904.91  980.91

<none>             903.04  981.04

+ employment   4   895.11  981.11

- telephone    1   905.28  981.28

- housing      2   907.58  981.58

+ residpeople  1   901.94  981.94

- residence    3   910.50  982.50

+ property     3   899.28  983.28

+ numcredits   3   900.64  984.64

- foreign      1   909.67  985.67

- debtors      2   912.24  986.24

+ job          3   902.89  986.89

- personal     3   915.04  987.04

- others       2   913.21  987.21

- duration     1   911.34  987.34

- credit       1   911.50  987.50

- installment  1   917.92  993.92

- savings      4   925.25  995.25

- history      4   925.74  995.74

- purpose      9   939.70  999.70

- check        3   975.57 1047.57

 

Step:  AIC=980.91

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

  installment + personal + debtors + residence + others + housing +

  telephone + foreign

 

Df Deviance     AIC

<none>             904.91  980.91

+ age          1   903.04  981.04

+ employment   4   897.04  981.04

- telephone    1   907.69  981.69

+ residpeople  1   903.95  981.95

- housing      2   910.11  982.11

- residence    3   912.96  982.96

+ property     3   901.18  983.18

+ numcredits   3   902.60  984.60

- foreign      1   911.56  985.56

- debtors      2   914.35  986.35

- others       2   914.61  986.61

+ job          3   904.63  986.63

- credit       1   913.18  987.18

- personal     3   917.50  987.50

- duration     1   914.06  988.06

- installment  1   919.35  993.35

- savings      4   927.70  995.70

- history      4   928.79  996.79

- purpose      9   940.82  998.82

- check        3   978.40 1048.40

 

> fit.step$anova #제거된 변수 보기

Step Df  Deviance Resid. Df Resid. Dev      AIC

1               NA        NA       947   887.4372 993.4372

2         - job  3 0.5588674       950   887.9960 987.9960

3  - numcredits  3 2.8582392       953   890.8543 984.8543

4    - property  3 3.1777611       956   894.0320 982.0320

5 - residpeople  1 1.0747973       957   895.1068 981.1068

6  - employment  4 7.9298736       961   903.0367 981.0367

7         - age  1 1.8704615       962   904.9072 980.9072

> summary(fit.step) #최종모델

 

Call:

  glm(formula = y ~ check + duration + history + purpose + credit +

      savings + installment + personal + debtors + residence +

      others + housing + telephone + foreign, family = binomial,

      data = german)

 

Deviance Residuals:

  Min       1Q   Median       3Q      Max 

-2.7904  -0.7290   0.3885   0.6911   2.1780 

 

Coefficients:

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

(Intercept)   -9.736e-01  7.032e-01  -1.385 0.166204   

checkA12       3.863e-01  2.136e-01   1.809 0.070468 . 

checkA13       1.055e+00  3.636e-01   2.902 0.003714 **

checkA14       1.782e+00  2.308e-01   7.721 1.15e-14 ***

duration      -2.726e-02  9.034e-03  -3.018 0.002546 **

historyA31     1.290e-01  5.297e-01   0.244 0.807596   

historyA32     8.608e-01  4.104e-01   2.097 0.035956 * 

historyA33     9.975e-01  4.675e-01   2.133 0.032889 * 

historyA34     1.564e+00  4.329e-01   3.612 0.000303 ***

purposeA41     1.591e+00  3.684e-01   4.320 1.56e-05 ***

purposeA410    1.397e+00  7.732e-01   1.806 0.070849 . 

purposeA42     6.766e-01  2.529e-01   2.675 0.007467 **

purposeA43     8.867e-01  2.443e-01   3.629 0.000284 ***

purposeA44     5.231e-01  7.546e-01   0.693 0.488206   

purposeA45     1.335e-01  5.388e-01   0.248 0.804301   

purposeA46    -2.006e-01  3.883e-01  -0.517 0.605426   

purposeA48     2.060e+00  1.202e+00   1.714 0.086523 . 

purposeA49     7.396e-01  3.318e-01   2.229 0.025820 * 

credit        -1.230e-04  4.314e-05  -2.852 0.004351 **

savingsA62     3.126e-01  2.805e-01   1.115 0.264984   

savingsA63     4.303e-01  3.887e-01   1.107 0.268291   

savingsA64     1.396e+00  5.184e-01   2.692 0.007106 **

savingsA65     1.004e+00  2.606e-01   3.852 0.000117 ***

installment   -3.218e-01  8.621e-02  -3.733 0.000189 ***

personalA92    2.159e-01  3.754e-01   0.575 0.565268   

personalA93    8.302e-01  3.672e-01   2.261 0.023766 * 

personalA94    3.551e-01  4.434e-01   0.801 0.423122   

debtorsA102   -4.978e-01  4.005e-01  -1.243 0.213967   

debtorsA103    1.074e+00  4.205e-01   2.555 0.010628 * 

residence2    -7.181e-01  2.796e-01  -2.568 0.010223 * 

residence3    -3.929e-01  3.246e-01  -1.210 0.226104   

residence4    -2.893e-01  2.806e-01  -1.031 0.302546   

othersA142     5.959e-02  4.061e-01   0.147 0.883344   

othersA143     6.787e-01  2.355e-01   2.882 0.003955 **

housingA152    5.098e-01  2.271e-01   2.245 0.024799 * 

housingA153    2.464e-01  3.288e-01   0.749 0.453710   

telephoneA192  3.051e-01  1.838e-01   1.660 0.096958 . 

foreignA202    1.439e+00  6.253e-01   2.301 0.021383 * 

  ---

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

 

(Dispersion parameter for binomial family taken to be 1)

 

Null deviance: 1221.73  on 999  degrees of freedom

Residual deviance:  904.91  on 962  degrees of freedom

AIC: 980.91

 

Number of Fisher Scoring iterations: 5



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

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

 

 

단계적선택법의 AIC 980.91

[참고] 후진소거법의 AIC 980.91

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

Step:  AIC=980.91

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

    installment + personal + debtors + residence + others + housing +

    telephone + foreign

 

              Df Deviance     AIC

<none>             904.91  980.91

- telephone    1   907.69  981.69

- housing      2   910.11  982.11

- residence    3   912.96  982.96

- foreign      1   911.56  985.56

- debtors      2   914.35  986.35

- others       2   914.61  986.61

- credit       1   913.18  987.18

- personal     3   917.50  987.50

- duration     1   914.06  988.06

- installment  1   919.35  993.35

- savings      4   927.70  995.70

- history      4   928.79  996.79

- purpose      9   940.82  998.82

- check        3   978.40 1048.40

 

> fit.step.back$anova #제거된 변수 보기

           Step Df  Deviance Resid. Df Resid. Dev      AIC

1               NA        NA       947   887.4372 993.4372

2         - job  3 0.5588674       950   887.9960 987.9960

3  - numcredits  3 2.8582392       953   890.8543 984.8543

4    - property  3 3.1777611       956   894.0320 982.0320

5 - residpeople  1 1.0747973       957   895.1068 981.1068

6  - employment  4 7.9298736       961   903.0367 981.0367

7         - age  1 1.8704615       962   904.9072 980.9072

 

 

[참고] 전진선택법 AIC : 993.44

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

Start:  AIC=993.44

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

    employment + installment + personal + debtors + residence +

    property + age + others + housing + numcredits + job + residpeople +

telephone + foreign

 

> fit.step.forward$anova #제거된 변수 보기

  Step Df Deviance Resid. Df Resid. Dev      AIC

1      NA       NA       947   887.4372 993.4372

 



 

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

> p = predict(fit.step, newdata=german,type='response')

> threshold = 0.5 #cutoff기준 0.5 정함

> yhat = ifelse(p>threshold,1,0)

> head(yhat)

  1 2 3 4 5 6

  1 0 1 1 0 1

> class.tab = table(german$y,yhat,dnn=c("Actual","Predicted"))#실값과 예측값 배열

> class.tab

       Predicted

Actual   0   1

     0 158 142

     1  82 618

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


 

4. 예측력 측도

> sum(german$y==yhat)/length(german$y) #Prediction Accuracy 예측정확도

[1] 0.776

> sum(german$y!=yhat)/length(german$y) #Misclassification Rate 오분류율

[1] 0.224

> class.tab[1,1]/apply(class.tab,1,sum)[1] #Specificity 특이도

0

0.5266667

> class.tab[2,2]/apply(class.tab,1,sum)[2] #Sensitivity 민감도

1

0.8828571

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


 

5. ROC 곡선 및 AUC 생성

 

> library(ROCR)

> pred<-prediction(p,german$y)

> perf<-performance(pred,'tpr','fpr') #민감도와 1-특이도 계산 과정

> plot(perf,lty=1,col=2,xlim=c(0,1),ylim=c(0,1),xlab='1-Specificity',ylab='Sensitivity',main='ROC Curve')

> lines(x=c(0,1),y=c(0,1),col='grey')


> performance(pred,'auc')@y.values #면적 계산

[[1]]

[1] 0.8312286


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

반응형
Posted by 마르띤
,
반응형

목표변수가 연속형인 경우 -> 선형 회귀모델, 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 대각선을 중심으로 모여 있으면 예측력이 좋다고 있다

 

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

반응형
Posted by 마르띤
,
반응형

문제: A제품과 B 제품이 있는데, B 제품의 제품 클레임이 A제품보다 많다고 한다. A, B 두 제품의 품질 간에 유의한 차이가 있는지 독립표본 두 모평균 비교를 수행하라.

 

1. 데이터를 불러오자.

> setwd('c:/Rwork')

> claim<-read.csv('claim.csv',header = T)

> boxplot(ppm~제품, data=claim) 

 

##A B 나누어 살펴보기

> library(dplyr)

> claim.A <- claim %>% filter(제품 == "A")

> claim.B <- claim %>% filter(제품 == "B")

> nrow(claim.A) ; nrow(claim.B)

[1] 105

[1] 40

 

> claim.A$제품 %>% table #질적-명목형 변수

A   B

105   0

> claim.A$ppm %>% summary #양적-이산형 변수

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

3.0    83.0   211.0   259.7   385.0   994.0

> claim.A$ppm %>% sd

[1] 217.4249


> claim.B$ppm %>% summary #양적-이산형 변수

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

29.00   71.25  145.50  272.60  345.20 1104.00

> claim.B$ppm %>% sd

[1] 279.8323

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

> claim.A$ppm %>% hist

> claim.B$ppm %>% hist


-> 해석: A의 평균 ppm 259, 표준편차는217, B는 평균이 272, 표준편차 279. BoxplotHistogram만으로는 확연히 B의 클레임 ppm이 높다고는 볼 수 없다.

 


2. R을 이용하여 통계 분석을 해보자.

두 집단의 평균이 같은지를 비교하여 두 집단의 차이를 평가하는 경우 two-sample t-test를 사용하는데, 쌍을 이룬 두 변수 간에 차이의 평균이 0인지 검정하는 paired t-test와는 달리 서로 독립적인 두 집단의 평균의 차이가 0인지를 검정한다. Two-sample t-test는 다음 순서를 따른다.

 

두 집단의 분산이 같은지 검정하자.

> var.test(ppm~제품, data=claim)

F test to compare two variances

 

data:  ppm by 제품

F = 0.6037, num df = 104, denom df = 39, p-value = 0.04555

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

95 percent confidence interval:

  0.3454412 0.9901513

sample estimates:

  ratio of variances

0.6037024

-> 해석:

귀무가설 H0: 두 모형의 분산이 같다.

대립가설 H1: 두 모형의 분산이 다르다.

F = 0.6037, num df = 104, denom df = 39, p-value = 0.04555

두 분산의 비율이 0.6배이고, p값이 0.05보다 작으므로 두 분산이 같다는 귀무가설을 기각, 즉 두 모집단이 다르다, 분산이 다른 경우 에는 welch t.test를 한다.

 

분산이 다른 경우 에는 welch t.test 실시한다.

> t.test(ppm~제품,data=claim)

 

Welch Two Sample t-test

 

data:  ppm by 제품

t = -0.26303, df = 57.854, p-value = 0.7935

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

95 percent confidence interval:

  -111.13697   85.32268

sample estimates:

  mean in group A mean in group B

259.7429        272.6500

-> 해석:

귀무가설 H0: μ1-μ2 = 0,

대립가설 H1: μ1-μ2 ≠ 0

검정통계량 T(df=57.854)=-0.26303

p-value 0.7935 > 0.05


평균이 같다는 귀무가설을 기각하지 못한다. 즉 제품 A와 제품B간에는 유의수준 5%에서 품질 차이가 없다고 결론지을 수 있다.



아직 R과 통계 공부 걸음마 단계인데 어서 익숙해지고 싶다.

반응형
Posted by 마르띤
,
반응형

회귀분석(regression analysis): 독립변수와 종속변수 간의 함수 관계를 규명하는 통계적인 분석방법


Ŷ=f(X)+ε


 - 독립변수(independent variable) 또는 설명변수(explanatory variable): 다른 변수에 영향을 주는 변수, 흔히 Y= β0 + β1X 공식에서 X

 - 종속변수(dependent variable) 또는 반응변수(response variable): 독립 변수에 의해 영향을 받는다는 변수, 흔히 Y= β0 + β1X 공식에서 Y

 

회귀(回歸)라는 말은다시 본디의 자리로 돌아온다라는 뜻으로 통계 분석에 처음 사용한 사람은 영국의 우생학자 Galton. 완두콩 실험을 통해 부모콩의 무게를 X, 자식콩의 무게를 Y축으로 산점도를 그리자, 이들의 관계식은 양이 관계이나 1보다 작아서 자식의 무게는 평균 무게로 회귀하려는 경향이 있다는 사실을 발견하고 이를 회귀(regression)으로 표현. 당시에 Galton 연구실에서 일하던 동료 연구원 Karl Pearson 이를 계량적으로 처음으로 분석하여 발표.

 

1. 데이터를 불러와서 산점도 그래프를 그리기

> market=read.table('market-1.txt',header=T)

> head(market,3)

NUMBER X  Y

1      1 4  9

2      2 8 20

3      3 9 22

> plot(market$X,market$Y,xlab='광고료',ylab='총판매액',pch=19)

> title('광고료와 판매액의 산점도')


 

2. 단순 회귀 분석 실시

> market.lm=lm(Y~X,data=market)

> summary(market.lm)

 

 

해석

추정값은 Coefficients: Estimate에서 확인. 추정된 회귀식은 Ŷ=-2.27 + 2.6 X

 


3. 산점도 위에 회귀직선을 그리자 - 회귀선의 추정

> abline(market.lm)

> identify(market$X,market$Y)

[1]  4  5 10

# Identify는 재미있는 함수인데, 본 함수를 입력하고 마우스로 점을 클릭하면 그림처럼 값을 알 수 있다.


> xbar = mean(market$X)

> ybar = mean(market$Y)

> xbar

[1] 8

> ybar

[1] 18.6

> points(xbar,ybar,pch=17,cex=2.0,col='RED')

> text(xbar,ybar,"(8,18.6)")

> fx <- "Y-hat = -2.27 + 2.6*X "

> text(locator(1),fx) #locator(1) 마우스로 클릭하면서 지정


 


4. 회귀식 특징

> names(market.lm)

[1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values"

[6] "assign"        "qr"            "df.residual"   "xlevels"       "call"        

[11] "terms"         "model"       

> market.lm$resid

1          2          3          4          5          6          7          8

0.8347826  1.4000000  0.7913043 -3.6000000 -1.6000000  0.9652174  4.6173913  1.1826087

9         10

-3.3826087 -1.2086957

> resid=market.lm$residual

> sum(resid) #특징1. 잔차의합은 0이다

[1] 0

> sum(market$X*resid) #특징2. 잔차들의 Xi 의한 가중합은 0이다

[1] 2.220446e-15

> sum(market.lm$fitted*resid) #특징3. 잔차들의 Yi 의한 가중합은 0이다

[1] -1.24345e-14

> names(market.lm)

 [1] "coefficients"  "residuals"     "effects"       "rank"         

 [5] "fitted.values" "assign"        "qr"            "df.residual"  

 [9] "xlevels"       "call"          "terms"         "model"        

> sum(market.lm$fitted.values)

[1] 186

> sum(market$Y)

[1] 186

#특징4. 추정값Yhat의 값과 관찰값Yi의 값은 같다.


5. 회귀모형의 정도

 - 산점도 위에 회귀직선을 그려 회귀선의 정도를 대략 짐작할 있으나, 이러한 경우는 독립변수가 하나인 경우에만 유용하게 쓰일 있다. 추정된 회귀선의 정도를 측정하는 여러 가지 측도(measure) 중에서 널리 이용되는 세가지를 알아보자

 ① 분산분석표에 의한 F-검정

 ② 결정계수

 ③ 추정값의 표준오차 

 ④ 상관계수와 결정계수


① 분산분석표에 의한 F-검정

요인

자유도

제곱합

평균제곱

F0

회귀

1

SSR(회귀제곱합)

MSR=SSR

MSR/MSE

잔차

n-2

SSE(잔차제곱합)

MES=SSE/n-2

 

n-1

SST( 제곱합)

 

 

- SST(Total sum of squares): 제곱합

- SSR(Sum of squares due to regression): 회귀제곱합, 설명되는 편차

- SSE(Sum of squares due to residual errors): 잔차제곱합, 설명되지 않는 편차

 

귀무가설 H0 : β1 = 0

대립가설 H1 : β1 0

 

F0 > F(1, n-2 ; α )이면 귀무가설 H0 : β1 = 0 기각하고, 회귀직선이 유의하다고 말한다. R분석 결과에서는 검정통계량 F0 대한 유의확률 p값이 제공된다. p < 유의확률 α이면 귀무가설 H0 : β1 = 0 기각한다

 

분산분석표

> anova(market.lm) #분산분석표

 

요인

자유도

제곱합

평균제곱

F0

회귀

1 = 1

SSR(회귀제곱합)

= 313.04

MSR=SSR

= 313.04

MSR/MSE

= 45.24

잔차

10-2 = 8

SSE(잔차제곱합)

= 55.36

MSE(잔차 평균제곱)  = SSE/n-2= 6.92

 

10-1 = 9

SST( 제곱합)

= 368.4

 

 

p값은0.0001487p < 유의확률 α 이므로 귀무가설  H0 : β1 = 0 기각한다, 따라서 회귀식은 유의하다.

 

> qf(0.95,1,8)

[1] 5.317655

유의수준 α =0.05에서 F-기각역 F(1,8;0.05)의 값은 5.32, " F0 = 45.24 > 5.32"이므로 귀무가설 기각, 회귀선은 유희하다

 

> 1-pf(45.25,1,8)

[1] 0.0001485485

유의확률 p값을 이용한 검정은 0.0001485485 이 값이 주어진 유의수준 α =0.05 보다 작을수록 귀무가설을 기각한다.

 


결정계수


R2=SSR/SST =1-SSE/SST


R2 결정계수(coefficient of determination)라고 부른다.


R2=SSR/SST = 313.04/368.4= 84.97%


이는 총변동 주에서 회귀직선에 의하여 설명되는 부분이 84.97%라는 의미로서, 추정된 회귀선의 정도가 높다는 것을 있다.

 

R 통해서도 있는데, Multiple R-squared 0.8497임을 있다.

> summary(market.lm)  



 

 추정값의 표준오차

선형회귀모형 Y= β0 + β1X + ε 을 표본의 자료로부터 적합시킬 때,


Y의 기댓값은 E(Y) = β0 + β1X, 분산은 σ2로 가정,


하였다. 따라서 Y의 측정값들이 회귀선 주위에 가깝게 있다면 σ의 추정값은 작아질 것이다. 

분산분석표에서 잔차평균제곱 MSE는 σ2의 불편추정량이 된다. 따라서 MSE의 제곱근을 추정값의 표준오차(standard error of estimate)라고 부르며, 다음과 같이 표현한다.



SY•X = SQRT(MSE) = SQRT(SSE/n-2) = 2.63


 

R 통해서도 있는데, Residual standard error 2.631임을 있다.

> summary(market.lm) 


 


상관계수와 결정계수

상관계수는 연속인 변수 간의 선형관계(linear relationship) 어느 정도인가를 재는 측도로서, 단순 회귀 분석에서는 상관계수 r 다음과 같이 구할 있다.


r = ±SQRT(R2 )


즉 상관계수는 결정계수 R2 제곱근이며, 만약 추정된 회귀선의 기울기 b1 양이면 양의 상관계수를 갖고, 기울기b1 음이면 음의 상관계수를 갖는다. 회귀식 Ŷ=-2.27 + 2.6 X에서 b1 2.6으로 양이므로 상관계수는 0.8497 sqrt , 0.9217임을 있다.

반응형
Posted by 마르띤
,
반응형

데이터 마이닝 기법의 구분

데이터마이닝에서 사용되는 기법은 크게 지도학습(supervised learning)과 자율학습(unsupervised learning)으로 나눌 수 있음

 

 - 지도학습의 목표는 입출력 간의 관계를 결정하는 시스템에 대한 유용한 근사 시스템을 구하는 것으로 정의할 수 있음. Y = aX + b

- 자율학습에서는교사의 역할에 해당하는 실제 출력값이 존재하지 않음. 따라서 데이터에 존재하는 여러 가지 형태의 특징을 찾는 데 그 목표를 둔다.

 

<데이터 마이닝의 기법>



 


독일신용평가 데이터 셋

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

> setwd('c:/Rwork')

> german<-read.table('germandata.txt')

> head(german,2) # 값들의 변수명이 없음.


> names<-c("check","duration","history","purpose","credit","savings","employment","installment",   "personal",   "debtors",    "residence",  "property",   "age", "others",       "housing",    "numcredits", "job", "residpeople",      "telephone",       "foreign"     ,"y")

> colnames(german)<-names

> head(german,2)

 

> german$y<-factor(german$y,levels=c(1,2),labels=c('good','bad'))

> head(german,2)

  

> summary(german)

#  residence,numcredits,residpeople 실제 범주형이지만 수치형으로 인식. 범주형으로 변환 필요

 

> class(german$residence) #integer 수치형

[1] "integer"

> class(german$check) #factor 범주형

[1] "factor"

> german$residence = factor(german$residence)

> german$numcredits = factor(german$numcredits)

> german$residpeople = factor(german$residpeople)

> class(german$residence) #factor 변환

[1] "factor"

> class(german$numcredits) #factor 변환

[1] "factor"

> class(german$residpeople) #factor 변환

[1] "factor"

> table(german$residence)

1   2   3   4

130 308 149 413 

 

출처: 데이터마이닝, 장영재

반응형
Posted by 마르띤
,
반응형

Stars plot 이용한 다변량 데이터의 시각화

> total<-read.table('성적.csv',header=T,sep=',')

> total

이름 국어 영어 수학 국사 화학 물리

1 박지영   90   85   55   88   91   79

2 김태함   70   65   80   75   76   89

3 김효섭   92   95   76   65   89   91

4 임경희   76   89   88   95  100   91

5 권혁진   97   87   83   91   86   91

6 하혜진   80   86   97   85   69   77

7 이준원   80   30   40   50   70   90

8 윤정웅   70   82   54   56   58   60

9 주시현   90   95  100   85   89   92

 

> row.names(total)<-total$이름 #학생 별로 성적을 불러오기 위한 작업

> total

이름 국어 영어 수학 국사 화학 물리

박지영 박지영   90   85   55   88   91   79

김태함 김태함   70   65   80   75   76   89

김효섭 김효섭   92   95   76   65   89   91

임경희 임경희   76   89   88   95  100   91

권혁진 권혁진   97   87   83   91   86   91

하혜진 하혜진   80   86   97   85   69   77

이준원 이준원   80   30   40   50   70   90

윤정웅 윤정웅   70   82   54   56   58   60

주시현 주시현   90   95  100   85   89   92

> total<-total[,2:7]

> total

국어 영어 수학 국사 화학 물리

박지영   90   85   55   88   91   79

김태함   70   65   80   75   76   89

김효섭   92   95   76   65   89   91

임경희   76   89   88   95  100   91

권혁진   97   87   83   91   86   91

하혜진   80   86   97   85   69   77

이준원   80   30   40   50   70   90

윤정웅   70   82   54   56   58   60

주시현   90   95  100   85   89   92

 

 

그래프 그리기

> stars(total,flip.labels = FALSE, draw.segments = FALSE, frameplot=TRUE, full=TRUE, main='성적 분석-Star Chart')

 


> stars(total,flip.labels = FALSE, draw.segments = TRUE, frameplot=TRUE, full=TRUE, main='성적 분석-Nightingale Chart')

 


> stars(total,flip.labels = FALSE, draw.segments = TRUE, frameplot=TRUE, full=FALSE, main='성적 분석-Nightingale Chart')

 


> stars(total,flip.labels = FALSE, draw.segments = TRUE, frameplot=TRUE, full=FALSE, main='성적 분석-Nightingale Chart',scale=T,key.loc=c(10,2))


반응형
Posted by 마르띤
,
반응형

1. 타이타닉 선실 등급 별 생존 여부

#Class=1st, 2nd, 3rd, Crew / Survived=Yes, No
> data("Titanic")
> str(Titanic)
 table [1:4, 1:2, 1:2, 1:2] 0 0 35 0 0 0 17 0 118 154 ...
 - attr(*, "dimnames")=List of 4
  ..$ Class   : chr [1:4] "1st" "2nd" "3rd" "Crew"
  ..$ Sex     : chr [1:2] "Male" "Female"
  ..$ Age     : chr [1:2] "Child" "Adult"
  ..$ Survived: chr [1:2] "No" "Yes"
> apply(Titanic,c(1,4),sum)
      Survived
Class   No Yes
  1st  122 203
  2nd  167 118
  3rd  528 178
  Crew 673 212
> addmargins(apply(Titanic,c(1,4),sum))
      Survived
Class    No Yes  Sum
  1st   122 203  325
  2nd   167 118  285
  3rd   528 178  706
  Crew  673 212  885
  Sum  1490 711 2201
> par(mfrow=c(1,2))
> mosaicplot(~Class+Survived,data=Titanic)
> mosaicplot(~Class+Survived,data=Titanic,color=c('grey','red'))


2. 타이타닉 성별/연령 별 생존 여부

> mosaicplot(~Sex+Survived,data=Titanic,color=c('grey','red'))
> mosaicplot(~Age+Survived,data=Titanic,color=c('grey','red'))


 3. 타이타닉 성인 남성 / 성인 여성 생존 여부

> mosaicplot(~Class+Survived,data=as.table(Titanic[,'Male','Adult',]),color=c('grey','red'), main="Male+Adult")
> mosaicplot(~Class+Survived,data=as.table(Titanic[,'Female','Adult',]),color=c('grey','red'), main="Female+Adult")


반응형

'KNOU > 1 데이터시각화' 카테고리의 다른 글

막대그래프  (0) 2016.03.08
원 그래프  (0) 2016.03.07
Posted by 마르띤
,

집합

KNOU/1 대학수학 2016. 3. 9. 09:56
반응형

1. 집합

> library(prob)
> U=c('a','b','c','d','e','f')
> A=c('a','b','c','d')
> B=c('a','d','e','f')
> C=c('a','b','c')
> union(A,B) #AUB
[1] "a" "b" "c" "d" "e" "f" > intersect(A,B) #A∩B [1] "a" "d" > intersect(union(A,B),C) #(AUB)∩C [1] "a" "b" "c" > setdiff(A,B) #A-B [1] "b" "c" > setdiff(union(A,B),C) #(AUB)-C [1] "d" "e" "f"


2. 벤 다이어그램

> library(VennDiagram)
> grid.newpage() #새페이지
> venn.plot<-draw.triple.venn(area1=65,area2=75,area3=85,n12=35,n23=15,n13=25,n123=5,
+ category=c('First','Second','Third'),
+ fill=c('blue','red','green'),lty='blank',
+ cex=2,cat.cex=2,cat.col=c('blue','red','green')) #cex=숫자크기, cat.cex=카테고리크기
> grid.draw(venn.plot)



반응형
Posted by 마르띤
,