결측값 다루기 : 왜 평균대체가 문제인가?
Missing
결측값을 포함한 자료 분석하기
결측값을 포함한 자료를 분석할 때에는 결측값의 존재에 대해 신경을 써야 한다. 가장 중요한 점은 결측값이 발생하는 원인을 파악하는 것이다.
루빈(Rubin)은 결측값이 발생하는 원인/과정 세 가지로 분류했다.
- Missing Completely At Random(MCAR) : 결측값은 완전히 무작위로 발생한다. 다른 말로 결측값의 분포는 관찰된 값의 분포와 같다. \(p(y|\textrm{observed}) = p(y|\textrm{missing})\)
- Missing At Random(MAR) : 결측값의 발생확률은 다른 변수에 의해 영향을 받는다. 하지만 다른 변수에 조건부로 결측 확률과 실제 값은 독립이다. \(p(y|\textrm{observed}, x_{1}, x_{2}, \cdots, x_{n}) = p(y|\textrm{missing}, x_{1}, x_{2}, \cdots, x_{n})\)
- Missing Not At Random(MNAR) : 위의 둘을 제외한 나머지
데이터
library(dplyr)## Warning: package 'dplyr' was built under R version 4.1.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.3# 표본 크기 100
n <- 100
# 키 평균 170, 표준편차 15
h <- rnorm(100, 170, 15)
# 키와 체중의 관계 0.48 * 키
# 키가 주어졌을 때, 체중의 분포 N(0,7^2)
w <- 0.48 * h + rnorm(n, 0, 7)
# 체중의 모평균
w_pop <- 170*0.48 # 81.6
h_pop <- 170
# 키가 클수록 체중의 결측 확률이 높아진다.
w_missing <- runif(n, 0, 1) < (h-min(h))/(max(h)-min(h))
dat = data.frame(h=h,
w=ifelse(w_missing, NA, w),
w_complete = w,
w_missing = w_missing)
#dat %>% gather()
ggplot(dat, aes(x=h, y=w_complete, col=factor(w_missing))) +
geom_point() +
scale_color_manual(values=c('black', 'grey'))
데이터 시각화
# 관찰값과 결측값의 비교
library(gghighlight)
library(cowplot)
ggp_complete <-
ggplot(dat, aes(x=h, y=w_complete)) +
geom_point() +
#facet_wrap(~w_missing) +
gghighlight(!w_missing) +
labs(title='Complete Data')## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.#install.packages('ggExtra')
cop1 <- ggExtra::ggMarginal(ggp_complete + theme_classic(), type='histogram')
ggp <-
ggplot(dat, aes(x=h, y=w)) +
geom_point() +
labs(title='Missing Removed')
# ggp <-
# ggplot(dat, aes(x=h, y=w_complete, col=w_missing)) +
# geom_point() +
# scale_color_manual(values=c("black", "white")) +
# labs(title='Missing Removed') +
# guides(col="none")
cop2 <- ggExtra::ggMarginal(ggp + theme_classic(), type='histogram')## Warning: Removed 43 rows containing missing values (geom_point).## Warning: Removed 43 rows containing missing values (geom_point).plot_grid(cop1, cop2)
데이터 분석
결측값 제외
체중이 모두 관측된 경우의 평균과 결측값을 제외한 경우의 체중 평균을 비교해보면, 결측값을 제외하면 체중 평균에서 편향이 관찰된다. 다시 말해, 결측값을 제외하고 계산된 평균은 실제 평균보다 대부분 작다.
## average weight?
mean(dat$w_complete)## [1] 82.23173wmean_est = mean(dat$w, na.rm = TRUE)
wmean_est## [1] 81.26467t.test(dat$w)##
## One Sample t-test
##
## data: dat$w
## t = 59.015, df = 56, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 78.50617 84.02317
## sample estimates:
## mean of x
## 81.26467결측값 대체(missing imputation): 평균대체법(mean substitution)
자주 거론되는 가장 간단한 방법은 평균대체법이다. 결측값을 평균으로 대체하는 방법이다. 하지만 이 방법은 추정평균값에 존재하는 편향을 제거하지 못한다.
w_mean <- mean(dat$w, na.rm=TRUE)
dat$w_imputed <- ifelse(dat$w_missing, w_mean, w)
wmean_est = mean(dat$w_imputed)
wmean_est## [1] 81.26467res <- t.test(dat$w_imputed)
res##
## One Sample t-test
##
## data: dat$w_imputed
## t = 103.93, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 79.71320 82.81614
## sample estimates:
## mean of x
## 81.26467체중모평균은 81.6이지만 추정된 평균값은 81.26이며 신뢰구간은 79.71-82.82이다.
결측값 대체(missing imputation): 회귀대체법
결측 체중값을 회귀모형을 통해 추정하는 방법이다. 다음의 코드는 결측 체중값을 회귀분석으로 추정한 후, 이 값을 사용하여 체중의 평균을 추정한다.
mod <- lm(w ~ h)
w_hat <- predict(mod, dat)
#w_hat <- coef(mod) %*% rbind(1,dat$h)
dat$w_imputed <- ifelse(dat$w_missing, w_hat, w)
wmean_est = mean(dat$w_imputed)
wmean_est## [1] 82.93874res <- t.test(dat$w_imputed)
res##
## One Sample t-test
##
## data: dat$w_imputed
## t = 91.145, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 81.13318 84.74429
## sample estimates:
## mean of x
## 82.93874추정된 체중 평균은 82.94, 95%-신뢰구간은 81.13-84.74이다. 다음의 결측값이 없는 경우의 체중 평균 추정값과 신뢰구간을 비교해보자.
mean(dat$w_complete)## [1] 82.23173t.test(dat$w_complete)##
## One Sample t-test
##
## data: dat$w_complete
## t = 83.954, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 80.28822 84.17523
## sample estimates:
## mean of x
## 82.23173한 가지 주목할 점은 표준오차가 위의 회귀대체법보다 크다는 점이다. 회귀대체법은 경우에 따라 평균대체법에 보다 정확하게 모수를 추정하지만 추정 오차가 작게 계산되는 문제가 있다. 이를 해결하는 방법의 하나는 다중대체(multiple imputation) 또는 베이지안 방법을 사용하는 것이다.
결론
만약 관심의 대상이 \(Y = \beta_0 + \beta_1 x\) 의 \(\beta_1\) 인 경우 결측값을 제외하고 회귀분석을 해도 \(\beta_1\) 의 추정값에 편향이 발생하지 않는다. 하지만 \(Y\) 의 평균을 추정할 때에는 편향이 발생한다. 이를 해결하는 가장 간단한 방법은 결측값 \(y\) 를 회귀대체법으로 대체하는 것이다.
PS
물론 \(Y = \beta_0 + \beta_1 x\) 에서 \(\beta_0\) 과 \(\beta_1\) 을 추정하고 \(\mathbb{E}[X]\) 도 추정한 후, \(\mathbb{E}[Y] = \beta_0 + \beta_1 \mathbb{E}[X]\) 의 관계식과 \(\hat{\beta_0}\) 과 \(\hat{\beta_1}\) , 그리고 \(\hat{\mathbb{E}[X]}\) 의 추정오차를 고려해서 추정하는 방법도 있다.
Leave a comment