Disclaimer

이 글은 언제나 ’진행형’입니다. 문제가 생기면 항상 수정되거나 폐기됩니다. 대신 수정될 경우 변경 역사는 아래와 같이 간단히 남길 예정입니다. 각종 문의는 로 메일 주시라.

  1. 20180611 오타 수정
  2. 20180612 샘플링 분포 그림 추가

최근의 논란

최근 최저임금 인상의 효과를 두고 몇가지 논란이 있었다. 일단 두 가지를 먼저 언급해야 겠다.

  1. (사회) 실험이 아니기 때문에 최저 임금 인상의 효과를 인과관계라는 맥락에서 파악하는 것은 불가능하다.
  2. 적용할 수 있는 계량경제학적 방법이 없는 것은 아니지만, 아직 관련한 본격적인 논의를 보지는 못했다.

이 와중에 조선비즈에 정부의 발표를 반박하는 기사가 실렸고, 이에 대해서 블로고스피어 상의 비판이 있었다. 나는 이 논란에 개입하지 않을 것이다. 아울러 무슨 대단한 학술적인 의도나 목적을 지니고 있지 않다. 업계나 학계의 상식에 어긋나거나 황당한 내용이 있다면 언제든 지도편달 바란다. 달게 받고 수정하겠다. 다만 한가지 아쉬움만 지적하자. 이 논란의 시작이 청와대의 90% 소득 증가 발언이다. 정부 혹은 국책 연구기관이 직접 보다 상세한 해설을 해주면 어떨까 싶다.

글의 목적

  1. 데이터에 의지해 서로 다른 논의를 확인헤보자.
  2. 불필요한 사항으로 다투지 말자.

분석의 전제

이 분석은 통계청에서 발표하는 “가계동향조사”에 바탕을 둔다. 이 자료의 여러 가지 한계는 잠시 접어두자. 그럼에도 불구하고 그나마 빠르게 확인할 수 있는 자료가 이 녀석이니 쓰는 것이겠지. 아래 두 가지를 가정하자.

  1. 가계동향조사를 믿는다.
  2. 조사의 가중치를 믿는다.

가중치를 어떻게 써먹을까?

가중치는 결국 해당 가구가 모집단에서 차지하는 비중을 나타낸다. 2017년 판 “이용자 가이드”에 보면 다음과 같이 되어 있다.

가계동향조사의 가중치는 설계가중치와 무응답 가중치를 적용하고 사후 가중한 가중치를 가구별로 부여하며, 가중치는 추출률의 역수로서 각 가구가 대표하는 가구 수를 나타낸다.

가중치에 관해서 교과서에 나오는 정도의 내용을 잘 반영했다는 이야기다. 대안적인 가중치를 구하고 싶어도 공개되지 않은 추가적인 자료가 필요하다. 그러니 일단 믿도록 하자. 가중치를 써 먹을 세가지 정도의 방법을 생각해보았다.

  1. 100 분위 선정에서 가중치를 무시하고, 대신 100 분위 블록 내 소득 계산에 활용한다.
  2. 가중치 만큼 해당 샘플을 반복해서 늘린 뒤, 확대된 샘플로 100분위 블록 계산한다.
  3. 가중치에 따른 확률에 따라서 100 개씩 뽑는 샘플링을 충분히 반복한다.

첫번째 방법은 가중치를 무시하고 각 연도 샘플에서 100 분위를 구한다. 대신 해당 분위 내에서 평균 소득 혹은 중위 소득을 구할 때 해당 가중치를 활용해서 가중 평균, 가중 중간값을 구한다. 뭐 단순한 방법이지만 일종의 기저(baseline)로서의 해볼만한 방법이다. 둘째 방법은 보다 직관적이다. 가중치는 샘플에 속한 각 가구가 대표하는 가구의 수다. 그 숫자대로 뻥튀기를 한 뒤 이 녀석으로 다시 100 분위를 내는 방법이다. 하지만, 이는 해당 가구의 크기를 기계적으로 반영할 뿐 결국 이 가구들이 모집단에서 뽑혔다는 점을 담지 못한다. 나중에 보겠지만 특정 분위에서 소득의 분포가 안정적이지 않을 수 있다. 샘플을 늘리는 방법으로는 표본이 지니고 있을지 모르는 미묘한 특성을 반영하기 힘들다.

일단 1,2의 방법을 구현해서 Shiny app으로 재현 가능하게 만들어 보았다. 해당 방법이 논란의 중심이기도 해서 공을 좀 들였다. 상세한 내용은 해당 앱의 “REAME” 페이지를 참고하면 되겠다.1 대체적인 논의의 양상은 아래 세번째 방법과 비슷하기 때문에 별도의 언급은 생략하도록 하자.

반복 샘플링을 통한 각 분위 소득 분포 생성

100 분위 소득이니 속편하게 샘플 수는 100 개로 하자. 100개의 크기를 지닌 표본 추출을 충분히 많이 반복한다. 그리고 이 반복 추출된 결과에서 각 분위 별로 소득의 평균(중간값)을 구한다. 엄밀한 방법이 아닐지 모르지만 제법 그럴 듯하지 않은가. 반복 샘플링을 구현한 코드를 일부분 옮기면 아래와 같다.

vars_basic <- c("V100", "V101", "V102", "V103")

pick_x <- function(year, x, df, what_V4 = c(1,2)){
  df %.>% 
    dplyr::filter(., V1 == year, V4 %in% what_V4) %.>% 
    sample_n(., size = x, replace = T, weight = V99)
}
sim_ndraw <- function(n_draw, what_V4, vars_sel = vars_basic, df=tbl){
#  
  list(year = c(rep(2017,n_draw), rep(2018,n_draw)),  round=c(seq(1:n_draw), seq(1:n_draw))) -> itr_sch
  
  itr_sch %.>% 
    pmap_df(., 
      function(year, round) {
        pick_x(year, 100, df, what_V4) %.>% 
        mutate(., rd = unlist(round)[1]) 
      } 
    ) %.>% 
    dplyr::select(., V1, V3, rd, one_of(vars_sel))
#
}
sim_summary <- function(var, funs, df){
  
  var_exp <- enquo(var)
  
  df %.>% 
    arrange(., V1, rd, -!!var_exp) %.>% 
    group_by(., V1, rd) %.>% 
    mutate(., rank_sample = seq(1:100)) %.>% 
    group_by(., V1, rank_sample) %.>% 
    summarise(., val = funs(!!var_exp)) %.>% 
    group_by(., rank_sample) %.>%  
    mutate(., V1 = factor(V1)) %.>% 
    spread(., key = V1, value = val) %.>% 
    rename(., 
      Y17 = `2017`,
      Y18 = `2018`
      )
}

#sim_ndraw(10000, c(1,2), c(vars_basic, "V106")) -> ttl1
#sim_ndraw(10000, c(1), c(vars_basic, "V106"))   -> ttl2

1만 번 샘플링한 위의 결과를 그림으로 나타내보자.

readRDS("./data/ttl1.rds") -> ttl1 
readRDS("./data/ttl2.rds") -> ttl2 

ttl1 %.>% 
  sim_summary(V101, mean, .) %.>% 
  mutate(., 
    group = "all households" 
  )-> ttl1_1

ttl2 %.>% 
  sim_summary(V101, mean, .)  %.>% 
  mutate(., 
    group = "labor-income households" 
  )-> ttl2_1

ttl1_1 %.>% bind_rows(., ttl2_1) -> tbl_smr 
tbl_smr %.>% ggplot(.) + 
  aes(x = rank_sample, y = Y18-Y17) + 
  geom_col(color = "black") + 
  facet_grid(.~group, scales = "free_y") + 
  ggtitle("Sampling Income Mean") + 
  xlab("rank(descending order)") + 
  ylab("income difference(2018-2017)") +
  theme_economist()

tbl_smr %.>% ggplot(.) + 
  aes(x = rank_sample, y = (Y18-Y17)*100/Y17) + 
  geom_col(color = "black") + 
  facet_wrap(~group, ncol=2) + 
  ggtitle("Sampling Income Change") + 
  xlab("rank(descending order)") + 
  ylab("Change(%)") +
  theme_economist() 

경상 소득의 경우

여기서 전체가구(all households)란 노동소득 가구(labor-income households)와 노동외소득 가구(nonlabor-income households)를 모두 포괄한다. 이렇게 보면 상위 54 분위까지 그리고 하위에 속한 98, 99, 100 분위에서 소득 증가가 관찰된다. 노동소득 가구만 보면 전체 분위에서 소득의 증가를 관찰할 수 있다.

변화율

소득 규모에 따른 효과를 보기 위해서 각 구간별로 변화율을 계산해 봤다. 그림에서 보듯이 V자 형태가 나타난다. 증가 혹은 감소라는 이분법적인 틀을 벗어나 보면, 두 기간 사이에 소득 증가는 상위층 혹은 하위층일수록 많이 관찰되고 중간층에서 가장 낮게 관찰된다.

샘플링 100 분위의 분포

이렇게 뽑힌 100 분위 소득의 경우 각 분위별로 1 만개의 관찰 수를 지니고 있다. 이 관찰 수들의 분포는 어떤 모양일까?

상위 1% 분위에서 분포의 안정성이 깨졌다. 보다 자세히 보기 위해서 2018년 전체 가구 자료에서 상위 10% 이상만 한번 그려봤다.

“가계동향조사”로 극상위층의 소득을 추정하는 작업에 문제가 있다는 그간의 지적을 확인할 수 있다. 4,3,2,1 분위(상위 4% 이상)에서 분포가 안정적이지 않다. 대략 상위 5% 이상은 제거하고 비교하는 게 타당하지 않을까 싶다.

The code

github.com/anarinsk/MDIS

이 문서를 포함해 모든 소스는 리포지토리에 있다. 마음껏 가져다 쓰시라! 그런데 MDIS 사용 규약에는 어긋나서, 혹시라도 태클이 들어오면 리포를 내려야 할 수도 있겠다.

정리

  • 샘플의 30% 정도만 (경상) 소득이 늘었다는 이야기가 일부 자료에서는 유사한 형태로 재현 가능하다. 하지만 가중치를 고려하면 전체 가구에 대해서는 50~60% 정도의 소득 증가가 적절한 이야기가 아닐까 싶다. 노동소득 가구는 거의 전 분위에서 경상 소득이 증가했다. 최근 경기를 고려하면 상식적으로 타당하다.
  • 퍼센트 변화를 보면 V 형태의 특징이 관찰된다. 관련한 해석은 각자 알아서.

  1. 원래 Shiny app을 무료로 호스팅할 수 있는 서비스가 있으나 1G 램의 제약 때문에 코드가 제대로 돌지 않았다. 별도로 사설 서비스를 마련해 서비스를 구현했다. 언제 없어질 것인지는…