한국복지패널데이터를 분석하여 봅시다 :) part2

이번에는 직업데이터를 활용해서 다양한 분석을 해보겠습니다. 먼저 직업별 월급 차이를 알아 볼 건데요, 먼저 변수를 검토하고 전처리를 해보겠습니다. 먼저 직업을 나타낸 변수에 대해 알아보겠습니다.

 

table(welfare$code_job)

 

결과를 보시면 직업코드들을 보실 수 있습니다. 직업분류코드를 이용해서 직업명칭 변수를 만들어야 합니다. 직업분류코드 목록을 불러와서 데이터프레임형태로 만들어 줍니다.

 

list_job <- read_excel("c:\\Rdata\\Data\\Koweps_Codebook.xlsx", col_names = T, sheet = 2)
head(list_job)
dim(list_job)  #149개의 직업

 

데이터테이블은 dim으로 형태를 알아볼 수 있습니다. 보시면 2개의 열에 149개의 행이 있는데 이것은 149개의 직업이 있다고 생각할 수 있습니다. 이제 이 직업명칭테이블을 저희가 가지고 있던 welfare의 직업코드들과 합쳐야 합니다. 바로 left join으로요.

 

#code_job은 양쪽 테이블에 존재하기 때문에 code_job을 키로 두개의 데이터프레임을 합칩니다.
welfare<- left_join(welfare, list_job, id = "code_job")

welfare %>%
  filter(!is.na(code_job)) %>%
  select(code_job, job) %>%
  head(10)

 

결과를 확인하면 잘 합쳐진 것을 알 수 있습니다. 이제 이 변수를 이용해서 직업별 월급차이를 분석해보겠습니다. 변수간의 관계를 확인해보죠. 먼저 직업별 월급 평균을 구해보겠습니다.

 

job_income <- welfare %>%
  filter(!is.na(job) & !is.na(income)) %>%
  group_by(job) %>%
  summarise(mean_income = mean(income))

head(job_income)

 

한눈에 알아보기 쉽도록 내림차순으로 정렬해보고 상위 10개의 직업을 뽑아보겠습니다.

 

top10 <- job_income %>%
  arrange(desc(mean_income)) %>%
  head(10)

ggplot(data = top10, aes(x = reorder(job, mean_income), y = mean_income)) +
  geom_col() +
  coord_flip() # x축의 이름들이 길어서 보기편하게 x축과 y축을 바꿔줍니다.

 

 

그래프를 확인하시면 월급 평균 상위 10가지 직업을 보실 수 있습니다. 그렇다면 반대로 어떤 직업들이 월급 평균이 적은지도 뽑아보겠습니다. 상위 10가지 직업과의 비교를 편하게 하기 위해서 y축의 범위를 같게 설정해주겠습니다.

 

#하위 10위 추출
bottom10 <- job_income %>%
  arrange(mean_income) %>%
  head(10)

ggplot(data = bottom10, aes(x = reorder(job, -mean_income), y = mean_income)) +
  geom_col() +
  coord_flip() +
  ylim(0, 850)

이렇게 월급평균이 낮은 하위 10개의 직업도 알아보았습니다.　중간에 약사 및 한약사도 있는데요.. 보통 저희는 월급이 높은 직업 중 하나로 알고 있는데요, 저렇게 나온 이유에 대해서는 나중에 파고들어 보겠습니다.

 

다음은 성별에 따른 직업을 알아보겠습니다. 이것으로 꼭 그렇다고는 할 수 없지만, 어느 정도 남녀의 선호직업 혹은 남초, 여초 직업들에 대해서도 알 수 있겠네요. 바로 변수 간의 관계를 분석해보겠습니다.

 

#남성 직업 빈도 상위 10개
job_male <- welfare %>%
  filter(!is.na(job) & sex == "male") %>%
  group_by(job) %>%
  summarise(n = n()) %>%
  arrange(desc(n)) %>%
  head(10)

ggplot(data = job_male, aes(x= reorder(job, n), y = n)) +
  geom_col() +
  coord_flip()

job_female <- welfare %>%
  filter(!is.na(job) & sex == "female") %>%
  group_by(job) %>%
  summarise(n = n()) %>%
  arrange(desc(n)) %>%
  head(10)

ggplot(data = job_female, aes(x= reorder(job, n), y = n)) +
  geom_col() +
  coord_flip()

이렇게 남성들이 많이 가지고 있는 직업과 여성들이 많이 가지고 있는 직업에 대해서 알아보았습니다 :)

 

다음은 조금 참신하게 종교유무에 따른 이혼율을 알아보겠습니다. 언뜻 보면 두 변수가 연관이 있을까 싶은데요. 한 번 분석해서 알아보겠습니다.

 

table(welfare$religion)

 

결과를 보시면 1과 2로 나뉘어져있습니다. 1은 있음, 2는 없음인데요, 이 또한 표시해주겠습니다.

 

welfare$religion <- ifelse(welfare$religion == 1, "있음", "없음")
table(welfare$religion)

qplot(welfare$religion)

종교가 있다고 답한 사람이 8047명 없는 사람이 8617명인 것을 확인했습니다. 이제 혼인상태 변수를 확인해보겠습니다.

 

table(welfare$marriage)

 

결과를 보시면 0에서 6까지 다양한 결과가 나오는데요, 이는 0은 비해당(18세 미만), 1은 유배우, 2는 사별, 3은 이혼, 4는 별거, 5는 미혼, 6은 기타(사망 등)입니다. 이 중 저희가 필요한 것은 배우자가 있는 경우와 이혼한 경우만 필요하므로 따로 뽑아줍니다.

 

welfare$group_marriage <- ifelse(welfare$marriage == 1, "결혼",

ifelse(welfare$marriage == 3, "이혼", NA))

#필요하지 않은 집단은 NA로 두고 나중에 분석에서 제외합니다.

table(welfare$group_marriage)

 

결과를 확인해보시면 결혼은 8431명, 이혼은 712명이 있다는 것을 알 수 있습니다. 이것을 이용해서 종교 유무에 따른 이횬율 표를 만들어 보겠습니다. 먼저 종교 유무 및 결혼 상태별로 나눠 빈도를 구한 뒤 각 종교 유무 집단의 전체빈도로 나눠 소수 첫째자리까지의 비율을 구합니다.

 

religion_marriage <- welfare %>%
  filter(!is.na(group_marriage)) %>%
  group_by(religion, group_marriage) %>%
  summarise(n = n()) %>%
  mutate(tot_group = sum(n)) %>%
  mutate(pct = round(n/tot_group*100, 1))

religion_marriage

 

이혼에 해당하는 값만 추출해 이횬율 표를 만듭니다.

 

divorce <- religion_marriage %>%
  filter(group_marriage == "이혼") %>%
  select(religion, pct)

ggplot(data = divorce, aes(x= religion, y = pct)) + geom_col()

표를 보면 종교가 없는 경우의 이혼율이 1.1% 더 높은 것으로 나왔는데요, 이것을 가지고는 종교가 없는 사람들이 이혼을 덜 한다고 보는 것은 조금 성급한 것 같습니다. 만약 이 차이에 의미를 둔다면 다른 데이터들로 근거를 더욱 보강해야 되겠네요.

 

지금까지는 전체를 대상으로 종교유무에 따른 이혼율을 분석했습니다. 이번에는 종교유무에 따른 이혼율이 연령대별로 다른지를 알아보겠습니다.

 

ageg_marriage <- welfare %>%
  filter(!is.na(group_marriage)) %>%
  group_by(ageg, group_marriage) %>%
  summarise(n = n()) %>%
  mutate(tot_group = sum(n)) %>%
  mutate(pct = round(n/tot_group*100, 1))

 

초년은 사례가 부족해서 다른 연령대와 비교하기에 적합하지 않으므로 중년층과 노년층만 비교해보겠습니다.

 

ageg_divorce <- ageg_marriage %>%
  filter(ageg != "young" & group_marriage == "이혼") %>%
  select(ageg, pct)

ggplot(data = ageg_divorce, aes(x = ageg, y = pct)) + geom_col()

중년층과 노년층의 이혼율을 알 수 있습니다. 이것을 종교의 유무를 통해 알아보겠습니다.

 

ageg_religion_marriage <- welfare %>%
  filter(!is.na(group_marriage) & ageg != "young") %>%
  group_by(ageg,religion, group_marriage) %>%
  summarise(n = n()) %>%
  mutate(tot_group = sum(n)) %>%
  mutate(pct = round(n/tot_group*100, 1))

 

이혼율만 추출하여 그래프를 그립니다.

 

df_divorce <- ageg_religion_marriage %>%
  filter(group_marriage == "이혼") %>%
  select(ageg, religion, pct)

ggplot(data = df_divorce, aes(x = ageg, y = pct, fill = religion)) +
  geom_col(position = "dodge")

결과를 보시면 중년층은 종교가 없는 사람들의 이혼율이 높고, 노년층은 그 반대입니다. 이것은 나이대에 따른 차이이므로 어떠한 유의미한 결론은 도출할 수도 있을 것 같네요.

 

이제 마지막으로 지역별 연령대 비율을 알아보겠습니다. 작은 프로젝트를 할 때 많이 사용하는 기본적인 분석 중 하나입니다. 이 내용을 많은 결과들을 도출할 수 있고, 근거가 될 수 있습니다. 앞에서 사용했던 연령대 변수를 이용하여 변수 간 관계를 분석해 보겠습니다.

 

table(welfare$code_region)

결과를 보시면 1~7까지의 숫자가 보이는데요, 이것은 지역코드입니다. 1번 서울, 2번 수도권(인천/경기), 3번 부산/경남/울산, 4번 대구/경북, 5번 대전/충남, 6번 강원/충북, 7번 광주/전남/전북/제주도입니다. 이것을 이용해서 지역명 변수를 추가하죠.

 

list_region <- data.frame(code_region = c(1:7),
                          region = c("서울",
                                     "수도권(인천/경기)",
                                     "부산/경남/울산",
                                     "대구/경북",
                                     "대전/충남",
                                     "강원/충북",
                                     "광주/전남/전북/제주"
                                     ))
list_region

 

결과를 보시면 코드에 맞게 지역명이 들어간 것을 확인할 수 있습니다. 이렇게 지역 코드 목록을 만들었다면 이제 테이블에 지역명 변수를 추가하겠습니다.

 

welfare <- left_join(welfare, list_region, id = "code_region")

welfare %>%
  select(code_region, region) %>%
  head

 

이렇게 지역명 변수를 추가하고 확인해보았습니다. 이제 변수 간 관계를 분석해봅시다. 먼저 지역별 연령대 비율표를 만들겠습니다. 지역 및 연령대별로 나눠 빈도를 구하고, 각 지역의 전체 빈도로 나눠 비율을 구합니다.

 

region_ageg <-  welfare %>%
  group_by(region, ageg) %>%
  summarise(n = n()) %>%
  mutate(tot_group = sum(n)) %>%
  mutate(pct = round(n/tot_group*100,2))

ggplot(data = region_ageg, aes(x = region, y = pct, fill = ageg)) +
  geom_col() +
  coord_flip()

이 표를 보기 좋게 노년층 비율이 많은 순으로 정렬해 보겠습니다.

 

list_order_old <- region_ageg %>%
  filter(ageg == "old") %>%
  arrange(pct)

list_order_old

노년층 비율을 내림차순 정렬해줍니다. 그다음 지역명 순서 변수를 만들어줍니다.

 

order <- list_order_old$region
order

 

이렇게 만든 변수를 활용해 그래프를 그려보겠습니다.

 

ggplot(data = region_ageg, aes(x = region, y = pct, fill = ageg)) +
  geom_col() +
  coord_flip() +
  scale_x_discrete(limits = order))

그런데 순서가 노년층이 중간에 있어서 보기에 좋지 않네요. 바꿔주도록 하겠습니다. 막대색을 순서대로 나열하려면 fill 파라미터에 지정할 변수의 범주(levels) 순서를 지정하면 됩니다.

 

levels(region_ageg$ageg)

 

현재 ageg변수는 levels가 없습니다. foctor()를 이용해 ageg변수를 foctor타입을 바꾼 후 , level파라미터를 이용해 순서를 정해줍니다.

 

region_ageg$ageg <- factor(region_ageg$ageg,
                           level = c("old", "middle", "young"))

levels(region_ageg$ageg)

 

결과를 보시면 원하는 순서로 바뀐 것을 알 수 있습니다. 이제 다시 그래프를 그려보죠.

 

ggplot(data = region_ageg, aes(x = region, y = pct, fill = ageg)) +
  geom_col() +
  coord_flip() +
  scale_x_discrete(limits = order)

이렇게 대구/경북지역이 노년층 비율이 가장 높다는 것을 알 수 있었습니다. 먼저 어떤 분석 프로젝트를 진행할 때에 이런식으로 기본적인 통계들을 내보고 그 결과들을 활용해 가설을 세우고, 변수들간의 상관관계를 분석하여 새로운 결론들을 도출해낼 수 있습니다. 그런 분석프로젝트의 아주 좋은 기본예제라고 생각합니다 :)