[Data Science] Boston 데이터 회귀 분석 - 부동산 가격 예측 문제

목표 : Boston 데이터를 활용하여 Boston 지역의 부동산 가격 예측 회귀분석을 해보자.

1. 함수 작성 및 환경 준비

rmse <- function(yi, yhat_i){   sqrt(mean((yi - yhat_i)^2)) } binomial_deviance <- function(y_obs, yhat){   epsilon = 0.0001   yhat = ifelse(yhat < epsilon, epsilon, yhat)   yhat = ifelse(yhat > 1-epsilon, 1-epsilon, yhat)   a = ifelse(y_obs==0, 0, y_obs * log(y_obs/yhat))   b = ifelse(y_obs==1, 0, (1-y_obs) * log((1-y_obs)/(1-yhat)))   return(2*sum(a + b)) } panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...){   usr <- par("usr"); on.exit(par(usr))   par(usr = c(0, 1, 0, 1))   r <- abs(cor(x, y))   txt <- format(c(r, 0.123456789), digits = digits)[1]   txt <- paste0(prefix, txt)   if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)   text(0.5, 0.5, txt, cex = cex.cor * r) } library(tidyverse) library(gridExtra) library(MASS) library(glmnet) library(randomForest) library(gbm) library(rpart) library(boot) library(data.table) library(ROCR) library(ggplot2) library(dplyr) library(gridExtra)

2. 데이터 Read

Boston Data는 MASS Package 안에 있는 R 내장 Data set이다. 따라서, Read 후 간단하게 data 구성 요소정도만 파악하자.

data <- Boston glimpse(data) summary(data)

3. 시각화

set.seed(1810) pairs(data %>% dplyr::sample_n(min(1000, nrow(data))),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor)

산점도를 통해 반응변수 medv와 lstat(저소득층 주민 비율), rm(집 안의 방의 개수)의 변수와의 상관관계가 높음을 알 수 있다.

4. 훈련, 검증, 테스트세트 데이터 구분

set.seed(1810) n <- nrow(data) idx <- 1:n training_idx <- sample(idx, n * .60) idx <- setdiff(idx, training_idx) validate_idx <- sample(idx, n * .20) test_idx <- setdiff(idx, validate_idx) training <- data[training_idx,] validation <- data[validate_idx,] test <- data[test_idx,]

5. 선형 회귀 모형

선형회귀분석 명령은 lm()이다. glm(, family='gaussian')과 유사한 결과를 준다.

> data_lm_full <- lm(medv ~ ., data=training) > summary(data_lm_full) Call: lm(formula = medv ~ ., data = training) Residuals:      Min       1Q   Median       3Q      Max  -14.6480  -2.9439  -0.6139   2.0945  23.8238  Coefficients:               Estimate Std. Error t value Pr(>|t|)     (Intercept)  47.570421   7.159132   6.645 1.51e-10 *** crim         -0.093204   0.046597  -2.000 0.046415 *   zn            0.050774   0.018712   2.713 0.007058 **  indus         0.018522   0.083637   0.221 0.824892     chas          3.737056   1.088668   3.433 0.000685 *** nox         -20.042307   5.193933  -3.859 0.000141 *** rm            2.493443   0.577246   4.320 2.15e-05 *** age           0.021823   0.018260   1.195 0.233007     dis          -1.398149   0.279924  -4.995 1.02e-06 *** rad           0.353518   0.089955   3.930 0.000106 *** tax          -0.012438   0.005101  -2.438 0.015361 *   ptratio      -1.020976   0.184690  -5.528 7.24e-08 *** black         0.007527   0.003922   1.919 0.055953 .   lstat        -0.692811   0.070103  -9.883  < 2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 5.025 on 289 degrees of freedom Multiple R-squared:  0.7177, Adjusted R-squared:  0.705  F-statistic: 56.51 on 13 and 289 DF,  p-value: < 2.2e-16

예측력이 가장 강한 변수들이 어떤 것인지 쉽게 알 수 있다.

가장 유의한(가장 작은 p-값을 갖는 변수) 설명 변수들은 lstat, rm 이외에 dis(보스턴 5대 고용중심으로부터의 가중 평균 거리)와 ptratio가 있다.

lm 모형의 결과를 예측에 사용하려면 predict.lm() 함수를 사용하면 된다.

> predict(data_lm_full, newdata = data[1:5,])        1        2        3        4        5  30.43083 25.26532 30.29149 28.54173 27.45441

6. 선형회귀 모형에서 변수 선택

분류분석에선 간단한 선형 모형을 적합하였지만, 여기선 모든 이차상호작용을 고려한 모형을 사용할 것이다.

R에서는 formula interface를 사용하여 쉽게 이차상호작용 모형을 적합할 수 있다.

> data_lm_full2 = lm(medv ~ .^2, data = training) > summary(data_lm_full2) Call: lm(formula = medv ~ .^2, data = training) Residuals:     Min      1Q  Median      3Q     Max  -7.0332 -1.3980 -0.0998  1.2489 15.8191  ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.946 on 211 degrees of freedom Multiple R-squared:  0.9291, Adjusted R-squared:  0.8986  F-statistic:  30.4 on 91 and 211 DF,  p-value: < 2.2e-16

이차상호작용 모형은 모수의 개수가 92개에 달한다. 직접 확인해보면 다음과 같이 확인해볼 수 있다.

> length(coef(data_lm_full2)) [1] 92

위의 결과를 간단하게 설명하면 다음과 같다.

1) R^2 = 0.9291로, 수많은 변수를 사용한 이 모형은 y변수의 변동의 대부분을 설명해준다.

2) Adjusted R_2 = 0.8986으로 더 작고, F 통계량의 P-값은 0에 가깝다.

즉, 사용된 설명변수의 적어도 일부는 주택 가격을 예측하는데 도움이 된다는 것이다.

이 '지나치게 복잡한' 모형의 92개 모수 중에서 가장 중요한 모수를 선택하기 위해서 MASS::stepAIC()함수를 사용할 수 있다.

스텝(stepwise) 변수 선택으로 선형회귀 모형에서 중요한 변수를 자동으로 선택해준다.

다음 명령을 사용하여 변수 선택을 실행하자.

library(MASS) data_step = stepAIC(data_lm_full, scope = list(upper = ~ .^2, lower = ~1))  > summary(data_step) Call: lm(formula = medv ~ crim + zn + indus + chas + nox + rm + age +      dis + rad + tax + ptratio + black + lstat + rm:lstat + rad:lstat +      rm:rad + indus:dis + crim:chas + chas:rad + indus:lstat +      rm:ptratio + chas:nox + chas:rm + chas:lstat + crim:rm +      crim:lstat + crim:nox + age:rad + rm:age + age:lstat + nox:lstat +      tax:ptratio + indus:ptratio + chas:age + age:black + age:tax +      indus:tax + zn:rad + rm:black, data = training) Residuals:     Min      1Q  Median      3Q     Max  -8.7432 -1.3013 -0.1589  1.3760 18.4420  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.936 on 263 degrees of freedom Multiple R-squared:  0.9123, Adjusted R-squared:  0.8993  F-statistic: 70.14 on 39 and 263 DF,  p-value: < 2.2e-16 > length(coef(data_step)) [1] 40

최종 모형은 모수의 개수가 40개로 줄어들었음을 확인 할 수 있다.

즉, stepwise 변수 선택을 사용하자!!!

7. 모형 평가

y_obs <- validation$medv yhat_lm <- predict(data_lm_full, newdata=validation) yhat_lm_2 <- predict(data_lm_full2, newdata=validation) yhat_step <- predict(data_step, newdata=validation) > rmse(y_obs, yhat_lm) [1] 4.468361 > rmse(y_obs, yhat_lm_2) [1] 3.662617 > rmse(y_obs, yhat_step) [1] 3.409973

세 모형들 중에서 stepwise 변수 선택을 행한 모형이 예측 능력이 가장 높다.

8. 라쏘 모형 적합(glmnet)

xx <- model.matrix(medv ~ .^2-1, data) x <- xx[training_idx, ] y <- training$medv > glimpse(x)  num [1:303, 1:91] 6.539 0.54 9.232 0.198 6.393 ...  - attr(*, "dimnames")=List of 2   ..$ : chr [1:303] "371" "261" "372" "216" ...   ..$ : chr [1:91] "crim" "zn" "indus" "chas" ...

mode.matrix() 함수에서 ^2-1을 사용한 것은 모든 이차상호작용을 포함하기 위해서다. 앞의 lm() 결과에서 보듯 상호작용 후의 변수 선택 / 모형 선택이 더 좋은 예측력을 줄 것이기 때문이다.

data_cvfit <- cv.glmnet(x, y) plot(data_cvfit)

각 lambda 값에서 선택된 변수의 개수는 72개, 58개이다.

coef(data_cvfit, s = c("lambda.1se")) coef(data_cvfit, s = c("lambda.min")) > predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,])            1 371 52.02525 261 34.78688 372 32.72658 216 24.18534 476 11.56524  > y_obs <- validation$medv > yhat_glmnet <- predict(data_cvfit, s="lambda.min", newx=xx[validate_idx,]) > yhat_glmnet <- yhat_glmnet[,1] # change to a vector from [n*1] matrix > rmse(y_obs, yhat_glmnet) [1] 3.176944

 stepwise 변수 선택을 행한 모형보다 glmnet이 좋은 예측력을 가지는 것을 보인다.

9. 나무 모형

> data_tr <- rpart(medv ~ ., data = training) > data_tr n= 303   node), split, n, deviance, yval       * denotes terminal node  1) root 303 25845.6900 22.43795      2) lstat>=9.95 174  4204.2320 17.27184        4) lstat>=18.905 44   906.5298 12.10227          8) nox>=0.603 33   281.8388 10.30606 *        9) nox< 0.603 11   198.8091 17.49091 *      5) lstat< 18.905 130  1723.8400 19.02154         10) lstat>=14.395 63   602.6889 17.17778 *       11) lstat< 14.395 67   705.6057 20.75522 *    3) lstat< 9.95 129 10733.8400 29.40620        6) rm< 7.0115 95  3882.6520 25.78000         12) age< 89.45 87  1339.7340 24.57701           24) rm< 6.543 50   297.7282 22.39400 *         25) rm>=6.543 37   481.7330 27.52703 *       13) age>=89.45 8  1047.7990 38.86250 *      7) rm>=7.0115 34  2111.6200 39.53824         14) rm< 7.443 17   375.7894 34.50588 *       15) rm>=7.443 17   874.7953 44.57059 * opar <- par(mfrow = c(1,1), xpd = NA) plot(data_tr) text(data_tr, use.n = TRUE) par(opar)

1) 첫 번째 분할에서 왼쪽가지(lstat>=9.95)는 저소득 지역, 오른쪽 지역(lstat<9.95)는 고소득 지역

2) 왼쪽가지(저소득 지역)에서는 nox로 더 분할되고, 고소득 지역은 rm(방의 개수)로 더 분할된다.

즉, 나무 모형은 데이터의 다른 영역에서 변수들 간의 변화하는 상관관계를 찾아낼 수있다.

위와 마찬가지로 나무 모형도 RMSE를 통해 검증해보면 다음과 같다.

> yhat_tr <- predict(data_tr, validation) > rmse(y_obs, yhat_tr) [1] 3.71948

stepwise 변수 선택을 통한 모형과 glmnet 모형보다 떨어지는 예측력을 보인다.

10. 랜덤 포레스트(RF)

set.seed(1810) data_rf <- randomForest(medv ~ ., training) > data_rf Call:  randomForest(formula = medv ~ ., data = training)                 Type of random forest: regression                      Number of trees: 500 No. of variables tried at each split: 4           Mean of squared residuals: 13.27373                     % Var explained: 84.44 opar <- par(mfrow=c(1,2)) plot(data_rf) varImpPlot(data_rf) par(opar)

RF에서도 알 수 있듯, rm과 lstat 변수가 높은 상관관계를 보인다. RMSE를 계산해보면 다음과 같다.

> yhat_rf <- predict(data_rf, newdata=validation) > rmse(y_obs, yhat_rf) [1] 2.943556

이때까지 살펴본 모형들중 RF의 RMSE 값이 가장 낮다. 즉, 가장 예측력이 좋다.

11. 부스팅(gbm)

set.seed(1810) data_gbm <- gbm(medv ~ ., data=training,                 n.trees=30000, cv.folds=3, verbose = TRUE) > (best_iter = gbm.perf(data_gbm, method="cv")) [1] 586

최적의 트리 개수는 586개이다. gbm의 RMSE 값을 구해보면 다음과 같다.

> yhat_gbm <- predict(data_gbm, n.trees=best_iter, newdata=validation) > rmse(y_obs, yhat_gbm) [1] 3.598836

예측력이 RF보다 떨어지는 것이 보인다.

12. 최종 모형 선택과 테스트세트 오차 계산

> data.frame(lm = rmse(y_obs, yhat_step), +            glmnet = rmse(y_obs, yhat_glmnet), +            rf = rmse(y_obs, yhat_rf), +            gbm = rmse(y_obs, yhat_gbm)) %>% +   reshape2::melt(value.name = 'rmse', variable.name = 'method') No id variables; using all as measure variables   method     rmse 1     lm 3.409973 2 glmnet 3.176944 3     rf 2.943556 4    gbm 3.598836

RF > glmnet > stepwise 변수 선택 모형 > gbm 순으로 예측력이 좋으므로 RF를 선택한다.

RF의 테스트세트에서의 오차는 다음과 같다.

> rmse(test$medv, predict(data_rf, newdata = test)) [1] 2.832759

여러 회귀분석 방법의 오차를 비교하는 시각화 방법 중 하나는 예측오차의 분포를 병렬상자그림으로 보여주는 것이다.

boxplot(list(lm = y_obs-yhat_step,              glmnet = y_obs-yhat_glmnet,              rf = y_obs-yhat_rf,              gbm = y_obs-yhat_gbm), ylab="Error in Validation Set") abline(h=0, lty=2, col='blue')

또한 모형들간의 예측값들끼리 산점도 행렬을 그려보면 다음과 같다.

pairs(data.frame(y_obs=y_obs,                  yhat_lm=yhat_step,                  yhat_glmnet=c(yhat_glmnet),                  yhat_rf=yhat_rf,                  yhat_gbm=yhat_gbm),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor)

그림에서 알 수 있듯 stepwise 변수 선택 모형과 glmnet 예측값끼리는 상관계수가 높고, RF와 gbm값끼리도 상관계수가 높은 것으로 보인다.

관측값과 상관계수가 가장 높은 방법은 RF이다.

* 전체 소스코드 : https://github.com/SeongkukCHO/PCOY/blob/master/data-science/Boston/Boston.R

* 참고문헌 : 실리콘밸리 데이터 과학자가 알려주는 따라하며 배우는 데이터 과학(저자 : 권재명)