반응형

보스턴 하우징 데이터 랜덤포레스트 방법의 회귀앙상블 모형

1) 데이터 읽기

> Boston$chas=factor(Boston$chas)

> Boston$rad=factor(Boston$rad)

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

Factor 함수를 이용하여 범주형으로 변경하고, medv 변수를 목표변수로 다른 변수를 입력변수로 사용한다.

 

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

>library(randomForest)

> rf.boston<-randomForest(medv~.,data=Boston,ntree=100,mtry=5,importance=T,na.action=na.omit)

> rf.boston

 

Call:

  randomForest(formula = medv ~ ., data = Boston, ntree = 100,      mtry = 5, importance = T, na.action = na.omit)

Type of random forest: regression

Number of trees: 100

No. of variables tried at each split: 5

 

Mean of squared residuals: 9.743395

% Var explained: 88.46

 

> summary(rf.boston)

Length Class  Mode    

call              7    -none- call    

type              1    -none- character

predicted       506    -none- numeric 

mse             100    -none- numeric 

rsq             100    -none- numeric 

oob.times       506    -none- numeric 

importance       26    -none- numeric 

importanceSD     13    -none- numeric 

localImportance   0    -none- NULL    

proximity         0    -none- NULL    

ntree             1    -none- numeric 

mtry              1    -none- numeric 

forest           11    -none- list    

coefs             0    -none- NULL    

y               506    -none- numeric 

test              0    -none- NULL    

inbag             0    -none- NULL    

terms             3    terms  call

 함수 설명

ntree=100, 분류기 개수. 디폴트는 500

mtry=5, 중간노드마다 랜덤하게 선택되는 변수들의 개수. 디폴트는 분류나무의 경우 sqrt(p), 회귀나무의 경우 p/3

importance=T, 변수의 중요도 계산, 디폴트는 F

na.action=na.omit, 결측치를 처리하는 방법, 변수의 중요도를 계산하게 하고, 결측치는 필요한 경우에만 삭제.

 

 

names()함수를 통해 rf.boston에 저장된 오브젝트의 리스트를 불러내어, $predicted를 이용하여 훈련 데이터의 예측 집단을 출력할 수 있다.

> names(rf.boston)

[1] "call"            "type"            "predicted"       "mse"             "rsq"           

[6] "oob.times"       "importance"      "importanceSD"    "localImportance" "proximity"     

[11] "ntree"           "mtry"            "forest"          "coefs"           "y"              

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

> head(rf.boston$predicted,30)

1        2        3        4        5        6        7        8        9       10

28.25382 22.55963 34.14192 35.45333 34.06798 26.81151 21.01950 16.78839 17.80599 19.23591

11       12       13       14       15       16       17       18       19       20

21.02440 21.23466 21.80889 20.05162 19.30557 20.21721 21.61349 18.46000 18.14724 19.96174

21       22       23       24       25       26       27       28       29       30

14.10136 18.55984 16.05801 15.04825 16.70996 15.70548 17.84748 14.82048 18.88633 20.64939

 

importance()함수를 통해 계산된 입력변수의 중요도를 알 수 있다.

> importance(rf.boston,type=1)

%IncMSE

crim     8.325232

zn       2.061869

indus    5.130483

chas     1.030915

nox      9.211906

rm      17.090802

age      5.229782

dis      8.322716

rad      5.342500

tax      4.604745

ptratio  7.102056

black    5.292651

lstat   14.652271

결과를 보면 rm변수과 lstat 변수의 중요도가 가장 높음을 알 수 있다.

함수 설명

 type=1,은 정분류율의 평균감소값, 2는 불순도의 평균감소값을 이용하여 계산

 

목표변수의 적합값을 구하고 평가하기 위해 평균오차제곱합(mse)를 계산.

> names(Boston)

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

[10] "tax"     "ptratio" "black"   "lstat"   "medv"  

> Boston$medv.hat = predict(rf.boston,newdata=Boston)

> mean((Boston$medv-Boston$medv.hat)^2) #mean square error(mse)

[1] 1.915207

 

기존 선형회귀 한 회귀 분류 나무 모형 결과의 평균오차제곱합 mean((Boston$medv-Boston$medv.hat)^2) = 10.8643  대비 랜덤포레스트의 평균오차제곱합이 1.915207로 설명력이 매우 증가되었음을 알 수 있다. 랜덤포레스트의 경우 부트스트랩을 이용하기 때문에 확률임의추출에 의한 변동성이 있을 수 있다. 따라서 모델링을 할 때 마다 결과가 다르기 때문에, 랜덤포레스트를 수차례 반복 시행하고 예측결과의 평균값을 취하는 경우도 있다.

(기존 보스턴 하우징 데이터 회귀나무모 사례 분석 링크 로가기)

 

랜덤 포레스트 회귀앙상블의 적합값과 실제값의 일치도를 보자. 예측일치도가 우수함을 알 수 있다.

> plot(Boston$medv,Boston$medv.hat,xlab='Observed Values',ylab='Fitted Values')

> abline(0,1)


 

 

 

기존 분류 회귀의 나무모형의 적합값과 실제값의 일치도와 비교해봐도 매우 우수함을 알 수 있다.

 

 

 

 

이 의사결정 나무를 활용하여 30% 검증 데이터에 적용시켜서 분류예측치를 구해보자. 그리고 그 예측치를 구해보자.

> set.seed(1234)

> nrow(Boston)

[1] 506

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

> rf.train.boston<-randomForest(medv~.,data=Boston.train,ntree=100,importance=T,na.action=na.omit)

#obtain the predicted values

> medv.hat.test<-predict(rf.train.boston,newdata=Boston.test)

> mean((Boston.test$medv-medv.hat.test)^2) #predicted mean square error

[1] 4.114596

검증 데이터에 대한 평균오차제곱합 mse 4.11로 계산되었다. 기존 회귀 분류 나무모형 의 검증 데이터 오분류율 13.95258와 비교해보면 상당히 향상된 결과임을 알 수 있다.

(기존 보스턴 하우징 데이터 회귀나무모 사례 분석 링크 로가기)


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

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

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

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


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

목표변수가 연속형인 경우 -> 선형 회귀모델, 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 마르띤
,