-
결측값(Missing Value) 처리하기R 2018. 11. 7. 10:10
데이터 과학
결측데이터
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월부터 트럼프 승리 예측
- 2016년 미국 대선에서
결측값 구분 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 FALSENA 갯수와 비율
> n_miss(missing_v)[1] 5> prop_miss(missing_v)[1] 0.52.2 결측값 연산
컴퓨터 과학에서 일반적으로 결측값(missing value)은 존재하지 않는 값(null)으로 컴퓨터에서 표현된다. 하지만 데이터 과학에서는 결측값이 무응답 혹은 단순히 자리만 차지하는
NA, 을 표현하는NaN(Not a Number), 무한을 표현하는inf가 있다.R에
NA와 같이 미리 예약된 예약어가 몇개 있다.NaN,NULL,Inf가 대표적인다.NULL: 값이 없다로 즉, 값이 존재하지 않는다는 의미로 사용NaN: Not a Number로 처럼 수학적으로 정의가 되지않는다는 의미.Inf: Infinite의 약자로 무한대를 의미하는데, 수학적으로 은 무한대가 된다.
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에 따라 값이 달라지고,NA와NaN를 연산하면 연산순서에 따라 결과가 달라진다.> NA | TRUE[1] TRUE> FALSE | NA[1] NA> NA + NaN[1] NA> NaN + NA[1] NaN3 결측값 식별과 현황파악 4 5
결측데이터를 처리하기 전에 결측데이터를 처리하는 프로세스는 다음과 같다.
- 결측값을 식별한다. 원본데이터에서 다양한 형태로 결측정보가 표현되어 있으니 우선 현황 파악이 먼저다.
- 파악된 현황정보를 바탕으로 결측값을 컴퓨터가 처리가능한 형태로 부호화한다.
- 마지막으로 파악된 결측정보와 적절히 인코딩되어 컴퓨터에 저장된 결측값을 자료형에 맞춰 알고리즘을 적용하여 결측값을 처리한다.
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"] <- NA3.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 0.000000 7.142857 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.000000naniar팩키지 함수> 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 129 2 14.3 2 168 2 14.3 3 293 2 14.3 4 479 2 14.3 5 496 2 14.3 6 8 1 7.14 7 25 1 7.14 8 26 1 7.14 9 28 1 7.14 10 30 1 7.14 # ... with 496 more rows3.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 431 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 45 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 25 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 5 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 24.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> 865 결측값 처리 사례 - 옛날 방식
결측값 처리 전략은 다음과 같은 4가지 전략이 존재한다.
- 관측점 제거
- 데이터가 상당히 많은 경우, 동시에 모집단을 대표하는데 무리가 없는 경우
na.action=na.omit설정을 적용한다. lm(medv ~ ptratio + rad, data=BostonHousing, na.action=na.omit)
- 데이터가 상당히 많은 경우, 동시에 모집단을 대표하는데 무리가 없는 경우
- 변수 제거
- 특정 변수에 상당한 값이 결측값인 경우, 다른 변수가 결측값을 많이 갖는 변수에 상응하는 정보량을 갖는 경우 제거한다.
- 평균/중위수/최빈값으로 대체(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)- 예측값으로 대체한다.
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 1.03617849 2.45337279 1.56632461 0.06567957knn대체 기법은 범주형 자료의 경우 적용에 한계가 있다. 이런 경우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] FALSE5.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.8524912 5.3668569 2.3166478 0.1104378> regr.eval(actuals_ptratio, ptratios_median_pred)mae mse rmse mape 1.7900000 5.7798000 2.4041215 0.1102176두번째
knn기법을 활용한 경우mape가 줄어든 것이 확인된다.> ptratios_knn_pred <- knnOutput[is.na(BostonHousing$ptratio), "ptratio"] > regr.eval(actuals_ptratio, ptratios_knn_pred)mae mse rmse mape 0.116221431 0.338332058 0.581663183 0.009223923세번째
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.2145987 6.8920295 2.6252675 0.1310058네번째는
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.4375000 9.1027500 3.0170764 0.14661765.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.82두번째로
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.38마지막으로
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.246 명시적, 묵시적 결측값
결측값이 명시적으로 나타난 경우와 암묵적으로 나타난 경우로 나눌 수 있다. 명시적, 묵시적 결측값 처리를 위한 사례로
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 9637 결측값 효과 분석
결측값 효과 분석을 위해서 결측값이 있는 상황과 없는 상황에 대한 명확한 구별이 필요하다.
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 !NA7.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.37.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 [7] 13.0000000> impute_median(missing_vv)[1] 1 3 5 7 9 11 13bind_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)
[출처] https://statkclee.github.io/data-science/ds-missing.html?fbclid=IwAR0DXSHM106j0-1M1pC51iuyDyaI-XryZlfYjXcT_93bltNOVc23pad1GZc
'R' 카테고리의 다른 글
dplyr로 DB 데이터를 다룰 때 여러 테이블을 union 해보자 (0) 2018.11.16 수학적 주석 넣기(Mathematical Annotation in R) (0) 2018.11.07 Frequency Table : 빈도표 (0) 2018.11.07 pu (0) 2018.11.05 정규표현식 : escape sequence (0) 2018.10.27 - 1916년 대선: 우편 설문조사를 바탕으로 민주당 우드로 윌슨 후보 당선 예측한