[Data Science] abalone 데이터 회귀 분석 - 전복 나이 예측

목표 : abalone 데이터를 통해 전복 나이를 예측하는 회귀 분석을 해보자.

* 데이터 설명 : http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.names

* 데이터 다운로드 : http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data

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

URL <- "http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data" download.file(URL, destfile = "data.csv") data <- tbl_df(read.table(file.choose(), strip.white = TRUE,                           sep=",", header = FALSE)) names(data) <- c('Sex','Length','Diameter','Height','Whole',                  'Shucked','Viscera','Shell','Rings') > glimpse(data) Observations: 4,177 Variables: 9 $ Sex      <fct> M, M, F, M, I, I, F, F, M, F, F, M, M, F, F, M, I, F, M, M, M, I, F, F, ... $ Length   <dbl> 0.455, 0.350, 0.530, 0.440, 0.330, 0.425, 0.530, 0.545, 0.475, 0.550, 0.... $ Diameter <dbl> 0.365, 0.265, 0.420, 0.365, 0.255, 0.300, 0.415, 0.425, 0.370, 0.440, 0.... $ Height   <dbl> 0.095, 0.090, 0.135, 0.125, 0.080, 0.095, 0.150, 0.125, 0.125, 0.150, 0.... $ Whole    <dbl> 0.5140, 0.2255, 0.6770, 0.5160, 0.2050, 0.3515, 0.7775, 0.7680, 0.5095, ... $ Shucked  <dbl> 0.2245, 0.0995, 0.2565, 0.2155, 0.0895, 0.1410, 0.2370, 0.2940, 0.2165, ... $ Viscera  <dbl> 0.1010, 0.0485, 0.1415, 0.1140, 0.0395, 0.0775, 0.1415, 0.1495, 0.1125, ... $ Shell    <dbl> 0.150, 0.070, 0.210, 0.155, 0.055, 0.120, 0.330, 0.260, 0.165, 0.320, 0.... $ Rings    <int> 15, 7, 9, 10, 7, 8, 20, 16, 9, 19, 14, 10, 11, 10, 10, 12, 7, 10, 7, 9, ...

결측치는 없다.

3. 시각화

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

1) 반응변수 RIngs에 대해 설명변수 Shell이 상관관계가 가장 높음을 알 수 있다.

2) Whole과 Shucked, Whole과 Shell의 상관관계가 가장 높음을 알 수 있다.

조금 더 자세히 살펴보면 다음과 같다.

p1 <- data %>% ggplot(aes(Rings)) + geom_bar() p2 <- data %>% ggplot(aes(factor(Rings), Shell)) + geom_boxplot() p3 <- data %>% ggplot(aes(factor(Rings), Height)) + geom_boxplot() p4 <- data %>% ggplot(aes(Shell, Whole)) + geom_point(alpha=.1) + geom_smooth() grid.arrange(p1, p2, p3, p4, ncol=2)

1) 전복이 5~15살에 많이 분포되어 있음을 알 수 있다.

2) 나이가 많을수록 말린 후의 무게가 많이 나감을 알 수 있다.

3) 무게와 나이가 양의 상관관계를 보임을 알 수 있다.

4) 당연한거겠지만 전체 무게와 말린 후의 무게가 강한 양의 상관관계를 보인다.

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 모형을 적합하면 다음과 같다.

> data_lm_full <- lm(Rings ~ ., data=training) > summary(data_lm_full) Call: lm(formula = Rings ~ ., data = training) Residuals:     Min      1Q  Median      3Q     Max  -8.3557 -1.3089 -0.3421  0.8809 14.2300  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.166 on 2496 degrees of freedom Multiple R-squared:  0.5411, Adjusted R-squared:  0.5394  F-statistic:   327 on 9 and 2496 DF,  p-value: < 2.2e-16 > predict(data_lm_full, newdata = data[1:5,])         1         2         3         4         5   8.904290  7.728173 10.936506  9.627705  6.644735

SexM과 Length를 제외한 모든 변수가 유의함을 보인다.

이차상호작용 모형을 적합하면 다음과 같다.

> data_lm_full2 = lm(Rings ~ .^2, data = training) > summary(data_lm_full2) Call: lm(formula = Rings ~ .^2, data = training) Residuals:      Min       1Q   Median       3Q      Max  -10.6899  -1.2168  -0.2542   0.8566  13.8243  ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.086 on 2461 degrees of freedom Multiple R-squared:  0.5805, Adjusted R-squared:  0.573  F-statistic:  77.4 on 44 and 2461 DF,  p-value: < 2.2e-16 > length(coef(data_lm_full2)) [1] 45

이차상호작용의 모수는 45개이다.

stepwise 변수 선택 모형을 적합하면 다음과 같다.

> data_step = stepAIC(data_lm_full, scope = list(upper = ~ .^2, lower = ~1)) > summary(data_step) Call: lm(formula = Rings ~ Sex + Length + Diameter + Height + Whole +      Shucked + Viscera + Shell + Sex:Shucked + Shucked:Viscera +      Sex:Diameter + Sex:Shell + Height:Shell + Height:Shucked +      Length:Viscera + Length:Whole + Length:Diameter + Diameter:Viscera +      Whole:Viscera + Shucked:Shell + Diameter:Shell, data = training) Residuals:      Min       1Q   Median       3Q      Max  -10.6159  -1.2343  -0.2731   0.8435  14.2758  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.084 on 2480 degrees of freedom Multiple R-squared:  0.578, Adjusted R-squared:  0.5737  F-statistic: 135.9 on 25 and 2480 DF,  p-value: < 2.2e-16 > length(coef(data_step)) [1] 26

stepwise 변수 선택 모형의 모수는 26개이다.

세 가지의 모형 평가를 실시해보면 다음과 같다.

> y_obs <- validation$Rings > 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] 2.169579 > rmse(y_obs, yhat_lm_2) [1] 2.133406 > rmse(y_obs, yhat_step) [1] 2.136988

이차상호작용 모형을 적용했을 때, RMSE값이 가장 낮은 것으로 보았을 때 예측력이 가장 좋음을 알 수 있다.

6. 라쏘 모형 적합

xx <- model.matrix(Rings ~ .^2-1, data) x <- xx[training_idx, ] y <- training$Rings glimpse(x) data_cvfit <- cv.glmnet(x, y) plot(data_cvfit)

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

모형 평가는 다음과 같다.

> y_obs <- validation$Rings > 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] 2.118328

이차상호작용 모형보다 glmnet 모형의 RMSE값이 더 작은 것으로 보았을 때, 예측력이 더 좋음을 알 수 있다.

7. 나무 모형 적합

> data_tr <- rpart(Rings ~ ., data = training) > data_tr n= 2506  node), split, n, deviance, yval       * denotes terminal node  1) root 2506 25519.5100  9.882682      2) Shell< 0.15375 765  3131.1950  7.308497        4) Diameter< 0.2225 121   180.3306  4.925620 *      5) Diameter>=0.2225 644  2134.7250  7.756211         10) Sex=I 415   872.8048  7.178313 *       11) Sex=F,M 229   872.1572  8.803493 *    3) Shell>=0.15375 1741 15091.6700 11.013790        6) Shell< 0.4095 1470 10021.9900 10.602720         12) Shell< 0.24925 594  2772.3520  9.796296 *       13) Shell>=0.24925 876  6601.4100 11.149540           26) Shucked>=0.42975 574  2172.0350 10.449480 *         27) Shucked< 0.42975 302  3613.3810 12.480130             54) Whole< 0.98875 231  2246.2940 11.891770 *           55) Whole>=0.98875 71  1026.9580 14.394370 *      7) Shell>=0.4095 271  3473.9260 13.243540         14) Shucked>=0.58525 194  2037.9850 12.448450           28) Shell< 0.568 160  1150.9750 11.862500 *         29) Shell>=0.568 34   573.5588 15.205880 *       15) Shucked< 0.58525 77  1004.3120 15.246750 * > opar <- par(mfrow = c(1,1), xpd = NA) > plot(data_tr) > text(data_tr, use.n = TRUE) > par(opar)

모형 평가는 다음과 같다.

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

다른 모형들에 비해 높은 RMSE값을 보이며, 그에 따라 다른 모형에 비해 좋지 않은 모형임을 알 수 있다.

8. 랜덤 포레스트 모형 적합

set.seed(1810) data_rf <- randomForest(Rings ~ ., training) data_rf opar <- par(mfrow=c(1,2)) plot(data_rf) varImpPlot(data_rf) par(opar)

위에서도 알 수 있다시피, Shell변수와의 상관관계가 높은 것을 다시 한번 알 수 있다.

모형 평가는 다음과 같다.

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

현재까지 적합된 모형들과 비교했을 때 가장 좋은 예측력을 보인다.

9. 부스팅 모형 적합

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

최적의 트리 개수는 659개이다.

모형 평가는 다음과 같다.

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

RF모형 보다 떨어지는 예측력을 보인다.

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

> data.frame(lm = rmse(y_obs, yhat_lm_2), +            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 2.133406 2 glmnet 2.118328 3     rf 2.103759 4    gbm 2.178855 > rmse(test$Rings, predict(data_rf, newdata = test)) [1] 2.225618

RF 모형을 적합할 때 가장 낮은 RMSE 값을 나타낸다. 즉, 가장 좋은 예측 모형이다.

따라서, RF모형을 통해 테스트세트 데이터를 적합시키면 RMSE=2.225를 얻을 수 있다.

예측오차의 분포를 병렬상자그림으로 오차를 비교해보면 다음과 같다.

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)

glmnet 모형과 rf모형의 예측력이 가장 높음을 알 수 있으며, glmnet과 lm모형이 높은 상관관계를 보임을 알 수 있다.

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

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