[Data Science] Winequality(red wine) 데이터 회귀 분석 - 와인 품질 예측

목표 : winequaliry-red 데이터를 통해 와인 품질을 예측하는 회귀 분석을 해보자.

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

* 데이터 다운로드 : http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv

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/wine-quality/winequality-red.csv" download.file(URL, destfile = "data.csv") data <- tbl_df(read.table(file.choose(), strip.white = TRUE,                           sep=";", header = TRUE)) glimpse(data) summary(data)

결측치는 없는 것으로 확인된다.

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) 반응변수 quality는 설명번수 alcohol과 가장 높은 상관관계를 보인다.

2) pH와 fixed.acidity 사이에서 가장 높은 상관관계를 보인다.

p1 <- data %>% ggplot(aes(quality)) + geom_bar() p2 <- data %>% ggplot(aes(factor(quality), alcohol)) + geom_boxplot() p3 <- data %>% ggplot(aes(factor(quality), volatile.acidity)) + geom_boxplot() p4 <- data %>% ggplot(aes(pH, fixed.acidity)) + geom_point(alpha=.1) + geom_smooth() grid.arrange(p1, p2, p3, p4, ncol=2)

1) 소수의 와인만이 품질이 좋고, 나쁨을 보인다.

2) 품질이 좋을수록 도수가 높은 편을 보인다.

3) 품질이 좋지 않을수록 volatile.acidity이 높은 모습이 보인다.

4) fixed.acidity와 pH는 음의 상관관계를 보인다.

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(quality ~ ., data=training) > summary(data_lm_full) Call: lm(formula = quality ~ ., data = training) Residuals:      Min       1Q   Median       3Q      Max  -2.61095 -0.37507 -0.03307  0.44355  1.96941  Coefficients:                        Estimate Std. Error t value Pr(>|t|)     (Intercept)           6.6046127 28.1081584   0.235  0.81428     fixed.acidity        -0.0179643  0.0344863  -0.521  0.60255     volatile.acidity     -0.9728225  0.1623529  -5.992 2.94e-09 *** citric.acid          -0.1412841  0.1894202  -0.746  0.45593     residual.sugar        0.0271264  0.0206856   1.311  0.19005     chlorides            -1.8275182  0.5769014  -3.168  0.00159 **  free.sulfur.dioxide   0.0036550  0.0027827   1.313  0.18935     total.sulfur.dioxide -0.0042632  0.0009551  -4.463 9.04e-06 *** density              -1.5722521 28.6800242  -0.055  0.95629     pH                   -0.6948722  0.2479114  -2.803  0.00517 **  sulphates             1.0519115  0.1489324   7.063 3.15e-12 *** alcohol               0.3029134  0.0355065   8.531  < 2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6554 on 947 degrees of freedom Multiple R-squared:  0.3703, Adjusted R-squared:  0.363  F-statistic: 50.62 on 11 and 947 DF,  p-value: < 2.2e-16 > predict(data_lm_full, newdata = data[1:5,])        1        2        3        4        5  5.027269 5.198621 5.257843 5.586480 5.027269

volatile.acidity, total.sulfur.dioxide 등의 변수가 유의한 것으로 나타났다.

다음은 이차상호작용 모형을 포함한 모형이다.

> data_lm_full2 = lm(quality ~ .^2, data = training) > summary(data_lm_full2) Call: lm(formula = quality ~ .^2, data = training) Residuals:      Min       1Q   Median       3Q      Max  -1.92756 -0.40104 -0.00544  0.40693  1.72930  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6336 on 892 degrees of freedom Multiple R-squared:  0.4457, Adjusted R-squared:  0.4047  F-statistic: 10.87 on 66 and 892 DF,  p-value: < 2.2e-16 > length(coef(data_lm_full2)) [1] 67

이차상호작용 모형은 모수는 67개이다.

다음은 stepwise 변수 선택 모형이다.

> data_step = stepAIC(data_lm_full, scope = list(upper = ~ .^2, lower = ~1)) > summary(data_step) Call: lm(formula = quality ~ fixed.acidity + volatile.acidity + citric.acid +      residual.sugar + chlorides + free.sulfur.dioxide + total.sulfur.dioxide +      density + pH + sulphates + alcohol + volatile.acidity:total.sulfur.dioxide +      residual.sugar:alcohol + chlorides:free.sulfur.dioxide +      fixed.acidity:chlorides + total.sulfur.dioxide:density +      citric.acid:total.sulfur.dioxide + citric.acid:pH + fixed.acidity:pH +      fixed.acidity:residual.sugar + density:sulphates + chlorides:density +      free.sulfur.dioxide:density + free.sulfur.dioxide:sulphates +      total.sulfur.dioxide:pH, data = training) Residuals:      Min       1Q   Median       3Q      Max  -2.10955 -0.37468 -0.02552  0.40886  1.72726  Coefficients:  ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6343 on 933 degrees of freedom Multiple R-squared:  0.4189, Adjusted R-squared:  0.4033  F-statistic:  26.9 on 25 and 933 DF,  p-value: < 2.2e-16 > length(coef(data_step)) [1] 26

stepwise 변수 선택의 모수는 26개이며, 많이 줄어든 모습을 볼 수 있다.

모든 모형에 대한 RMSE 값을 구해 평가해보면 다음과 같다.

> y_obs <- validation$quality > 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] 0.6471014 > rmse(y_obs, yhat_lm_2) [1] 0.663733 > rmse(y_obs, yhat_step) [1] 0.6441662

stepwise 변수 선택 모형의 RMSE 값이 제일 낮은걸로 보아, 예측력이 가장 좋은 것을 알 수 있다.

6. 라쏘 모형 적합

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

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

> predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,])             1 1172 5.733698 825  5.659219 1178 6.353398 684  5.684794 1513 5.474132 > y_obs <- validation$quality > 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] 0.6283625

stepwise모형보다 RMSE값이 낮은 것으로 보아, glmnet 모형의 예측력이 더 좋음을 알 수 있다.

7. 나무 모형 적합

data_tr <- rpart(quality ~ ., data = training) data_tr 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] 0.6777646

glmnet 모형보다 많이 예측력이 떨어짐을 알 수 있다.

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

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

위에서 보았듯이, alcohol과의 상관관계가 가장 높음을 알 수 있다.

모형 평가는 다음과 같다.

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

glmnet보다 RMSE 값이 낮은 것으로 보아, 예측력이 더 좋음을 알 수 있다.

9. 부스팅 모형 적합

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

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

모형 평가는 다음과 같다.

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

RF보다 예측력이 떨어짐을 보인다.

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

> 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 0.6441662 2 glmnet 0.6283625 3     rf 0.5978746 4    gbm 0.6627304 > rmse(test$quality, predict(data_rf, newdata = test)) [1] 0.5998013

RF의 예측력이 가장 좋으며, 테스트세트 데이터를 RF 모형에 적합시켰을 때의 RMSE값은 0.599이다.

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

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')

gbm 모형의 예측오차율이 가장 적음을 보인다.

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

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)

이때까지 평가했던대로 RF 모형의 예측력이 가장 좋음을 보인다.

또한, lm과 glmnet 모형의 상관관계가 가장 높은 것도 알 수 있다.

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

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