[R] 중고 자동차 데이터를 활용한 탐색적 자료분석 실전 예시

1. 탐색적 자료분석

 이번 포스팅에서는 중고자동차가격 데이터에 대해 EDA(탐색적 자료분석)를 진행해보도록 하겠습니다. 탐색적 자료분석(EDA, Exploratory Data Analysis)은 데이터 분석에 있어서 매우 중요한 step입니다. 절차는 본 분석에 진행되기에 앞서 어떻게 데이터를 분석할건지 탐색하는 것입니다. 사회에서 가장 통계분석을 못하는 유형은 탐색적 분석 절차를 거치지 않고 바로 모델링부터 도전하려는 사람들입니다. 이 글을 읽으시는 분들은 그런 실수를 안하시기를 바랍니다.

탐색적 자료분석은 정답이 있는 것이 아닙니다. 그저 짧은 시간을 투자해 최대한의 정보를 뽑아, 데이터에 대해 이해를 하는 단계라고 보시면 됩니다. 다양한 방법들이 존재하지만, 제가 주로 하는 방법은 ’시각화’입니다. 시각화만큼 데이터를 바로 이해하기 좋은 방법은 없습니다. 시각화를 주로 사용하지만, 필요에 따라 간단하게 Linear Model들 또한 활용을 합니다. Linear Model이라 하면 t-test, 분산분석, 회귀분석 등이 해당이 됩니다.

2. 데이터 및 패키지 로딩

library(ggplot2)
library(dplyr)
library(reshape)
CAR = read.csv("D:\\Dropbox\\DATA SET(Dropbox)\\car-price-predictions\\1.04. Real-life example.csv")
colnames(CAR)[1] = "Brand"
str(CAR)

'data.frame':   4345 obs. of  9 variables:
 $ Brand       : Factor w/ 7 levels "Audi","BMW","Mercedes-Benz",..: 2 3 3 1 6 3 2 1 5 7 ...
 $ Price       : num  4200 7900 13300 23000 18300 ...
 $ Body        : Factor w/ 6 levels "crossover","hatch",..: 4 6 4 1 1 1 4 5 5 3 ...
 $ Mileage     : int  277 427 358 240 120 0 438 200 193 212 ...
 $ EngineV     : num  2 2.9 5 4.2 2 5.5 2 2.7 1.5 1.8 ...
 $ Engine.Type : Factor w/ 4 levels "Diesel","Gas",..: 4 1 2 4 4 4 2 1 1 2 ...
 $ Registration: Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 1 ...
 $ Year        : int  1991 1999 2003 2007 2011 2016 1997 2006 2012 1999 ...
 $ Model       : Factor w/ 312 levels "1 Series","100",..: 19 267 239 225 228 146 19 67 190 150 ...

Brand	Price	Body	Mileage	EngineV	Engine.Type	Registration	Year	Model
BMW	4200	sedan	277	2.0	Petrol	yes	1991	320
Mercedes-Benz	7900	van	427	2.9	Diesel	yes	1999	Sprinter 212
Mercedes-Benz	13300	sedan	358	5.0	Gas	yes	2003	S 500
Audi	23000	crossover	240	4.2	Petrol	yes	2007	Q7
Toyota	18300	crossover	120	2.0	Petrol	yes	2011	Rav 4
Mercedes-Benz	199999	crossover	0	5.5	Petrol	yes	2016	GLS 63

colSums(is.na(CAR))

       Brand        Price         Body      Mileage      EngineV 
           0          172            0            0          150 
 Engine.Type Registration         Year        Model 
           0            0            0            0 

CAR = na.omit(CAR)

데이터는 꽤 간단합니다. 우리의 목적은 차량 가격을 예측하는 것입니다. 차량가격을 나타내는 변수인 Price를 중심으로 EDA를 진행해보도록 하겠습니다. 또한 결측치가 존재하기 때문에 결측치가 포함된 데이터는 삭제를 해주겠습니다.

3. 브랜드에 따른 차량가격 시각화

 먼저 가장 관심있는 브랜드에 따라 차량가격의 분포가 다른지 확인해보겠습니다. 브랜드는 범주형 자료, 차량 가격은 연속형 자료이기 때문에, 이에 적합한 시각화 방법을 선택해야됩니다. 이럴 때, 박스플롯은 탐색적분석에서 매우 유용한 그래프입니다.

ggplot(CAR) +
    geom_boxplot(aes(x = Brand, y = Price), outlier.colour = "red") +
    theme_bw() +
    theme(text = element_text(size = 15, face = "bold"),
              axis.text.x = element_text(angle = 90))

빨간 점은 outlier(이상점)을 의미합니다. 일반적인 범위에서 벗어나는 차량 가격들입니다 메르세데즈 밴츠의 경우가 상당히 이상점이 많은 것을 확인할 수가 있습니다.

ggplot(CAR) +
    geom_boxplot(aes(x = Brand, y = Price), outlier.colour = "red") +
    theme_bw() +
    theme(text = element_text(size = 15, face = "bold"),
              axis.text.x = element_text(angle = 90))  +
    facet_wrap(~Body)

이번에는 차량의 Body에 따라 어떤지 살펴보겠습니다. 어떤 Body의 차량들이 많이 고가인지 살펴보니 주로 crossover, sedan, van, other등에 위치한 것을 볼 수 있습니다. 근데, 전체적으로 가격이 높은 Body는 crossover입니다. 다른 Body들은 이상점이 높은 건데, crossover은 그냥 상자자체가 매우 큰 것을 확인할 수가 있습니다.

여기서 하나 확인할 수는 있는 것이 기본적으로 차량가격에 이상점이 매우 많다는 것입니다. 이상점이 많을 경우, 통계적인 분석에 많은 error가 발생할 수가 있습니다. 그러므로 다음 단계에서는 Price에 대한 분포를 살펴보도록 하겠습니다.

4. 차량가격 분포 탐색

ggplot(CAR) +
  geom_histogram(aes(x = Price, y = ..density..), 
                 fill = 'royalblue',col ='black',alpha = 0.4) +
  geom_density(aes(x = Price),col = 'red', size = 1.2) +
  xlab("Price") + ylab("Density") +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"))

차량가격에 대한 분포를 살펴보시면, 매우 오른쪽으로 꼬리가 치우친 분포인 것을 확인할 수가 있습니다. 이런 데이터는 정규분포와 매우 거리가 먼 상황입니다. 모든 데이터가 정규분포일 필요는 없지만, 적어도 분석하고자 하는 Target변수는 정규분포를 따라주는게 분석의 피로도를 줄일 수 있습니다.

qqnorm(CAR$Price)

qqplot을 그려보면 역시나, 정규분포와 거리가 멀다는 것을 확인할 수가 있습니다. 근데 분포가 오른쪽으로 꼬리가 치우친 경우는 사실 매우 간단하게 처리를 할 수가 있습니다. log변환을 주면 오른쪽으로 매우 치우친 분포를 중앙으로 집결시킬 수가 있습니다.

## Log Transform

ggplot(CAR) +
  geom_histogram(aes(x = log(Price), y = ..density..), 
                 fill = 'royalblue',col ='black',alpha = 0.4,binwidth = 0.1) +
  geom_density(aes(x = log(Price)),col = 'red', size = 1.2) +
  xlab("Price") + ylab("Density") +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"))

qqnorm(log(CAR$Price))

로그변환 한번 했을 뿐인데, 매우 정규분포와 유사해졌습니다. log변환된 값이기는 하지만, 분석을 할 때 이 log변환 된 값을 중심으로 분석을 하는 것이 더 통계적으로 유의한 결과들을 만들어낼 수 있는 가능성이 높아집니다. 가능성이 높아진다고 하는 얘기는, 통계분석에는 정해진 정답이 없기 때문입니다. “이 방법이 제일 좋다”라는 것은 통계분석에서 존재하지 않습니다. 즉 상황마다 적합한 분석방법이 모두 존재하는 것입니다.

5. Brand와 Body의 Price에 대한 Interaction effect 분석

 교호효과는 두 범주형 변수들의 수준(levels)조합에 따라 Target 변수의 기댓값이 달라지는가 나타내는 효과입니다. 주로 하나의 범주형 변수 수준(levels)에 따라 Target 변수들의 기댓값차이를 보는데, 이럴 때는 주효과(Main Effect)라고 합니다. 교호효과를 보는 방법은 간단합니다. 변수의 조합에 따라 기댓값(평균)의 변화를 나타내는 그래프를 그려주면 됩니다.

# Log
CAR %>%
    group_by(Brand, Body) %>%
    summarise(Mean = mean(log(Price))) %>%
    ggplot() +
    geom_point(aes(x = Brand, y = Mean, col = Body)) + 
    geom_line(aes(x = Brand, y = Mean, col = Body, group = Body)) + 
    theme_bw() + 
    theme(text = element_text(size = 15, face = "bold"),
              axis.text.x = element_text(angle = 90))

교호효과를 확인하는 좋은 방법은 선들이 만나는가 안만나는가입니다. 즉, 선들이 변수에 따라 평행관계를 유지하고 있으면 딱히 교호효과가 없다고 볼 수 있습니다. 예를 들어 crossover, other에 대한 수준은 Brand별로 비슷한 패턴을 보이고 있습니다. 하지만, hatch를 살펴보시면 밴츠와 미츠비씨에서 아주 꼴아박고 있다는 것을 확인할 수가 있습니다. 꼴아박으면서 다른 선들과 겹치게 되죠. 이런 상황이 교호효과가 존재할 수 있다는 합리적 의심을 가질 수 있게 도와줍니다. 교호효과가 정말 있는지 통계적으로 검정하는 방법은 분산분석을 활용하면 됩니다.

분산분석 이론 : https://mustlearning.tistory.com/14

## Anova

ANOVA = aov(log(Price) ~ Brand + Body + Brand:Body, data = CAR)
summary(ANOVA)

              Df Sum Sq Mean Sq F value Pr(>F)    
Brand          6  331.4   55.23  93.232 <2e-16 ***
Body           5  715.5  143.11 241.560 <2e-16 ***
Brand:Body    29   96.7    3.33   5.628 <2e-16 ***
Residuals   3984 2360.2    0.59                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Two Way ANOVA를 진행해본 결과, 역시나 주효과와 교호효과가 모두 존재하는 것을 확인할 수가 있습니다. 분산분석표에 대해서는 따로 설명을 적지는 않겠습니다.
분산분석의 결과가 유의할 경우 사후검정을 같이 해주는 것이 국룰이지만, 지금은 범주형 자료들의 수준이 너무 많기 때문에, 해봤자 머리만 아파집니다. 효율이 없으니 넘어가도록 합니다.

6. 자동차 모델 및 주행거리에 따른 Price 분석

300여개의 자동차 모델에 대해 4025개의 중고자동차 가격이 있는 데이터입니다. 그런 관계로 우리는 자동차 모델에 따라 값을 보기위해서는 데이터를 요약해주어야 합니다.

CAR %>%
    group_by(Model, Brand) %>%
    summarise(Mean_Price = mean(log(Price)),
                     Mean_Mile = mean(Mileage)) %>%
    ggplot() +
    geom_text(aes(x = Mean_Mile, y = Mean_Price, label = Model, col = Brand)) +
    theme_bw() + 
    theme(text = element_text(size = 15, face = "bold"),
          legend.position = c(0.8,0.8))

ggplot(CAR) +
    geom_point(aes(x = Mileage, y = log(Price))) +
    geom_smooth(aes(x = Mileage, y = log(Price))) + 
    theme_bw() + 
    theme(text = element_text(size = 15, face = "bold"))

보시면 Mileage(주행거리)에 따라 가격이 하락하고 있는 것을 확인할 수가 있습니다. 이런 경우 우리는 회귀분석을 적용해볼 수가 있습니다. 참고로 만약 log변환을 하지 않은 데이터로 산점도를 그려볼 경우, 다음과 같은 그래프가 나옵니다.

ggplot(CAR) +
    geom_point(aes(x = Mileage, y = Price)) +
    geom_smooth(aes(x = Mileage, y = Price)) + 
    theme_bw() + 
    theme(text = element_text(size = 15, face = "bold"))

만약 로그변환을 안했다면, 유의한 관계를 찾기 힘들 가능성이 높았을 겁니다.

\[ log(Price)=\beta_{0}+\beta{1} Mileage\\ log(Price)=\beta_{0}+\beta{1} Mileage + \beta{2}Mileage^{2} \]

그래프 상에서 비선형관계도 보이므로 일반선형회귀와 다항회귀분석을 진행해보겠습니다.

회귀분석 설명 : https://mustlearning.tistory.com/14

# Simple Linear
Reg = lm(log(Price) ~ Mileage, data = CAR)
summary(Reg)

Call:
lm(formula = log(Price) ~ Mileage, data = CAR)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8327 -0.4511  0.0226  0.4749  6.0869 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.2830105  0.0221983  463.23   <2e-16 ***
Mileage     -0.0053392  0.0001147  -46.54   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7524 on 4023 degrees of freedom
Multiple R-squared:   0.35, Adjusted R-squared:  0.3498 
F-statistic:  2166 on 1 and 4023 DF,  p-value: < 2.2e-16

\(R^2=34.98\%\)가 나왔습니다. 변수 하나만을 투입했는데 이정도면 나쁜 것은 아니라고 봅니다.

# Polynomial

Poly_Reg = lm(log(Price) ~ poly(Mileage,2), data = CAR)
summary(Poly_Reg)

Call:
lm(formula = log(Price) ~ poly(Mileage, 2), data = CAR)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0717 -0.4301  0.0498  0.4844  2.6382 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         9.40967    0.01146  821.01   <2e-16 ***
poly(Mileage, 2)1 -35.01901    0.72713  -48.16   <2e-16 ***
poly(Mileage, 2)2  12.29046    0.72713   16.90   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7271 on 4022 degrees of freedom
Multiple R-squared:  0.3931,    Adjusted R-squared:  0.3928 
F-statistic:  1303 on 2 and 4022 DF,  p-value: < 2.2e-16

\(R^2=39.28\%\)가 나왔습니다. 역시 선형회귀보다는 더 적합을 잘하는 것을 확인할 수가 있습니다. 다항회귀 항차를 더 높여볼 수도 있지만, 그래프 상에서 2차함수꼴을 보였으니 굳이 더 올려서 분석해볼 필요는 없을 것 같습니다. 원하시는 분은 한번 항차를 추가해서 값이 어떻게 변하는지 확인을 해보셔도 될 것 같습니다.

7. 엔진과 차량 가격 분석

 이번에는 엔진 변수를 활용해서 차량가격을 분석해보도록 하겠습니다.

ggplot(CAR) +
    geom_point(aes(x = EngineV, y = log(Price))) + 
    theme_bw() +
    theme(text = element_text(size = 15, face = "bold"))

EngineV와 Price간의 산점도인데, 매우 괴상하게 생겼습니다. 이런 그래프를 만나면 한가지 팁이 있습니다. x축의 범위를 조절해서 그래프를 다시 그려보는 것입니다ㅣ. 그 이유는 x범위는 매우 넓은 것에 비해 대부분의 데이터는 한쪽에 몰려있기 때문입니다.

ggplot(CAR) +
    geom_point(aes(x = EngineV, y = log(Price))) + 
    theme_bw() +
    theme(text = element_text(size = 15, face = "bold")) +
    scale_x_continuous(limits = c(0,10))

그래프를 잘라서 보니 다른 관계가 보입니다. 탐색적 자료분석의 핵심은 이런 부분에 있습니다. 데이터를 쪼개서도 보고 하나하나 살펴보며, 어떤 함정이 숨어있는지 발견하는 것, 이게 바로 데이터 분석을 정확하게 할 수 있는 실력의 지름길이 됩니다.

ggplot(CAR) +
    geom_boxplot(aes(x = Engine.Type, y = log(Price), fill = Engine.Type), alpha = 0.4, outlier.colour = 'red') +
    geom_jitter(aes(x = Engine.Type, y = log(Price), col = Engine.Type), alpha = 0.4)  +
    theme_bw() +
    theme(text = element_text(size = 15, face = "bold"))

다음으로는 엔진 타입에 대해 그래프를 그렸습니다. 역시 엔진 타입에 따라서도 가격이 다른 것을 확인할 수가 있습니다. 범주형 변수의 수준에 따른 기댓값(평균)의 차이가 있는지 확인하기 위하여 다시 한번 분산분석을 진행해보도록 하겠습니다.

ANOVA = aov(log(Price) ~ Engine.Type, data = CAR)
summary(ANOVA)

              Df Sum Sq Mean Sq F value   Pr(>F)    
Engine.Type    3     44   14.78   17.18 4.41e-11 ***
Residuals   4021   3460    0.86                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

분산분석 결과 차이가 매우 유의하다는 것을 확인할 수가 있습니다. 이번에는 수준 수가 4개이기 때문에 사후검정을 진행해보도록 하겠습니다.

TUKEY = TukeyHSD(ANOVA)
TUKEY 

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = log(Price) ~ Engine.Type, data = CAR)

$Engine.Type
                     diff         lwr         upr     p adj
Gas-Diesel    -0.24113612 -0.35376896 -0.12850328 0.0000003
Other-Diesel  -0.25238630 -0.49043654 -0.01433606 0.0327292
Petrol-Diesel  0.05838048 -0.02483674  0.14159771 0.2719751
Other-Gas     -0.01125018 -0.26273875  0.24023839 0.9994567
Petrol-Gas     0.29951661  0.18331123  0.41572198 0.0000000
Petrol-Other   0.31076679  0.07100555  0.55052802 0.0048353

plot(TUKEY)

사후검정결과가 나왔습니다. 좀 자세히 정리해볼 필요가 있겠습니다. 각 차이에 대한 95% 신뢰구간이 0을 포함하는지 안하는지만 확인하시면 됩니다. 그러면 다음과 같이 정리가 가능합니다.

\[ Gas-Diesel < 0 \\ Other-Diesel < 0 \\ Petrol-Diesel = 0 \\ Other-Gas = 0 \\ Petrol-Gas > 0 \\ Petrol-Other > 0 \]

이번에는 엔진과 엔진타입을 포함하여 가변수 회귀분석을 진행해보도록 하겠습니다.

Reg2 = lm(log(Price) ~ Engine.Type + EngineV, data = CAR)
summary(Reg2)

Call:
lm(formula = log(Price) ~ Engine.Type + EngineV, data = CAR)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.08079 -0.55506 -0.05422  0.57677  3.10458 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        9.408607   0.022757 413.431  < 2e-16 ***
Engine.TypeGas    -0.243758   0.043793  -5.566 2.77e-08 ***
Engine.TypeOther  -0.272324   0.092793  -2.935  0.00336 ** 
Engine.TypePetrol  0.055358   0.032366   1.710  0.08727 .  
EngineV            0.008598   0.002969   2.896  0.00379 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9267 on 4020 degrees of freedom
Multiple R-squared:  0.01471,   Adjusted R-squared:  0.01373 
F-statistic:    15 on 4 and 4020 DF,  p-value: 3.54e-12

만약 EngineV에 데이터 범위를 조정하지 않을 경우 회귀분석의 결과는 정말 처참하게 나옵니다. 그렇기 때문에 EngineV를 활용해 데이터를 추출해보고 다시 돌려보도록 하겠습니다.

CAR2 = subset(CAR, EngineV < 10)
Reg2s = lm(log(Price) ~ Engine.Type + EngineV, data = CAR2)
summary(Reg2s)

Call:
lm(formula = log(Price) ~ Engine.Type + EngineV, data = CAR2)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.11037 -0.47377  0.01656  0.53087  2.39005 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        8.40716    0.03649 230.402   <2e-16 ***
Engine.TypeGas    -0.39544    0.03916 -10.098   <2e-16 ***
Engine.TypeOther  -0.19139    0.08363  -2.289   0.0222 *  
Engine.TypePetrol -0.05580    0.02898  -1.925   0.0542 .  
EngineV            0.44104    0.01336  33.008   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8222 on 4000 degrees of freedom
Multiple R-squared:  0.224, Adjusted R-squared:  0.2232 
F-statistic: 288.6 on 4 and 4000 DF,  p-value: < 2.2e-16

\(R^{2}\)가 나름 22%까지 상승한 것을 볼 수가 있습니다. 이정도면 엔진 타입도 효과를 준다고 볼 수가 있습니다. Petrol Type에 대한 회귀계수 검정의 p-value가 0.0542로 0.05보다 높지만, 아주 높은 것도 아니며, 탐색적 자료분석은 Rough하게 분석하는 것이 주요 point입니다. 처음부터 너무 잔혹하게 볼 필요는 또 없는 것입니다.

여기서 우리에게 주어진 과제는 다음과 같습니다.

• EngineV를 어떻게 처리를 해야되는 것인가?
• Mileage를 어떻게 처리를 해야하는 것인가?

EngineV는 바로 위에서 사례를 봤으며, Mileage는 Mileage가 큰 경우에 한해 선형 데이터가 비선형으로 꺾여버리는 문제점이 존재하기 때문입니다. 이런 상황은 사실 데이터를 예측 모델링 혹은 관계를 설명하기 위한 모델링을 진행하는데 매우 방해가 됩니다. 그렇기 때문에, 값을 자르고 갈 것인지, 안고 갈 것인지 선택을 해야합니다. 이건 데이터 분석가들의 생각에 따라 갈리게 됩니다. 즉 자신이 맞다고 생각하는 방향으로 진행하면 됩니다. 대신에 그 이유는 충분히 설명할 수는 있어야합니다(“분석은 설득입니다”).

이렇게 차량의 특징들이 차량 가격에 어떤 관계를 가지고 있는지 확인을 하는 탐색적 자료분석을 진행했습니다. 많은 도움이 되셨으면 바라고, 다음 포스팅은 예측 modeling을 어떻게 진행해야하는지에 대해 다루어보도록 하겠습니다.