데이터 과학

결측데이터 - Base R 시각화

1 결측값 문제점

결측값이 존재하는 경우 적절한 결측값 처리를 하지 않고 모형을 만들고 시각화 등 후속작업을 이어나갈 경우 추론한 결론은 물론이고 전혀 다른 분석과 모형을 개발하는 결과를 초래할 수 있다. (Van Buuren 2018) NHANES “US National Health and Nutrition Examination Study”에서 도출된 데이터를 바탕으로 사례로 살펴보자. 먼저, 각 변수별로 결측값을 살펴보자. 이를 위해서 is.na() 함수로 전체 데이터프레임의 결측값을 판별하는데 출력값이 행렬(matrix)이라 이를 티블 데이터프레임으로 변환시킨 후에 신규로 dplyr 팩키지 across() 함수를 사용해서 결측값을 계산한다.

library(tidyverse)
library(NHANES)
library(broom)

set.seed(777)

nhanes <- NHANES %>% 
  janitor::clean_names() %>% 
  select(poverty, age, height, weight, tot_chol, education, gender, race1,  home_rooms, bmi) %>% 
  sample_frac(0.1)

nhanes %>% 
  is.na() %>% as_tibble() %>% 
  summarise(across(everything(), sum))

# A tibble: 1 x 10
  poverty   age height weight tot_chol education gender race1 home_rooms   bmi
    <int> <int>  <int>  <int>    <int>     <int>  <int> <int>      <int> <int>
1      83     0     42      6      165       282      0     0          1    42

lm() 선형회귀모형에서 결측치에 대해 어떤 작업도 하지 않았기 때문에 실제 base_lm 모형과 edu_lm 모형에서 사용된 회귀모형이 다른 데이터셋을 가지고 학습된 것을 확인할 수 있다.

기본 회귀모형

base_lm <- lm(poverty ~ age + weight, data = nhanes)

base_lm_summary <- summary(base_lm)

base_lm_summary$na.action %>% length

[1] 87

교육을 고려한 모형

edu_lm <- lm(poverty ~ age + weight + education, data = nhanes)

edu_lm_summary <- summary(edu_lm)

edu_lm_summary$na.action %>% length

[1] 345

2 R 팩키지

• 시각화: VIM
• 모형: simputation
• 나무모형: missForest

3 결측값 분류

결측값은 다음 세가지 경우로 구분되고 분류된 결측값 범주에 따라 적절한 조치를 취한다.

• MCAR(Missing Completely At Random): 모든 정보가 데이터에 담겨있어 결측값이 문제가 되지 않는 경우.
  • 결측값이 \(X_{mis}\), \(X_{obs}\) 모두 영향을 받지 않는 경우
  • 데이터셋의 결측값의 위치가 다른 어떤 데이터에도 의존하지 않고 순전히 랜덤인 경우.
  • 기상관측에 사용되는 온도 풍향 센서가 기계적인 고장으로 데이터를 보내지 못하는 경우.
• MAR(Missing At Random): 더 이해하기 쉬운 명칭은 Missing Conditionally at Random으로 결측조건이 다른 변수에 따라 조건부로 발생되는 경우. 결측값이 관측된 데이터가 아닌 관측되지 않는 데이터에 따라 결정.
  • 결측값이 \(X_{mis}\)에 영향을 받지는 않으나 다른 \(X_{obs}\)에 영향을 받는 경우
  • 예를 들어, 남성이 여성보다 더 솔찍하게 나이, 몸무게 등을 밝힐 듯 싶다.
  • 기상관측에 사용되는 온도 풍향 센서가 점검을 위해 동작을 멈출 때, 대부분 주말에는 작업자가 휴일이라 결측값은 특별한 사유가 아니면 주중에 정기점검으로 발생된다.
• MNAR(Missing Not At Random): 결측값이 무작위가 아니라서 주도면밀한 추가 조사가 필요한 경우.
  • 결측값이 \(X_{mis}\)에 영향을 받는 경우.
  • 예를 들어, 작업장에서 결측된 사람들 대부분은 아마도 몸이 아파 조사에서 결측되고, 교육수준이 낮은 사람이 교육조사에 결측될 듯 싶다.
  • 기상관측에 사용되는 온도 풍향 센서의 경우 온도가 매우 낮거나 매우 높은 경우 데이터가 수집되지 않는 경우가 발생된다. 이런 경우는 온도 센서 그 값 자체문제로 결측이 발생된다.

만약 결측값을 관측점별로 삭제를 하게 되면 다음과 같은 문제가 발생된다.

• MCAR: 정보 손실
• MAR, MNAR: 제거된 관측점에 기반하여 모형이나 시각화를 하게 되면 편향이 발생
• 대부분의 치환방법은 MAR을 가정하고 있어 이를 사전에 탐지하는 것이 중요!!!

4 결측값 시각화

결측값을 시각화하는 것이 통계적 검증을 수행할 경우 p-값 해킹하는 여러 문제점도 있어 시각적으로 빠르게 확인하는 것도 큰 도움이 된다.

• aggr(): 데이터프레임 전체 현황을 파악
• spineMiss()
• mosaicMiss()

4.1 DF 수준 시각화

aggr() 함수를 사용해서 데이터프레임 전체 현황을 파악하는데 aggr() 함수를 사용해서 각 변수별 결측값의 분포와 변수를 조합했을 때 발생된 결측값 비율을 동시에 확인한다.

library(VIM)

nhanes_aggr <- nhanes  %>%  
  aggr(combined = FALSE, numbers = TRUE, col = "skyblue")

summary(nhanes_aggr)

 Missings per variable: 
   Variable Count
    poverty    83
        age     0
     height    42
     weight     6
   tot_chol   165
  education   282
     gender     0
      race1     0
 home_rooms     1
        bmi    42

 Missings in combinations of variables: 
        Combinations Count Percent
 0:0:0:0:0:0:0:0:0:0   617    61.7
 0:0:0:0:0:1:0:0:0:0   151    15.1
 0:0:0:0:1:0:0:0:0:0    38     3.8
 0:0:0:0:1:1:0:0:0:0    74     7.4
 0:0:1:0:1:1:0:0:0:1    33     3.3
 0:0:1:1:0:0:0:0:0:1     2     0.2
 0:0:1:1:0:1:0:0:0:1     1     0.1
 0:0:1:1:1:1:0:0:0:1     1     0.1
 1:0:0:0:0:0:0:0:0:0    47     4.7
 1:0:0:0:0:0:0:0:1:0     1     0.1
 1:0:0:0:0:1:0:0:0:0    13     1.3
 1:0:0:0:1:0:0:0:0:0    11     1.1
 1:0:0:0:1:1:0:0:0:0     6     0.6
 1:0:1:0:0:0:0:0:0:1     1     0.1
 1:0:1:0:1:1:0:0:0:1     2     0.2
 1:0:1:1:0:0:0:0:0:1     1     0.1
 1:0:1:1:0:1:0:0:0:1     1     0.1

4.2 변수 수준 시각화

한단계 더 들어가서 변수별 결측값 함유정도를 파악할 수 있는데 이런 경우 spineMiss() 함수를 사용해서 시각화한다.

nhanes %>% 
  select(education, tot_chol) %>% 
  as.data.frame() %>% 
  VIM::spineMiss()

Click in in the left margin to switch to the previous variable or in the right margin to switch to the next variable.
To regain use of the VIM GUI and the R console, click anywhere else in the graphics window.

4.3 mosaicMiss 모자이크 플롯

두 범주형 변수(gender, education)를 고려한 상태에서 연속형 변수(tot_chol)에 대한 결측값 정보를 파악할 수 있다.

nhanes %>% 
  VIM::mosaicMiss(highlight = "tot_chol", plotvars = c("gender", "education"))

4.4 histMiss 모자이크 플롯

범주형 변수(education)를 고려한 상태에서 연속형 변수(tot_chol)에 대한 결측값 정보를 파악할 수 있다.

nhanes %>% 
  select(tot_chol, education) %>% 
  as.data.frame() %>% 
  VIM::histMiss(only.miss = FALSE)

Click in in the left margin to switch to the previous variable or in the right margin to switch to the next variable.
To regain use of the VIM GUI and the R console, click anywhere else in the graphics window.

4.5 bartMiss 모자이크 플롯

범주형 변수(education)를 고려한 상태에서 연속형 변수(tot_chol)에 대한 결측값 정보를 파악할 수 있다.

nhanes %>% 
  select(education, tot_chol) %>% 
  as.data.frame() %>% 
  VIM::barMiss(only.miss = FALSE)

Click in in the left margin to switch to the previous variable or in the right margin to switch to the next variable.
To regain use of the VIM GUI and the R console, click anywhere else in the graphics window.

5 MAR - t-검정

특정 변수에 따라 결측값이 패턴을 가지는지 여부를 파악하기 위해서 t-검정을 시도해본다. 이를 위해서 특정변수 콜레스테롤값(tot_chol)의 결측여부를 표식하는 가변수(dummy variable)를 생성하고 나서 이를 남성 여성 벡터로 만들어 t.test() 함수에 넣어 t-검정작업을 수행한다.

Base R 문법

# 콜레스테롤 결측값 표식 
nhanes <- nhanes %>% 
  mutate(chol_dummy = ifelse(is.na(tot_chol), TRUE, FALSE))

# 남성 
nhanes_male <- nhanes %>% 
  filter(gender == "male") %>% 
  pull(chol_dummy)
  
# 여성
nhanes_female <- nhanes %>% 
  filter(gender == "female") %>% 
  pull(chol_dummy)

# t-검정
t.test(nhanes_male, nhanes_female)

    Welch Two Sample t-test

data:  nhanes_male and nhanes_female
t = -0.90232, df = 983.92, p-value = 0.3671
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.06744744  0.02495832
sample estimates:
mean of x mean of y 
0.1547389 0.1759834 

dplyr 티블 방식

# 콜레스테롤 결측값 표식 
nhanes_ttest_df <- nhanes %>% 
  mutate(chol_dummy = ifelse(is.na(tot_chol), TRUE, FALSE)) %>% 
  select(gender, chol_dummy) %>% 
  group_by(gender) %>% 
  summarise(na_value = list(chol_dummy)) %>% 
  spread(gender, na_value)

nhanes_ttest_df %>% 
  mutate(ttest = list( t.test(unlist(female), unlist(male)))) %>% 
  mutate(p_value = map_dbl(ttest, "p.value"),
         t_stat = map_dbl(ttest, "statistic"))

# A tibble: 1 x 5
  female      male        ttest   p_value t_stat
  <list>      <list>      <list>    <dbl>  <dbl>
1 <lgl [483]> <lgl [517]> <htest>   0.367  0.902

6 결측값 치환

결측값을 치환(imputation)한다는 것은 나름 합리적 근거를 가지고 결측값을 유의미하게 채워넣는 것을 의미하는 것으로 결측값 치환은 크게 두가지 방식으로 나눌 수 있다.

• 기부 방식 (donor-based imputation): 다른 관측점을 활용해서 결측값을 합리적으로 채워넣음
  • 평균, 핫덱(hot-deck), kNN 등
• 모형 방식 (model-based imputation): 통계 혹은 기계학습 모형을 활용하여 결측값을 합리적으로 예측하여 채워넣음.
  • GLM, CART 등

6.1 기부방식: 평균/중위수

평균을 기부받아 결측값을 채워넣게 되면 다음과 같이 ifelse() 함수를 사용해서 가능하다. 결측값인 경우 각 변수별로 결측값을 빼고 평균을 계산해서 치환한다. 이를 시각화하는데 marginplot()을 사용하여 평균이나 중위수를 치환하는데 tot_chol 변수 관점에서 보면 특정 값으로 치환되고, 마찬가지로 bmi 변수 관점에서 봐도 동일하게 평균이나 중위수로 치환된다.

nhanes_mean_imp <- nhanes %>% 
  mutate(bmi_imp = ifelse(is.na(bmi), TRUE, FALSE),
         tot_chol_imp = ifelse(is.na(tot_chol), TRUE, FALSE)) %>% 
  mutate(bmi = ifelse(is.na(bmi), mean(bmi, na.rm = TRUE), bmi),
         tot_chol = ifelse(is.na(tot_chol), mean(tot_chol, na.rm = TRUE), tot_chol)) %>% 
  select(matches("bmi|tot_chol"))

nhanes_mean_imp %>% 
  as.data.frame() %>% 
  marginplot(delimiter = "imp")

6.2 기부방식: 핫덱(hot-deck)

핫덱(hot-deck)은 가장 마지막 관측값을 활용해서 모든 결측값을 치환하는 기법이다. 결측값 패턴이 MCAR일 경우 유의미한 결측값 치환기법이 되고 단순한 알고리즘으로 앞선 평균/중위스 결측값 치환방법과 비교하여 동일한 값을 넣는 것이 아니라 MCAR 패턴에서 결측값과 연관된 변수를 찾아 핫덱 방법을 접목시킬 경우 나름 좋은 효과를 얻을 수 있다.

nhanes <- nhanes %>% 
  mutate(weight_imp = ifelse(is.na(weight), TRUE, FALSE),
         height_imp = ifelse(is.na(height), TRUE, FALSE))

nhanes_hotdeck <- hotdeck(nhanes, 
        variable = "weight",
        ord_var = "height")

nhanes_hotdeck %>% 
  select(height, weight, weight_imp, height_imp) %>% 
  marginplot(delimiter = "imp")

6.3 기부방식: kNN

k개 인접 이웃 관측점을 활용하여 가장 거리가 가까운 것으로 파악되는 관측점들을 결합시켜 결측값을 채워넣는다. 이를 위해서 몇가지 중요한 작업을 사전 선행하는데 다음을 들 수 있다.

1. 변수별로 관측점갯수를 파악하여 결측값이 적은 변수가 우선 반영되도록 작업
2. weightDist = TRUE로 설정하여 거리에 역수를 취해 가중평균으로 근접 관측점을 총합계산한다.

num_of_missing_var <- nhanes %>% 
  select(-matches("imp")) %>% 
  is.na() %>% 
  colSums() %>% 
  sort() %>% 
  names()

nhanes_knn <- nhanes %>% 
  select(num_of_missing_var) %>% 
  kNN(k=5)

nhanes_knn %>% 
  select(height, weight, weight_imp, height_imp) %>% 
  marginplot(delimiter = "imp")

7 모형방식

관측점 대신 변수를 활용하여 모형을 만들고 모형 예측값을 결측값에 치환하는 방식이다. 앞서 살펴본바와 같이 결측값이 있는 경우 선형회귀모형을 통한 치환은 여전히 결측값을 담고 있다.

library(simputation)

nhanes <- nhanes %>% 
  select(-matches("imp")) %>% 
  as.data.frame()

nhanes_lm_imp <- impute_lm(nhanes, height + weight ~ .)

회귀로 치환전

nhanes %>% is.na() %>% colSums() %>% as.data.frame()

             .
poverty     83
age          0
height      42
weight       6
tot_chol   165
education  282
gender       0
race1        0
home_rooms   1
bmi         42
chol_dummy   0

회귀로 치환후

nhanes_lm_imp %>% is.na() %>% colSums() %>% as.data.frame()

             .
poverty     83
age          0
height      42
weight       6
tot_chol   165
education  282
gender       0
race1        0
home_rooms   1
bmi         42
chol_dummy   0

7.1 모형원리

회귀모형을 통해 치환하는 방식은 A변수 결측값을 회귀모형으로 추정하여 채워넣고 B변수 결측값을 회귀모형으로 추정하여 채워넣고 더이상 변동이 없거나 매우 적은 경우 계산을 멈추는 방식이다.

# 핫덱으로 초기화
nhanes_imp <- hotdeck(nhanes)

# 결측값 마스크
missing_tot_chol <- nhanes_imp$tot_chol_imp

# Step 1: 적절한 횟수만큼 반복 !!!
# for (i in 1:3) {
#   # 원래 결측값을 NA 치환 후 회귀모형 적합
#   nhanes_imp$tot_chol[missing_tot_chol] <- NA
#   nhanes_imp <- impute_lm(nhanes_imp, tot_chol ~ age + gender + bmi)
# }

calc_improvement <- function(prev, curr) {
  mean(abs(curr - prev)/prev, na.rm=TRUE)
}

# Step 2: 치환값 추적 !!!
# for (i in 1:3) {
#   # 기준값
#   prev_nhanes <- nhanes_imp
#   # 원래 결측값을 NA 치환 후 회귀모형 적합
#   nhanes_imp$tot_chol[missing_tot_chol] <- NA
#   nhanes_imp <- impute_lm(nhanes_imp, tot_chol ~ age + gender + bmi)
#   # 전후 향상도 출력
#   cat("processing: ", calc_improvement(nhanes_imp$tot_chol, prev_nhanes$tot_chol), "\n")
# }

# Step 3: 치환값 자동추적 !!!
eps <- 0

while(eps < 0.01) {
  # 기준값
  prev_nhanes <- nhanes_imp
  # 원래 결측값을 NA 치환 후 회귀모형 적합
  nhanes_imp$tot_chol[missing_tot_chol] <- NA
  nhanes_imp <- impute_lm(nhanes_imp, tot_chol ~ age + gender + bmi)
  # 전후 향상도 출력
  eps <- calc_improvement(nhanes_imp$tot_chol, prev_nhanes$tot_chol)
  cat("processing: ", eps, "\n")
}

processing:  0.03408077 

7.2 나무모형

missForest 팩키지를 활용하여 nhanes 데이터셋의 결측값을 채워넣는다.

library(missForest)

nhanes_cart <- missForest(nhanes) 

  missForest iteration 1 in progress...done!
  missForest iteration 2 in progress...done!
  missForest iteration 3 in progress...done!
  missForest iteration 4 in progress...done!
  missForest iteration 5 in progress...done!
  missForest iteration 6 in progress...done!
  missForest iteration 7 in progress...done!

nhanes_cart_imp <- nhanes_cart$ximp 

nhanes_cart_imp %>% is.na() %>% colSums() %>% as.data.frame()

           .
poverty    0
age        0
height     0
weight     0
tot_chol   0
education  0
gender     0
race1      0
home_rooms 0
bmi        0
chol_dummy 0

성능평가는 NRMSE, PFC 값이 0에 가까우면 무난히 결측값이 채워진 것을 확인할 수 있다.

nhanes_cart$OOBerror

     NRMSE        PFC 
0.05676222 0.19498607 

missForest() 함수에 나무모형 초모수(hyper-parameter)를 조정해서 모형 적합도와 조정과 variablewise = TRUE을 통해 변수별 적합도도 함께 판단할 수 있다.

nhanes_custom_cart <- missForest(nhanes, variablewise = TRUE, mtry =3, ntree = 100) 

  missForest iteration 1 in progress...done!
  missForest iteration 2 in progress...done!
  missForest iteration 3 in progress...done!
  missForest iteration 4 in progress...done!
  missForest iteration 5 in progress...done!
  missForest iteration 6 in progress...done!

nhanes_colnames <- nhanes_cart_imp %>% is.na() %>% colSums() %>% names

names(nhanes_custom_cart$OOBerror) <- nhanes_colnames
nhanes_custom_cart$OOBerror

   poverty        age     height     weight   tot_chol  education     gender 
 1.6669358  0.0000000 23.5983990 34.2195191  0.8853788  0.5905292  0.0000000 
     race1 home_rooms        bmi chol_dummy 
 0.0000000  3.9816343  4.8679875  0.0000000 

8 mice와 부츠트래핑

mice - with - pool 작업흐름을 갖는다.

library(mice)

nhanes_df <- nhanes %>% 
  select(tot_chol, poverty, height, gender)

# 변수 선정
pred_mat <- quickpred(nhanes_df, mincor = 0.3)

nhanes_mice <- mice(nhanes_df, m = 10, defaultMethod = "pmm", predictorMatrix = pred_mat)

 iter imp variable
  1   1  tot_chol  poverty  height
  1   2  tot_chol  poverty  height
  1   3  tot_chol  poverty  height
  1   4  tot_chol  poverty  height
  1   5  tot_chol  poverty  height
  1   6  tot_chol  poverty  height
  1   7  tot_chol  poverty  height
  1   8  tot_chol  poverty  height
  1   9  tot_chol  poverty  height
  1   10  tot_chol  poverty  height
  2   1  tot_chol  poverty  height
  2   2  tot_chol  poverty  height
  2   3  tot_chol  poverty  height
  2   4  tot_chol  poverty  height
  2   5  tot_chol  poverty  height
  2   6  tot_chol  poverty  height
  2   7  tot_chol  poverty  height
  2   8  tot_chol  poverty  height
  2   9  tot_chol  poverty  height
  2   10  tot_chol  poverty  height
  3   1  tot_chol  poverty  height
  3   2  tot_chol  poverty  height
  3   3  tot_chol  poverty  height
  3   4  tot_chol  poverty  height
  3   5  tot_chol  poverty  height
  3   6  tot_chol  poverty  height
  3   7  tot_chol  poverty  height
  3   8  tot_chol  poverty  height
  3   9  tot_chol  poverty  height
  3   10  tot_chol  poverty  height
  4   1  tot_chol  poverty  height
  4   2  tot_chol  poverty  height
  4   3  tot_chol  poverty  height
  4   4  tot_chol  poverty  height
  4   5  tot_chol  poverty  height
  4   6  tot_chol  poverty  height
  4   7  tot_chol  poverty  height
  4   8  tot_chol  poverty  height
  4   9  tot_chol  poverty  height
  4   10  tot_chol  poverty  height
  5   1  tot_chol  poverty  height
  5   2  tot_chol  poverty  height
  5   3  tot_chol  poverty  height
  5   4  tot_chol  poverty  height
  5   5  tot_chol  poverty  height
  5   6  tot_chol  poverty  height
  5   7  tot_chol  poverty  height
  5   8  tot_chol  poverty  height
  5   9  tot_chol  poverty  height
  5   10  tot_chol  poverty  height

nhanes_mice_lm <- with(nhanes_mice, lm(tot_chol ~ poverty + height))

nhanes_mice_pool <- pool(nhanes_mice_lm)

summary(nhanes_mice_pool, conf.int = TRUE, conf.level = 0.95)

         term   estimate   std.error statistic        df      p.value
1 (Intercept) 3.56043333 0.347413555 10.248401  51.88538 4.529710e-14
2     poverty 0.03467002 0.023478689  1.476659 116.34857 1.424685e-01
3      height 0.00739669 0.002092375  3.535070  56.43826 8.220906e-04
         2.5 %     97.5 %
1  2.863260371 4.25760629
2 -0.011831009 0.08117106
3  0.003205878 0.01158750

시각화

stripplot(nhanes_mice,
          tot_chol ~ height | .imp,
          pch = 20, cex = 2)

참고문헌

Van Buuren, Stef. 2018. Flexible Imputation of Missing Data. CRC press.