1 미국대선 여론조사 1 2

제45대 미국 대통령을 선출하는 선거에서 많은 언론에서 예측한 바와 다르게 공화당 도널드 트럼프 후보가 민주당 힐러리 클린턴 후보를 누르고 당선되었다. 이를 두고 트럼프 후보 당선을 정확히 예측한 인공지능(AI) 사례를 들며 기존 예측기법에 대한 문제점 제기 및 새로운 시대의 도래를 언급하고 있다. 이에 앞서 지금과 동일한 상황이 미국에서 여러번 있어와서 새삼스러운 것은 아니다.

  • 1916년 대선: 우편 설문조사를 바탕으로 민주당 우드로 윌슨 후보 당선 예측한
    • 리터리리 다이제스트는 1936년 대통령 예측에서 실패
    • 1천만명 유권자를 대상으로 설문조사를 실시했으나 표본의 대표성에 문제가 있음
    • 주소 확보가 용이한 구독자, 자동차 등록부, 전화번호부 등 공화당 성향 유권자에 집중
  • 1936년 대선: 조기 갤럽의 ’미국공공여론연구소’에서 5만명 표본으로 프랭클린 루스벨트 당선 예측
    • 인구학적 분포에 대응되는 표본 추출 기법 차용
    • 1948년 선거일보다 3주 앞서 당선 후보 예측한 것이 화근이 되어, ’Dewey defeats truman’이라는 ’역사적 오보’를 남김.
  • 패널 조사: 장기간의 추적분석과 정치 성향이 표본 추출에 내재된 보다 정교한 여론조사 기법이 등장
    • 2016년 미국 대선에서 non-ignorable nonresponse 문제를 해결하지 못해 예측에 실패함.
    • 인공지능 인도 벤처기업 제니크.ai(Genic.ai) 인공지능 모그IA(MogIA)는 빅데이터 분석을 통해 10월부터 트럼프 승리 예측

결측값 구분 3

  • MCAR(Missing Completely At Random): 모든 정보가 데이터에 담겨있어 결측값이 문제가 되지 않는 경우.
  • MAR(Missing At Random): 더 이해하기 쉬운 명칭은 Missing Conditionally at Random으로 결측조건이 다른 변수에 따라 조건부로 발생되는 경우. 결측값이 관측된 데이터가 아닌 관측되지 않는 데이터에 따라 결정.
    • 예를 들어, 남성이 여성보다 더 솔찍하게 나이, 몸무게 등을 밝힐 듯 싶다.
  • MNAR(Missing Not At Random): 결측값이 무작위가 아니라서 주도면밀한 추가 조사가 필요한 경우.
    • 예를 들어, 작업장에서 결측된 사람들 대부분은 아마도 몸이 아파 조사에서 결측되고, 교육수준이 낮은 사람이 교육조사에 결측될 듯 싶다.

데이터 생성과정을 이해하고 각 단계별로 왜 데이터에 누락이 발생했는지 따진다. 예를 들어, 미국 대선에서 왜 일부 유권자가 설문조사에 응답을 거부했는지 파악하고, 설문 항목에 문제가 없는지, 설문 문항에 불명확한 점이 있는지 다각도로 조사한다.

2 R에서 NA

결측값은 기록되어야만 하지만 기록되지 않는 값을 의미한다. NA**N**ot **A**vailable의 약자다. naniar 팩키지는 요인형 자료형 처리를 위해 개발된 forcats와 유사하게 NA 결측값 데이터 처리를 위해 개발된 팩키지로 자리잡고 있다.

2.1 결측값 파악

먼저 특정 벡터에 결측값이 있는지 여부는 any_na 함수로 파악하고, 해당 벡터의 원소가 결측값인지는 are_na 함수로 파악한다.

그리고 n_miss(), prop_miss() 함수로 결측갯수와 결측비율도 산출이 가능하다.

NA 파악

library(naniar)
missing_v <- c(7, NA, 17, NA, NA, 25, 36, NA, NA, 100)

any_na(missing_v)
[1] TRUE
are_na(missing_v)
 [1] FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE

NA 갯수와 비율

n_miss(missing_v)
[1] 5
prop_miss(missing_v)
[1] 0.5

2.2 결측값 연산

컴퓨터 과학에서 일반적으로 결측값(missing value)은 존재하지 않는 값(null)으로 컴퓨터에서 표현된다. 하지만 데이터 과학에서는 결측값이 무응답 혹은 단순히 자리만 차지하는 NA, \(\frac {0}{0}\)을 표현하는 NaN (Not a Number), 무한을 표현하는 inf가 있다.

R에 NA와 같이 미리 예약된 예약어가 몇개 있다. NaN, NULL, Inf가 대표적인다.

  • NULL: 값이 없다로 즉, 값이 존재하지 않는다는 의미로 사용
  • NaN: Not a Number로 \(\frac{0}{0}\)처럼 수학적으로 정의가 되지않는다는 의미.
  • Inf: Infinite의 약자로 무한대를 의미하는데, 수학적으로 \(\frac{10}{0}\)은 무한대가 된다.

NaN, NULL, Inf 각각이 결측값(NA)인지 살펴보자. any_na(NULL)은 TRUE로 결측값이 되고, 나머지는 결측값이 아니다. 결측값의 정의에서 기록되어야 하지만 기록되지 않은 값이 결측값이기 때문에 정의에 맞춰 잘 계산됨이 확인된다.

any_na(NaN)
[1] TRUE
any_na(NULL)
[1] FALSE
any_na(Inf)
[1] FALSE

다음으로 R에서 특별히 정의된 예약어에 대해 연산을 수행하면 어떤 값이 되는지 확인해 보자.

NA에 TRUE/FALSE를 적용하면 TRUE/FALSE에 따라 값이 달라지고, NANaN를 연산하면 연산순서에 따라 결과가 달라진다.

NA | TRUE
[1] TRUE
FALSE | NA
[1] NA
NA + NaN
[1] NA
NaN + NA
[1] NaN

3 결측값 식별과 현황파악 4 5

결측데이터를 처리하기 전에 결측데이터를 처리하는 프로세스는 다음과 같다.

  1. 결측값을 식별한다. 원본데이터에서 다양한 형태로 결측정보가 표현되어 있으니 우선 현황 파악이 먼저다.
  2. 파악된 현황정보를 바탕으로 결측값을 컴퓨터가 처리가능한 형태로 부호화한다.
  3. 마지막으로 파악된 결측정보와 적절히 인코딩되어 컴퓨터에 저장된 결측값을 자료형에 맞춰 알고리즘을 적용하여 결측값을 처리한다.

결측값 대체 프로세스

3.1 결측데이터 생성

mlbench 보스톤 주택가격 데이터셋을 기본으로 결측값 관련 학습 내용에 대한 실습을 진행한다. 이를 위해서 연속형 변수 “ptratio”에 30개 NA 결측값을 주입하고, 범주형 변수 “rad”에 50개 NA 결측값를 주입한다.

## 실습데이터

library(mlbench)
library(tidyverse)
data("BostonHousing")

# 원데이터를 나중에 복구하기 위해 잠시 `original` 데이터프레임으로 저장
original <- BostonHousing

# "ptratio 연속형, "rad"범주형 
BostonHousing[sample(1:nrow(BostonHousing), 30), "ptratio"] <- NA
BostonHousing[sample(1:nrow(BostonHousing), 50), "rad"] <- NA

3.2 결측값 현황 파악

결측값 현황 파악을 위해서 변수기준으로 결측값이 얼마나 되는지 파악하는 것과, 관측점 기준으로 얼마나 결측값이 발생되었는지 파악하는 두가지 경우가 있다.

먼저, apply 계열 함수를 사용해서 변수별, 관측점별 결측값이 얼마나 포함되어 있는지 파악한다. check_missing_value() 함수를 활용하여 변수별, 관측점별 결측값 백분율을 계산한다. 만약 5%이상 관측점 혹은 변수에 결측값이 포함되어 있으면 변수를 버리거나 혹은 관측점을 제거하는 전략을 택할 수도 있다.

다음으로 naniar 팩키지에 내장된 함수를 사용하는 방법이 있다.

  • 전반적인 결측값 현황 파악
    • n_miss()
    • n_complete()
  • 변수와 관측점별 결측값 파악
    • miss_var_summary()
    • miss_case_summary()
  • 변수와 관측점 교차표
    • miss_var_table()
    • miss_case_table()
  • 시계열
    • miss_var_span()
    • miss_var_run()

apply 계열함수 사용

check_missing_value <- function(x){sum(is.na(x))/length(x)*100}
apply(BostonHousing, 2, check_missing_value)
    crim       zn    indus     chas      nox       rm      age      dis 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
     rad      tax  ptratio        b    lstat     medv 
9.881423 0.000000 5.928854 0.000000 0.000000 0.000000 
apply(BostonHousing, 1, check_missing_value) %>% head(20)
       1        2        3        4        5        6        7        8 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.142857 0.000000 
       9       10       11       12       13       14       15       16 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
      17       18       19       20 
0.000000 0.000000 0.000000 0.000000 

naniar 팩키지 함수

n_miss(BostonHousing)
[1] 80
n_complete(BostonHousing)
[1] 7004
miss_var_summary(BostonHousing)
# A tibble: 14 x 3
   variable n_miss pct_miss
   <chr>     <int>    <dbl>
 1 rad          50     9.88
 2 ptratio      30     5.93
 3 crim          0     0   
 4 zn            0     0   
 5 indus         0     0   
 6 chas          0     0   
 7 nox           0     0   
 8 rm            0     0   
 9 age           0     0   
10 dis           0     0   
11 tax           0     0   
12 b             0     0   
13 lstat         0     0   
14 medv          0     0   
miss_case_summary(BostonHousing)
# A tibble: 506 x 3
    case n_miss pct_miss
   <int>  <int>    <dbl>
 1    63      2    14.3 
 2   174      2    14.3 
 3   255      2    14.3 
 4   314      2    14.3 
 5     7      1     7.14
 6    23      1     7.14
 7    26      1     7.14
 8    34      1     7.14
 9    39      1     7.14
10    46      1     7.14
# … with 496 more rows

3.3 결측값 시각화

데이터프레임 행과 열에 걸쳐 분포되어 있는 결측값 현황을 새의 눈 관점(birds eye view)으로 살펴보자.

3.3.1 mice, VIM 팩키지

mice 팩키지를 사용해서 md.pattern 함수를 사용해서 “ptratio” 변수에 40개 주입한 NA 결측값을 확인해본다.

library(mice)
md.pattern(BostonHousing)

    crim zn indus chas nox rm age dis tax b lstat medv ptratio rad   
430    1  1     1    1   1  1   1   1   1 1     1    1       1   1  0
46     1  1     1    1   1  1   1   1   1 1     1    1       1   0  1
26     1  1     1    1   1  1   1   1   1 1     1    1       0   1  1
4      1  1     1    1   1  1   1   1   1 1     1    1       0   0  2
       0  0     0    0   0  0   0   0   0 0     0    0      30  50 80

혹은 VIM 팩키지를 활용하여 결측값을 시각화하는 것도 좋다. 결측값은 빨강색(red)로 설정하여 각변수별로 결측값에 대한 현황을 파악하기 쉽게 도식화했다.

library(VIM)
aggr_plot <- aggr(BostonHousing, col=c('gray','red'), numbers=TRUE, 
    sortVars=TRUE, labels=names(BostonHousing), cex.axis=.7, gap=3, ylab=c("Histogram of missing data", "Pattern"))


 Variables sorted by number of missings: 
 Variable      Count
      rad 0.09881423
  ptratio 0.05928854
     crim 0.00000000
       zn 0.00000000
    indus 0.00000000
     chas 0.00000000
      nox 0.00000000
       rm 0.00000000
      age 0.00000000
      dis 0.00000000
      tax 0.00000000
        b 0.00000000
    lstat 0.00000000
     medv 0.00000000
marginplot(BostonHousing[, c('ptratio','medv')])

marginplot(BostonHousing[, c('ptratio','medv')]) 그래프를 통해 변수 2개만 가능하지만, 결측값이 포함된 경우와 결측값이 없는 경우를 비교해 볼 수도 있다.

3.3.2 naniar 팩키지

naniar 팩키지 viz_miss() 함수로 데이터프레임 전체에 대한 시각화가 가능합니다. cluster=TRUE 선택옵션을 지정하면 전반적인 현황을 한눈에 파악할 수 있다.

bh_orig_v <- vis_miss(BostonHousing)
bh_cluster_v <- vis_miss(BostonHousing, cluster=TRUE)

cowplot::plot_grid(bh_orig_v, bh_cluster_v)

naniar 팩키지 gg_miss_var(), gg_miss_case() 함수를 사용하면 변수별, 관측점별로 결측값에 대한 정보를 더 자세히 파악할 수 있다.

bh_var_v  <- gg_miss_var(BostonHousing)
bh_case_v <- gg_miss_case(BostonHousing, show_pct = FALSE) + labs(y="결측 건수")

cowplot::plot_grid(bh_var_v, bh_case_v)

facet을 지정하여 범주형 변수 수준별로 결측변수별 결측값을 파악하는 것도 가능하다.

gg_miss_var(BostonHousing, facet = chas)

gg_miss_upset() 함수를 통해 변수와 관측점간 결측값도 동시에 시각화를 할 수 있다.

gg_miss_upset(BostonHousing)

그외에도 gg_miss_fct(), gg_miss_span() 함수로 범주형 변수와 시계열 변수들에 대한 결측값도 시각화를 할 수 있다.

4 결측값 처리

결측값을 지정하여 이를 처리할 경우 먼저 결측값을 지정해야 한다. 그리고 나서 결측값을 NA 채워넣는다.

4.1 결측값 지정

난잡한 데이터프레임을 하나 만들고 나서, 특정값 예를 들어 -99를 결측값으로 지정한 경우 -99 결측값이 포함된 변수가 어떤 것이며 결측값이 몇개나 되는지 파악하는 것이 miss_scan_count() 함수를 통해 가능하다. 두가지 이상 지정하는 것도 가능한데 search = list("-99", "-98")와 같이 저정하면 다수 결측값 패턴을 명시할 수 있다.

miss_df <- tribble(~alpha,  ~beta  , ~gamma,
              -99  ,  "E"    ,    97,
              4    ,  "missing",  95,
              -99  ,  "na"     ,  92,
              7    ,  "n/a"    ,  -98,
              10   ,  "missing",  "",
              3    ,  "N/A"   ,   -99,
              12   ,  "."     ,   88,
              16   ,  "."     ,   "",
              9    ,  "N/a"   ,   86)


miss_df %>% 
    miss_scan_count(search = list("-99", "-98"))
# A tibble: 3 x 2
  Variable     n
  <chr>    <int>
1 alpha        2
2 beta         0
3 gamma        2

4.2 결측값 채워넣기

결측값 채워넣을 때 dplyr 팩키지 변수 선택과 유사하게 다양한 방식으로 NA로 채워넣을 있다.

replace_with_na(): 특정변수 값을 지정하여 결측값으로 채워넣기 replace_with_na_all(): 모든 변수 특정 조건을 만족하는 값을 결측값으로 채워넣기 replace_with_na_at(): 변수 일부만 골라 특정 조건을 만족하는 값을 결측값으로 채워넣기 replace_with_na_if(): 특정 조건을 만족하는 변수만 골라 특정 조건을 만족하는 값을 결측값으로 채워넣기

-99, -98, "N/a" 값을 결측값으로 지정해서 replace_with_na_all() 함수로 모든 변수에 대해 결측값으로 채워넣는다.

결측값 채워 넣기 전

miss_df
# A tibble: 9 x 3
  alpha beta    gamma
  <dbl> <chr>   <chr>
1   -99 E       "97" 
2     4 missing "95" 
3   -99 na      "92" 
4     7 n/a     "-98"
5    10 missing ""   
6     3 N/A     "-99"
7    12 .       "88" 
8    16 .       ""   
9     9 N/a     "86" 

결측값 채워 넣은 후

miss_df %>% 
    replace_with_na_all(condition = ~.x %in% c(-99, -98, "N/a"))
# A tibble: 9 x 3
  alpha beta    gamma
  <dbl> <chr>   <chr>
1    NA E       "97" 
2     4 missing "95" 
3    NA na      "92" 
4     7 n/a      <NA>
5    10 missing ""   
6     3 N/A      <NA>
7    12 .       "88" 
8    16 .       ""   
9     9 <NA>    "86" 

5 결측값 처리 사례 - 옛날 방식

결측값 처리 전략은 다음과 같은 4가지 전략이 존재한다.

  1. 관측점 제거
    • 데이터가 상당히 많은 경우, 동시에 모집단을 대표하는데 무리가 없는 경우 na.action=na.omit 설정을 적용한다.
    • lm(medv ~ ptratio + rad, data=BostonHousing, na.action=na.omit)
  2. 변수 제거
    • 특정 변수에 상당한 값이 결측값인 경우, 다른 변수가 결측값을 많이 갖는 변수에 상응하는 정보량을 갖는 경우 제거한다.
  3. 평균/중위수/최빈값으로 대체(impute)
    • 결측값을 변수가 연속형 숫자형인 경우 평균/중위수로 대체하고, 범주형 요인형인 경우 최빈값으로 대체한다.
library(Hmisc)
library(DMwR)
library(dplyr)
impute(BostonHousing$ptratio, mean) %>% head # 평균으로 대체
   1    2    3    4    5    6 
15.3 17.8 17.8 18.7 18.7 18.7 
impute(BostonHousing$ptratio, median) %>% head  # 중위수로 대체
   1    2    3    4    5    6 
15.3 17.8 17.8 18.7 18.7 18.7 
impute(BostonHousing$ptratio, 20) %>% head  # 특정 값으로 대체
   1    2    3    4    5    6 
15.3 17.8 17.8 18.7 18.7 18.7 
# 혹은 팩키지를 사용하지 않고 직접 코드를 작성해서 작업
BostonHousing$ptratio[is.na(BostonHousing$ptratio)] <- mean(BostonHousing$ptratio, na.rm = T) 
  1. 예측값으로 대체한다.
  • DMwR 팩키지 knnImputation() 함수를 사용해서 k-nn (k-인접 군집분석)을 사용한다. k-nn 대체법은 인접한 최대 k 관측점 유클리드 거리를 계산하여 가장 근접된 값으로 대체한다.
library(DMwR)
data(BostonHousing)
BostonHousing[sample(1:nrow(BostonHousing), 40), "ptratio"] <- NA

knnOutput <- knnImputation(BostonHousing[, !names(BostonHousing) %in% "medv"]) 

# 효과 분석
actuals <- original$ptratio[is.na(BostonHousing$ptratio)]
predicteds <- knnOutput[is.na(BostonHousing$ptratio), "ptratio"]
regr.eval(actuals, predicteds)
       mae        mse       rmse       mape 
0.72317896 1.16829389 1.08087644 0.04104449 
  • knn 대체 기법은 범주형 자료의 경우 적용에 한계가 있다. 이런 경우 rpart, mice 팩키지를 활용한다.
  • 먼저, rpart를 사용하는 경우 연속형 변수, 범주형 변수 모두 의사결정나무 모형을 순차적으로 적합시킨다.
library(rpart)
data(BostonHousing)
BostonHousing[sample(1:nrow(BostonHousing), 40), "ptratio"] <- NA
BostonHousing[sample(1:nrow(BostonHousing), 40), "rad"] <- NA

class_mod <- rpart(rad ~ . - medv, data=BostonHousing[!is.na(BostonHousing$rad), ], method="class", na.action=na.omit)  # rad 변수가 범주형
anova_mod <- rpart(ptratio ~ . - medv, data=BostonHousing[!is.na(BostonHousing$ptratio), ], method="anova", na.action=na.omit)  # ptratio 변수는 숫자형
rad_pred <- predict(class_mod, BostonHousing[is.na(BostonHousing$rad), ])
ptratio_pred <- predict(anova_mod, BostonHousing[is.na(BostonHousing$ptratio), ])    
  • mice를 사용하는 경우: mice() 함수를 사용해서 먼저 모형을 생성시키고 나서, complete() 함수를 사용해서 결측값을 채워넣는 2단계 과정을 거침.
library(mice)
data(BostonHousing)
BostonHousing[sample(1:nrow(BostonHousing), 40), "ptratio"] <- NA
miceMod <- mice(BostonHousing[, !names(BostonHousing) %in% "medv"], method="rf")  # 확률숲(random forest) 모형으로 결측모형 생성.

 iter imp variable
  1   1  ptratio
  1   2  ptratio
  1   3  ptratio
  1   4  ptratio
  1   5  ptratio
  2   1  ptratio
  2   2  ptratio
  2   3  ptratio
  2   4  ptratio
  2   5  ptratio
  3   1  ptratio
  3   2  ptratio
  3   3  ptratio
  3   4  ptratio
  3   5  ptratio
  4   1  ptratio
  4   2  ptratio
  4   3  ptratio
  4   4  ptratio
  4   5  ptratio
  5   1  ptratio
  5   2  ptratio
  5   3  ptratio
  5   4  ptratio
  5   5  ptratio
miceOutput <- complete(miceMod)  # 생성된 데이터를 채워 넣음.
anyNA(miceOutput)    
[1] FALSE

5.1 성능 평가 – 연속형 변수 ptratio

다양한 결측값 처리 방법에 따른 성능 차이를 비교하는 것이 왜 고급 결측값 처리 방법을 활용해야 하는 근거도 된다.

먼저, 연속형 변수의 경우 다양한 결측값 처리 방법에 따른 성능의 차이를 비교해보자. “ptratio” 변수는 연속형 변수로 506개 변수중 50개 즉 10%를 결측값, NA로 치환한다.

평균과 중위수를 결측값 10%를 채워넣을 경우 중위수를 채워넣은 것이 mape를 봤을 때 대동소이하다. 6

data(BostonHousing)
original <- BostonHousing
BostonHousing[sample(1:nrow(BostonHousing), 50), "ptratio"] <- NA

actuals_ptratio <- original$ptratio[is.na(BostonHousing$ptratio)]
ptratios_mean_pred <- rep(mean(BostonHousing$ptratio, na.rm=T), length(actuals_ptratio))
ptratios_median_pred <- rep(median(BostonHousing$ptratio, na.rm=T), length(actuals_ptratio))
regr.eval(actuals_ptratio, ptratios_mean_pred)
      mae       mse      rmse      mape 
1.8482018 5.3799500 2.3194719 0.1121075 
regr.eval(actuals_ptratio, ptratios_median_pred)
      mae       mse      rmse      mape 
1.9160000 6.2288000 2.4957564 0.1194082 

두번째 knn 기법을 활용한 경우 mape가 줄어든 것이 확인된다.

ptratios_knn_pred <- knnOutput[is.na(BostonHousing$ptratio), "ptratio"]
regr.eval(actuals_ptratio, ptratios_knn_pred)
        mae         mse        rmse        mape 
0.084236878 0.185174328 0.430318867 0.005192704 

세번째 rpart 기법을 활용한 경우 mape가 줄어든 것이 확인된다.

ptratios_anova_mod <- rpart(ptratio ~ . - medv, 
    data=BostonHousing[!is.na(BostonHousing$ptratio), ], method="anova", na.action=na.omit)
ptratio_anova_pred <- predict(ptratios_anova_mod, BostonHousing[is.na(BostonHousing$ptratio), ])
regr.eval(actuals, ptratio_anova_pred)
       mae        mse       rmse       mape 
 2.8569334 10.1558011  3.1868168  0.1602098 

네번째는 rf 확률숲 모형을 적용하는데 결측값 처리 전용 mice 팩키지를 활용한다. mape 값이 하향된 것이 관측된다.

library(mice)
mice_mod <- mice(BostonHousing[, !names(BostonHousing) %in% "medv"], method="rf") # 1단계 모형 생성

 iter imp variable
  1   1  ptratio
  1   2  ptratio
  1   3  ptratio
  1   4  ptratio
  1   5  ptratio
  2   1  ptratio
  2   2  ptratio
  2   3  ptratio
  2   4  ptratio
  2   5  ptratio
  3   1  ptratio
  3   2  ptratio
  3   3  ptratio
  3   4  ptratio
  3   5  ptratio
  4   1  ptratio
  4   2  ptratio
  4   3  ptratio
  4   4  ptratio
  4   5  ptratio
  5   1  ptratio
  5   2  ptratio
  5   3  ptratio
  5   4  ptratio
  5   5  ptratio
mice_output <- complete(mice_mod)  # 2단계 결측값 채워넣기

ptratio_rf_pred <- mice_output[is.na(BostonHousing$ptratio), "ptratio"]
regr.eval(actuals, ptratio_rf_pred)
       mae        mse       rmse       mape 
 2.7175000 10.4917500  3.2390971  0.1530353 

5.2 성능 평가 – 범주형 변수 rad

범주형 변수의 경우, 다양한 결측값 처리 방법에 따른 성능의 차이를 비교해보자. “rad” 변수는 연속형 변수로 506개 변수중 50개 즉 10%를 결측값, NA로 치환한다.

먼저, names(sort(-table(BostonHousing$rad)))[1] 명령어를 통해 최빈값을 파악한다. 그리고 이를 결측값에 꽂아 넣는다.

data(BostonHousing)
original <- BostonHousing
BostonHousing[sample(1:nrow(BostonHousing), 50), "rad"] <- NA

actuals_rad <- original$rad[is.na(BostonHousing$rad)]
rad_mode_pred <- rep(names(sort(-table(BostonHousing$rad)))[1], length(actuals_rad))

mean(actuals_rad != rad_mode_pred) 
[1] 0.66

두번째로 rpart 의사결정나무 모형을 활용하여 결측값을 채워넣는다.

data(BostonHousing)
original <- BostonHousing
BostonHousing[sample(1:nrow(BostonHousing), 50), "rad"] <- NA

library(rpart)
class_mod <- rpart(rad ~ . - medv, data=BostonHousing[!is.na(BostonHousing$rad), ], method="class", na.action=na.omit)
rad_pred <- predict(class_mod, BostonHousing[is.na(BostonHousing$rad), ])

actuals_rad <- original$rad[is.na(BostonHousing$rad)]
rad_rpart_pred <- as.numeric(colnames(rad_pred)[apply(rad_pred, 1, which.max)])

mean(actuals_rad != rad_rpart_pred)
[1] 0.3

마지막으로 mice 확률숲 rf 모형을 사용해서 결측값을 채워넣는다.

data(BostonHousing)
original <- BostonHousing
BostonHousing[sample(1:nrow(BostonHousing), 50), "rad"] <- NA

library(mice)
mice_mod <- mice(BostonHousing[, !names(BostonHousing) %in% "medv"], method="rf") # 1단계 모형 생성

 iter imp variable
  1   1  rad
  1   2  rad
  1   3  rad
  1   4  rad
  1   5  rad
  2   1  rad
  2   2  rad
  2   3  rad
  2   4  rad
  2   5  rad
  3   1  rad
  3   2  rad
  3   3  rad
  3   4  rad
  3   5  rad
  4   1  rad
  4   2  rad
  4   3  rad
  4   4  rad
  4   5  rad
  5   1  rad
  5   2  rad
  5   3  rad
  5   4  rad
  5   5  rad
mice_output <- complete(mice_mod)  # 2단계 결측값 채워넣기

actuals_rad <- original$rad[is.na(BostonHousing$rad)]
rad_rf_pred <- mice_output[is.na(BostonHousing$rad), "rad"]
mean(actuals_rad != rad_rf_pred)
[1] 0.18

6 명시적, 묵시적 결측값

결측값이 명시적으로 나타난 경우와 암묵적으로 나타난 경우로 나눌 수 있다. 명시적, 묵시적 결측값 처리를 위한 사례로 tetris_df 데이터프레임을 생각해보자. 광춘, 성호, 충현이 테트리스 게임을 모두 했는데, 데이터에 보면 성호만 저녁에 게임을 한 기록이 누락된 것이 파악된다. 이런 경우 데이터는 모두 완전한 데이터지만, 묵시적인 결측값이 포함된 경우다.

spread() 함수로 긴형태 데이터프레임을 넓은 형태 데이터프레임으로 변환시키게 되면 결측값이 확인히 나타나게 된다.

묵시적 결측값 사례

tetris_df <- tribble (
    ~name,    ~time,     ~game_score,
    "광춘",  "morning",   358,
    "광춘",  "afternoon", 534,
    "광춘",  "evening",   100,
    "성호",  "morning",   139,
    "성호",  "afternoon", 177,
    "충현",  "morning",   963,
    "충현",  "afternoon", 962,
    "충현",  "evening",   929)

tetris_df
# A tibble: 8 x 3
  name  time      game_score
  <chr> <chr>          <dbl>
1 광춘  morning          358
2 광춘  afternoon        534
3 광춘  evening          100
4 성호  morning          139
5 성호  afternoon        177
6 충현  morning          963
7 충현  afternoon        962
8 충현  evening          929

명시적으로 결측값 보이게 하기

tetris_df %>% 
    spread(time, game_score)
# A tibble: 3 x 4
  name  afternoon evening morning
  <chr>     <dbl>   <dbl>   <dbl>
1 광춘        534     100     358
2 성호        177      NA     139
3 충현        962     929     963

이와 같은 결측값을 채워넣는 경우 tidyr::complete() 함수가 제격이다.

tetris_df %>%
    tidyr::complete(name, time)
# A tibble: 9 x 3
  name  time      game_score
  <chr> <chr>          <dbl>
1 광춘  afternoon        534
2 광춘  evening          100
3 광춘  morning          358
4 성호  afternoon        177
5 성호  evening           NA
6 성호  morning          139
7 충현  afternoon        962
8 충현  evening          929
9 충현  morning          963

7 결측값 효과 분석

결측값 효과 분석을 위해서 결측값이 있는 상황과 없는 상황에 대한 명확한 구별이 필요하다.

7.1 잠재 행렬(Shadow Matrix) 자료구조

먼저 잠재 행렬(Shadow Matrix)는 행렬의 모든 값을 NA, !NA로 바꿔주는 행렬이다. 이를 통해서 결측값이 있는 경우와 그렇지 않는 경우를 나눠서 분석할 수 있게 된다.

결측값이 포함된 원본 데이터

library(tidyverse)
library(naniar)

sm_df <- tribble(~income , ~education,
               48.69087,    NA,
               40.93218,    NA,
               52.69245,    "high_school",
               31.33808,    NA,
               89.35671,    "university",
               103.8727,    "university")

잠재행렬 반영한 데이터

bind_shadow(sm_df)
# A tibble: 6 x 4
  income education   income_NA education_NA
   <dbl> <chr>       <fct>     <fct>       
1   48.7 <NA>        !NA       NA          
2   40.9 <NA>        !NA       NA          
3   52.7 high_school !NA       !NA         
4   31.3 <NA>        !NA       NA          
5   89.4 university  !NA       !NA         
6  104.  university  !NA       !NA         

7.2 요약통계

잠재행렬(shadow matrix)을 통해 bind_shadow() 함수로 구성하게 되면, NA표(nabular) 데이터가 만들어지는데 이를 통해 결측값에 대한 후속 분석이 수월해진다. 즉, 결측값이 존재하는 경우와 결측값이 존재하지 않는 것을 비교하는 것이 가능하다.

airquality %>%
    bind_shadow() %>%
    group_by(Ozone_NA) %>%
    summarise(wind_mean = mean(Wind))
# A tibble: 2 x 2
  Ozone_NA wind_mean
  <fct>        <dbl>
1 !NA           9.86
2 NA           10.3 

7.3 시각화

잠재행렬로 데이터프레임이 구성되면 ggplot에 넣어 일반적인 범주형 변수로 결측값을 처리하는 방법도 있지만, naniar 팩키지 geom_miss_point()와 같은 함수를 사용하여 시각화하는 것도 가능하게 되었다.

단변량

airquality %>% 
  bind_shadow() %>%
    ggplot(aes(x = Wind, color = Ozone_NA)) +
      geom_density()

다변량

airquality %>% 
  bind_shadow() %>%
    ggplot(aes(x = Temp, y = Wind, color = Ozone_NA)) +
      geom_point()

앞서 잠재행렬 데이터프레임을 바탕으로 단변량, 다변량 시각화를 결측값 여부에 따라 진행하는 방법을 살펴봤다. 결측값에 대해 특별히 시각화를 지원하는 geom_miss_point() 함수를 사용하면 더 많은 정보를 시각적으로 확인할 수 있다. 특히, facet과 함께 섞어 사용할 경우 보지 못하고 놓칠 수 있는 점에 대해서도 식별이 가능하다.

geom_miss_g <- ggplot(airquality, aes(x = Ozone, y = Solar.R)) +
    geom_miss_point() +
    theme(legend.position = "none")

geom_miss_facet_g <- ggplot(airquality, aes(x = Wind, y = Ozone)) +
    geom_miss_point() +
    facet_wrap(~Month)

cowplot::plot_grid(geom_miss_g, geom_miss_facet_g)

8 채워넣기 효과 분석

결측값이 존재하는 경우와 존재하지 않는 경우 시각화를 통해 효과를 확인했다. 다음 단계로 결측값을 채워넣은 경우 효과를 파악하는 방법을 살펴보자.

8.1 impute_* 계열

naniar 팩키지 impute_* 계열 함수를 사용해서 결측값을 채워넣을 수 있다.

missing_vv <- c(1,3,5,NA,9,11,13)

impute_mean(missing_vv)
[1]  1  3  5  7  9 11 13
impute_below(missing_vv)
[1]  1.0000000  3.0000000  5.0000000 -0.4334966  9.0000000 11.0000000 13.0000000
impute_median(missing_vv)
[1]  1  3  5  7  9 11 13

bind_shadow() 함수를 사용해서 잠재행렬로 만든 후에 impute_mean_all() 함수를 사용해서 결측값을 모두 평균으로 채워넣고 이를 시각화한다. 모든 결측값이 평균값으로 대체된 것이 확인된다.

airquality %>%
    bind_shadow(only_miss = TRUE) %>%
    # add_label_missings() %>% 
    impute_mean_all() %>% 
    ggplot(aes(x = Ozone, fill = Ozone_NA)) +
      geom_histogram()

8.2 다변량 변수 결측값 채워넣기

변수 하나가 아니라 다변량인 경우 shadow_long() 함수를 통해서 결측값이 있는 변수만 다수 추려 이에 대해 결측값을 채워넣은 효과를 시각적으로 파악이 가능하다.

airquality %>%
    bind_shadow(only_miss = TRUE) %>%
    impute_mean_all() %>% 
    shadow_long(Ozone, Solar.R) %>% 
    ggplot(aes(x = value, fill = value_NA)) +
      geom_density(alpha = 0.3) +
      facet_wrap(~variable)

8.3 예측모형으로 채워넣기

simpuation 팩키지의 다양한 결측값 채워넣기 기능을 사용해서 결측값을 채워넣는 것도 가능하다.

naniar 팩키지 내장된 airquality 데이터셋에는 Ozone, Solar.R 두변수에 결측값이 내장되어 있다. 이를 simpuation 팩키지 RF 모형으로 결측값을 채워넣는다.

library(simputation)

miss_rf_df <- airquality %>%
    bind_shadow(only_miss = TRUE) %>%
    add_label_shadow() %>%
    impute_rf(Solar.R ~ Wind + Temp + Month) %>%
    impute_rf(Ozone ~ Wind + Temp + Month)

miss_mean_df <- airquality %>%
    bind_shadow(only_miss = TRUE) %>%
    impute_mean_all() %>% 
    add_label_shadow()

miss_rf_g <- miss_rf_df %>% 
    ggplot(aes(x=Solar.R, y= Ozone, color = any_missing)) +
        geom_point()

miss_mean_g <- miss_mean_df %>% 
    ggplot(aes(x=Solar.R, y= Ozone, color = any_missing)) +
        geom_point()

cowplot::plot_grid(miss_rf_g, miss_mean_g)