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

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

* 데이터 설명 : 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-white.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-white.csv" download.file(URL, destfile = "data4.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) density와 alcohol, density와 residual.sugar의 상관관계가 높음을 알 수 있다.

조금 더 세부적으로 확인해보면 다음과 같다.

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

1) 소수의 와인만이 아주 좋거나 아주 나쁨을 알 수 있다.

2) alcohol의 도수가 높을수록 좋은 품질을 가지는 영향이 있다.

3) 품질의 값이 6을 가지는 와인에서 이상치를 확인 할 수 있다.

4) alcohol의 도수가 높을수록 낮은 밀도를 보이는 음의 상관관계를 가진다.

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  -3.6200 -0.4998 -0.0518  0.4659  3.0616  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7503 on 2926 degrees of freedom Multiple R-squared:  0.2898, Adjusted R-squared:  0.2871  F-statistic: 108.5 on 11 and 2926 DF,  p-value: < 2.2e-16 > predict(data_lm_full, newdata = data[1:5,])        1        2        3        4        5  5.504005 5.200614 5.833953 5.798921 5.798921

citric.acid, chlorides, 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  -3.2777 -0.4821 -0.0286  0.4414  3.0779  Coefficients: Residuals:      Min       1Q   Median       3Q      Max  -3.13596 -0.48543 -0.03981  0.44439  3.10824  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7149 on 2871 degrees of freedom Multiple R-squared:  0.3674, Adjusted R-squared:  0.3529  F-statistic: 25.26 on 66 and 2871 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 + free.sulfur.dioxide:total.sulfur.dioxide +      volatile.acidity:alcohol + residual.sugar:pH + fixed.acidity:citric.acid +      chlorides:alcohol + free.sulfur.dioxide:alcohol + fixed.acidity:free.sulfur.dioxide +      pH:sulphates + total.sulfur.dioxide:sulphates + free.sulfur.dioxide:sulphates +      total.sulfur.dioxide:alcohol + residual.sugar:chlorides +      chlorides:pH + volatile.acidity:pH + citric.acid:pH + volatile.acidity:total.sulfur.dioxide +      sulphates:alcohol + residual.sugar:density + residual.sugar:alcohol +      chlorides:sulphates, data = training) Residuals:      Min       1Q   Median       3Q      Max  -3.13596 -0.48543 -0.03981  0.44439  3.10824  Coefficients: ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7137 on 2906 degrees of freedom Multiple R-squared:  0.3618, Adjusted R-squared:  0.355  F-statistic: 53.15 on 31 and 2906 DF,  p-value: < 2.2e-16 > length(coef(data_step)) [1] 32

stepwise 변수 선택 모형의 모수는 32개로써, 이전보다 많이 줄어든 모습이 보인다.

모든 모형에 대해 평가를 하면 다음과 같다.

> 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.7647687 > rmse(y_obs, yhat_lm_2) [1] 1.00542 > rmse(y_obs, yhat_step) [1] 0.8777487

lm() 모형의 RMSE값이 가장 낮은걸로 보면 예측력이 가장 높은 것을 알 수 있다.

6. 라쏘 모형 적합

> xx <- model.matrix(quality ~ .^2-1, data) > x <- xx[training_idx, ] > y <- training$quality > glimpse(x)  num [1:2938, 1:66] 6.8 7.3 7 6.7 6 5.7 6 6.9 6.6 6.6 ...  - attr(*, "dimnames")=List of 2   ..$ : chr [1:2938] "3590" "2528" "3610" "2098" ...   ..$ : chr [1:66] "fixed.acidity" "volatile.acidity" "citric.acid" "residual.sugar" ... > data_cvfit <- cv.glmnet(x, y) > plot(data_cvfit)

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

> predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,])             1 3590 5.863399 2528 5.906761 3610 6.012621 2098 5.625118 4641 6.324870 > 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.7616714

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

7. 나무 모형 적합

> data_tr <- rpart(quality ~ ., data = training) > data_tr n= 2938  node), split, n, deviance, yval       * denotes terminal node  1) root 2938 2319.62000 5.878489      2) alcohol< 10.91667 1896 1163.61800 5.608650        4) volatile.acidity>=0.2875 682  292.62320 5.277126 *      5) volatile.acidity< 0.2875 1214  753.92830 5.794893         10) volatile.acidity>=0.2075 761  419.32190 5.649146 *       11) volatile.acidity< 0.2075 453  291.28480 6.039735 *    3) alcohol>=10.91667 1042  766.74950 6.369482        6) free.sulfur.dioxide< 11.5 59   54.64407 5.542373 *      7) free.sulfur.dioxide>=11.5 983  669.32040 6.419125         14) alcohol< 11.775 456  273.20830 6.208333 *       15) alcohol>=11.775 527  358.31880 6.601518 * > 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.7677215

나무 모형의 예측력은 lm과 glmnet 두 모형보다 다 떨어지는 모습을 보인다.

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

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

위에서 봤듯이, alcohol의 변수와의 상관관계가 높음을 보인다.

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

이때까지 적합된 모형들의 예측력에 비해 RF의 예측력이 더 뛰어남을 보인다.

9. 부스팅 모형 적합

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

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

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

RMSE값이 RF보다 높은 것으로 보아, RF보다 예측력이 떨어짐을 알 수 있다.

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

> data.frame(lm = rmse(y_obs, yhat_lm), +            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.7647687 2 glmnet 0.7616714 3     rf 0.6427712 4    gbm 0.7269684

RF 모형의 예측력이 가장 좋다. 따라서 테스트세트를 RF모형으로 적합시켜보자.

> rmse(test$quality, predict(data_rf, newdata = test)) [1] 0.5979055

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

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)

1) 가장 좋은 예측력을 가지는 모형은 RF인것을 알 수 있다.

2) lm과 glmnet, rf와 gbm 끼리의 상관계수가 높음을 보인다.

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

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