티스토리 뷰
[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
결측치는 없다.
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
* 참고문헌 : 실리콘밸리 데이터 과학자가 알려주는 따라하며 배우는 데이터 과학(저자 : 권재명)
from http://kuklife.tistory.com/60 by ccl(A)