[Data Science] spambase 데이터 분류 분석 - 스펨 메일 예측 문제

목표 : 스펨베이스 데이터를 활용하여 스펨 메일 분류 분석을 해보자.

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

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

데이터 기본 구성

word_freq_make:         continuous.

word_freq_address:      continuous.

word_freq_all:          continuous.

word_freq_3d:           continuous.

word_freq_our:          continuous.

word_freq_over:         continuous.

word_freq_remove:       continuous.

word_freq_internet:     continuous.

word_freq_order:        continuous.

word_freq_mail:         continuous.

word_freq_receive:      continuous.

word_freq_will:         continuous.

word_freq_people:       continuous.

word_freq_report:       continuous.

word_freq_addresses:    continuous.

word_freq_free:         continuous.

word_freq_business:     continuous.

word_freq_email:        continuous.

word_freq_you:          continuous.

word_freq_credit:       continuous.

word_freq_your:         continuous.

word_freq_font:         continuous.

word_freq_000:          continuous.

word_freq_money:        continuous.

word_freq_hp:           continuous.

word_freq_hpl:          continuous.

word_freq_george:       continuous.

word_freq_650:          continuous.

word_freq_lab:          continuous.

word_freq_labs:         continuous.

word_freq_telnet:       continuous.

word_freq_857:          continuous.

word_freq_data:         continuous.

word_freq_415:          continuous.

word_freq_85:           continuous.

word_freq_technology:   continuous.

word_freq_1999:         continuous.

word_freq_parts:        continuous.

word_freq_pm:           continuous.

word_freq_direct:       continuous.

word_freq_cs:           continuous.

word_freq_meeting:      continuous.

word_freq_original:     continuous.

word_freq_project:      continuous.

word_freq_re:           continuous.

word_freq_edu:          continuous.

word_freq_table:        continuous.

word_freq_conference:   continuous.

char_freq_;:            continuous.

char_freq_(:            continuous.

char_freq_[:            continuous.

char_freq_!:            continuous.

char_freq_$:            continuous.

char_freq_#:            continuous.

capital_run_length_average: continuous.

capital_run_length_longest: continuous.

capital_run_length_total:   continuous.

1. 함수 작성 및 환경 준비

rm(list=ls()) 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)

2. 데이터 다운로드 및 Read

URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.data" download.file(URL, destfile = "data3.csv") data <- tbl_df(read.table(file.choose(), strip.white = TRUE,                           sep=",", header = FALSE)) glimpse(data) names(data) <- c('word_freq_make', 'word_freq_address', 'word_freq_all', 'word_freq_3d', 'word_freq_our',                  'word_freq_over', 'word_freq_remove', 'word_freq_internet', 'word_freq_order', 'word_freq_mail',                  'word_freq_receive', 'word_freq_will', 'word_freq_people', 'word_freq_report', 'word_freq_addresses',                  'word_freq_free', 'word_freq_business', 'word_freq_email', 'word_freq_you', 'word_freq_credit',                  'word_freq_your', 'word_freq_font', 'word_freq_000', 'word_freq_money', 'word_freq_hp',                  'word_freq_hpl', 'word_freq_george', 'word_freq_650', 'word_freq_lab', 'word_freq_labs',                  'word_freq_telnet', 'word_freq_857', 'word_freq_data', 'word_freq_415', 'word_freq_85',                  'word_freq_technology', 'word_freq_1999', 'word_freq_parts', 'word_freq_pm', 'word_freq_direct',                  'word_freq_cs', 'word_freq_meeting', 'word_freq_original', 'word_freq_project', 'word_freq_re',                  'word_freq_edu', 'word_freq_table', 'word_freq_conference', 'char_freq_;', 'char_freq_(',                  'char_freq_[', 'char_freq_!', 'char_freq_$', 'char_freq_#', 'capital_run_length_average',                  'capital_run_length_longest', 'capital_run_length_total',                  # 'spam'                  'class' names(data)[58] <- 'class' data$class <- factor(data$class) glimpse(data)

결측치는 없다.

3. 데이터 시각화

변수가 너무 많으므로 처음 10개와 class의 관계, 마지막 10개와 class의 관계만 살펴보자.

set.seed(1810) #1~10번째의 설명변수와 class간의 관계 pairs(data %>% dplyr::select(1:10, 58) %>%         sample_n(min(1000, nrow(data))),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor) #48~57번째의 설명변수와 class간의 관계 set.seed(1810) pairs(data %>% dplyr::select(48:57, 58) %>%         sample_n(min(1000, nrow(data))),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor)

1) 첫 번째 산점도

make, address, all, 3d, our 등의 처음 10개의 단어의 발생빈도퍼센트와 반응변수 class와의 상관관계를 보여준다.

이들 중 make 단어에 해당하는 word_freq_make 변수를 보면 본문에 make란 단어가 많이 등장하면 스팸일 가능성이 높다.

2) 두 번째 산점도

char_freq_$ 변수가 스팸인지 아닌지와 상관관계가 가장 높다. 따라서 본문에 $ 표시가 많이 나타나면 스팸일 가능성이 높다.

그렇다면 57개의 예측변수 중 반응변수 class와 가장 상관관계가 높은 예측변수는 어떤 것일까?

직접 상관관계를 계산한 후, 변수를 상간 관계 값의 크기로 정렬하여 시각화를 해보자.

tmp <- as.data.frame(cor(data[,-58], as.numeric(data$class))) tmp <- tmp %>% rename(cor=V1) tmp$var <- rownames(tmp) tmp %>%   ggplot(aes(reorder(var, cor), cor)) +   geom_point() +   coord_flip()

그림을 보면 your, 000, remove, $,you, free, business 등의 단어와 글자들이 스팸 여부와 가장 상관관계가 높음을 알 수 있다.

다음 코드로 스팸 여부와 가장 상관관계가 높은 몇 가지 변수를 시각화해보자.

library(ggplot2) library(dplyr) library(gridExtra) glimpse(data) p1 <- data %>% ggplot(aes(class)) + geom_bar() p2 <- data %>% ggplot(aes(class, `char_freq_$`)) +   geom_jitter(col='gray') +   geom_boxplot(alpha=.5) +   scale_y_sqrt() p3 <- data %>% ggplot(aes(`char_freq_$`, group=class, fill=class)) +   geom_density(alpha=.5) +   scale_x_sqrt() + scale_y_sqrt() p4 <- data %>% ggplot(aes(class, capital_run_length_longest)) +   geom_jitter(col='gray') +   geom_boxplot(alpha=.5) +   scale_y_log10() grid.arrange(p1, p2, p3, p4, ncol=2)

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

데이터를 구분하기 전, 일부 함수는 입력 데이터의 변수명에 특수문자가 포함되면 에러를 일으킨다. 이에 따라, 변수명에 포함된 특수문자를 일반 문자로 바꾸자.

# 변수명의 특수문자 처리 old_names <- names(data) new_names <- make.names(names(data), unique = TRUE) cbind(old_names, new_names) [old_names!=new_names, ] names(data) <- new_names

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. 로지스틱 회귀모형(glm)

> data_lm_full <- glm(class ~ ., data=training, family=binomial) Warning message: glm.fit: 적합된 확률값들이 0 또는 1 입니다  > summary(data_lm_full) Call: glm(formula = class ~ ., family = binomial, data = training) Deviance Residuals:      Min       1Q   Median       3Q      Max   -3.9292  -0.1975   0.0000   0.1110   4.5890   Coefficients:                              Estimate Std. Error z value Pr(>|z|)     (Intercept)                -1.642e+00  1.914e-01  -8.582  < 2e-16 *** word_freq_make             -1.528e-01  2.950e-01  -0.518 0.604488     ...생략... > predict(data_lm_full, newdata = data[1:5,], type='response')         1         2         3         4         5  0.7413183 0.9855739 0.9999811 0.7588621 0.7586058 #lm 모형 평가 > y_obs <- as.numeric(as.character(validation$class)) > yhat_lm <- predict(data_lm_full, newdata = validation, type='response') > pred_lm <- prediction(yhat_lm, y_obs) > performance(pred_lm, "auc")@y.values[[1]] [1] 0.9698248 > binomial_deviance(y_obs, yhat_lm) [1] 446.264

lm의 AUC = 0.9698, 이항편차 = 446.264이다.

6. 라쏘 모형 적합(glmnet)

> xx <- model.matrix(class ~ .-1, data) > x <- xx[training_idx, ] > y <- as.numeric(training$class) > glimpse(x)  num [1:2760, 1:57] 0 0 0 0 0 0.71 0 0 0 0 ...  - attr(*, "dimnames")=List of 2   ..$ : chr [1:2760] "3372" "2374" "3391" "1971" ...   ..$ : chr [1:57] "word_freq_make" "word_freq_address" "word_freq_all" "word_freq_3d" ... > data_cvfit <- cv.glmnet(x, y, family = "binomial") > plot(data_cvfit)

왼쪽의 세로 점선은 가장 정확한 예측값을 낳는 lambda.min(52개), 오른쪽 점선은 해석 가능한 모형을 위한 lambda.1se(48개) 값을 나타낸다.

coef.cv.glmnet() 함수를 이용하여 각각의 lambda 값에서 선택된 변수들을 보여준다.

> coef(data_cvfit, s = c("lambda.1se")) ...생략... > coef(data_cvfit, s = c("lambda.min"))  ...생략...

glmnet 모형에 대한 평가는 다음과 같다.

> predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,])               1 3372  -6.885280 2374  -5.531659 3391 -18.915527 1971  -3.342803 4359 -10.318129 > yhat_glmnet <- predict(data_cvfit, s="lambda.min", newx=xx[validate_idx,], type='response') > yhat_glmnet <- yhat_glmnet[,1] # change to a vector from [n*1] matrix > pred_glmnet <- prediction(yhat_glmnet, y_obs) > performance(pred_glmnet, "auc")@y.values[[1]] [1] 0.9764865 > binomial_deviance(y_obs, yhat_glmnet) [1] 381.4467

glmnet의 AUC = 0.976, 이항편차 = 381.446이다.

7. 랜덤 포레스트(RF)

> set.seed(1810) > data_rf <- randomForest(class ~ ., training) > data_rf Call:  randomForest(formula = class ~ ., data = training)                 Type of random forest: classification                      Number of trees: 500 No. of variables tried at each split: 7         OOB estimate of  error rate: 5.04% Confusion matrix:      0    1 class.error 0 1612   50  0.03008424 1   89 1009  0.08105647 > opar <- par(mfrow=c(1,2)) > plot(data_rf) > varImpPlot(data_rf) > par(opar)

약 n.trees = 500이면 충분히 작은 오차율을 달성 할 수 있으며,

평균지니지수 감소량으로부터는, RF 모형에서는 char_freq_!, $ 표시의 개수와 remove, free 등의 단어 빈도수 등이 중요한 예측 변수임을 알 수 있다.

RF에 대한 모형평가는 다음과 같다.

> yhat_rf <- predict(data_rf, newdata=validation, type='prob')[,'1'] > pred_rf <- prediction(yhat_rf, y_obs) > performance(pred_rf, "auc")@y.values[[1]] [1] 0.983322 > binomial_deviance(y_obs, yhat_rf) [1] 298.6843

RF의 AUC = 0.983, 이항편차 = 0.298이다.

8. 부스팅(gbm)

set.seed(1810) #만약 best_iter값 계산시 오류가 발생되면, factor가 아닌 int형으로 바꿔서 해볼 것 #training$class <- as.integer(ifelse(training$class == 0, 0, 1)) glimpse(training) data_gbm <- gbm(class ~ ., data=training, distribution="bernoulli",                 n.trees=30000, cv.folds=3, verbose = TRUE) > (best_iter = gbm.perf(data_gbm, method="cv"))  [1] 1038

가장 적당한 트리의 개수는 1038개이다.

gbm의 모형 평가는 다음과 같다.

> yhat_gbm <- predict(data_gbm, n.trees=best_iter, newdata=validation, type='response') > pred_gbm <- prediction(yhat_gbm, y_obs) > performance(pred_gbm, "auc")@y.values[[1]] [1] 0.9827549 > binomial_deviance(y_obs, yhat_gbm) [1] 283.9427

gbm의 AUC = 0.982, 이항편차 = 283.9427 이다.

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

> data.frame(method=c('lm', 'glmnet', 'rf', 'gbm'), +            auc = c(performance(pred_lm, "auc")@y.values[[1]], +                    performance(pred_glmnet, "auc")@y.values[[1]], +                    performance(pred_rf, "auc")@y.values[[1]], +                    performance(pred_gbm, "auc")@y.values[[1]]), +            bin_dev = c(binomial_deviance(y_obs, yhat_lm), +                        binomial_deviance(y_obs, yhat_glmnet), +                        binomial_deviance(y_obs, yhat_rf), +                        binomial_deviance(y_obs, yhat_gbm)))   method       auc  bin_dev 1     lm 0.9698248 446.2640 2 glmnet 0.9764865 381.4467 3     rf 0.9833220 298.6843 4    gbm 0.9827549 283.9427

gbm이 가장 좋은 성능을 보인다. 조금 더 자세히 보기 위해 ROC 곡선을 그려보면 다음과 같다.

perf_lm <- performance(pred_lm, measure = "tpr", x.measure = "fpr") perf_glmnet <- performance(pred_glmnet, measure="tpr", x.measure="fpr") perf_rf <- performance(pred_rf, measure="tpr", x.measure="fpr") perf_gbm <- performance(pred_gbm, measure="tpr", x.measure="fpr") plot(perf_lm, col='black', main="ROC Curve") plot(perf_glmnet, add=TRUE, col='blue') plot(perf_rf, add=TRUE, col='red') plot(perf_gbm, add=TRUE, col='cyan') abline(0,1) legend('bottomright', inset=.1,        legend=c("GLM", "glmnet", "RF", "GBM"),        col=c('black', 'blue', 'red', 'cyan'), lty=1, lwd=2)

최종 모형을 gbm으로 선택 후, 테스트셋을 통해 이항편차와 AUC를 구해보면 다음과 같다.

> y_obs_test <- as.numeric(as.character(test$class)) > yhat_gbm_test <- predict(data_gbm, n.trees=best_iter, newdata=test, type='response') > pred_gbm_test <- prediction(yhat_gbm_test, y_obs_test) > performance(pred_gbm_test, "auc")@y.values[[1]] [1] 0.9846569 > binomial_deviance(y_obs_test, yhat_gbm_test) [1] 269.3343

회귀 분석의 오차의 시각화는 다음과 같다.

boxplot(list(lm = y_obs-yhat_lm,              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 모형의 오차가 가장 적다.

10. 예측값 시각화

pairs(data.frame(y_obs=y_obs,                  yhat_lm=yhat_lm,                  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과 lm모형의 예측 확률값들의 상관관계가 높다.

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

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