[재능상품] R을 통한 지니계수 계산 및 시각화

 정보

• 업무명     : R을 통한 지니계수 계산 및 시각화

• 작성자     : 박진만

• 작성일     : 2020-12-06

• 설   명      :

• 수정이력 :

 

 내용

[개요]

• 안녕하세요? 웹 개발 및 연구 개발을 담당하고 있는 해솔입니다.

• 다년간 축적된 경험 (기상학 학술 보고서 및 국/영문 학술 논문 게재, 블로그 운영, IT 회사 웹 개발 담당) 및 노하우를 바탕으로 개개인에게 맞춤형 솔루션을 수행할 수 있습니다.

• 특히 재능 플랫폼 (크몽, 오투잡, 해피캠퍼스, 레포트 월드)에서 누구보다도 경쟁력 있는 가격으로 양질의 서비스를 제공하고 있습니다.

  • 아스키 형식의 텍스트 (text) 파일부터 과학자료 형식 (HDF, H5, NetCDF, Grib, Grb) 및 Data Base (DB) 자료까지 다양한 형태의 자료를 이용하여 수집, 전처리, 분석, 시각화해 드립니다.

  • 또한 웹 사이트에 대한 정보를 이용한 웹 크롤링 및 그에 따른 엑셀 및 DB 구축도 가능합니다.

  • 아울러 기초 통계 (빈도분포, Prired t-test, Wilcoxn 등)에서 지도/비지도 학습을 통한 회귀모형 구축에 이르기 까지 효율적인 통계 정보를 제공합니다.

  • 최근 대한민국의 후속위성인 천리안위성 2A호 웹 서비스 서브시스템 및 환경위성 2B호 통합 자료처리 서브시스템에 대한 웹 개발을 수행하였습니다.

• 그리고 해솔 블로그에서는 다양한 기상학/천문학 정보와 더불어 사무 자동화/프로그래밍 언어를 소개하오니 방문 부탁드립니다.
  

• 좋은 하루 보내세요.

 

[재능플랫폼] 오투잡

[IT개발 - 응용프로그래밍] 통계 분석, 데이터 분석, 시각화를 성실하게 해 드립니다. - 재능마켓 �

 

[재능플랫폼] 크몽

데이터수집, 파싱, 크롤링 해 드립니다. | 50,000원부터 시작 가능한 총 평점 0점의 IT·프로그래밍,

 

 요청

[세부 사항]

• R 스크립트를 이용하여 지니계수 산출 함수 작성

• 아래와 같은 그래프 생성

 

 

 완료

[사용 OS]

• Windows 10

 

[사용 언어]

• R v4.0.2

 

[명세]

• 지니계수를 계산 및 시각화 함수 작성

giniValue <- function(data = data,income_col = income_col,plotFlag = FALSE,plotName = "result",plotPath = "./",...){
    
  addArg <- list(...)

  if(length(addArg) >= 1) {
    
   for (j in 1:length(addArg)) {
  
     nn <- as.data.frame(addArg[j])
     params <- data.frame(dum = c())

     if(nn[2,1] == "numeric"){ 
       
     for (k in 3:nrow(nn)) {
       param <- as.numeric(nn[k,1])
       params <- rbind(params,param)
     }
       
     } else if (nn[2,1] == "character") {
       
       for (k in 3:nrow(nn)) {
         param <- as.character(nn[k,1])
         params <- rbind(params,param)
       }
       
     } else {
       return("ERROR : 자료형의 입력이 잘못 되었습니다.")
     }
     
     # col name 재설정
     newColName = nn[1,1]
     colnames(params)[1] <- newColName
     # col name 재설정
     
     ### 데이터 필터링 ###
     data <- dplyr::inner_join(data,params,by=newColName)
     ### 데이터 필터링 ###
     
     # print(data)
     
     if(nrow(data) == 0) {
       return("조건에 맞는 데이터가 없습니다.")
     }

   }
   
  }

  # print(head(data))
  
  # 필터링된 데이터 선택 #
  x <- data[[income_col]]
  x_L1 <- sort(x)
  x_sum <- sum(x_L1)
  x_len <- length(x_L1)
  # 필터링된 데이터 선택 #
  
  # 지니계수 계산
  g = REAT::gini(x, lc = FALSE, lcx = "% of objects", 
                 lcy = "% of regarded variable", lctitle = "Lorenz curve", lcg = FALSE, 
                 lcgn = FALSE)
  
  
  # plotflag 여부에 따라 그림을 출력할지 말지 결정 #
  if(isTRUE(plotFlag)) {
    
  xx <- c()
  xy <- c()
  xw <- c()
  
  # x축 : 인구누적비율
  # y축 : 소득누적비율
  
  # 로렌츠 커브의 넓이를 계산하기 위한 코드
  tt <- 0
  
  for (i in 1:x_len) {
    
    isx <- i/x_len
    
    tt = tt+x_L1[i]
    
    isy <- tt/x_sum
    
    w <- (1/x_len) * isy
    
    xx <- append(xx,isx)
    xy <- append(xy,isy)
    xw <- append(xw,w)
    
  }
  # 로렌츠 커브의 넓이를 계산하기 위한 코드
  
  ############## w1 플로팅을 위한 위치좌표 설정  ############## 
  nis = seq(1,x_len,1)/x_len
  which_w1 = which.max(nis-xy)
  w1 = c(xx[which_w1],(nis[which_w1]+xy[which_w1])/2,0.5-round(sum(xw),4))
  ############## w1 플로팅을 위한 위치좌표 설정  ############## 
  
  # 결과를 데이터프레임에 저장
  res <- data.frame(x = xx, y = xy)
  
  ggplot(data = res) +
    theme_bw() +
    theme(text = element_text(size=14,face = "bold"),
          plot.title= element_text(size=18,face = "bold")) +
    geom_line(aes(x=xx,y=xy),size = 1.5,color = "blue") + 
    geom_text(aes(0.07, 0.97, label = paste0("G : ",round(g,4))),size = 5) +
    geom_text(aes(0.07, 0.92, label = paste0("N : ",length(x))),size = 5) +
    geom_text(aes(w1[1], w1[2], label = paste0("W1 = ",w1[3])),size = 5) +
    geom_abline(size = 1.5, color = "black") +
    # xlim(0,1) +
    # ylim(0,1) +
    xlab("인구 누적 비율") +
    ylab("소득 누적 비율") +
    labs(title = "gini value") +
    scale_x_continuous(expand = expansion(mult = c(0.01, 0.01)),breaks = c(0,0.25,0.5,0.75,1.0),limits = c(0,1)) +
    scale_y_continuous(expand = expansion(mult = c(0, 0)),breaks = c(0,0.25,0.5,0.75,1.0)) +
    ggsave(paste0(plotPath,plotName,".png"),dpi = 300,width = 8, height = 8)
  
  }
  
    # 지니계수 리턴
    return(g)

}

 

• 데이터 로드 및 함수 사용

######## 최초 실행시 한번만 수행 ###########
# install.packages("haven")
# install.packages("dplyr")
# install.packages("tidyr")
# install.packages("ggplot2")
# install.packages("labelled")
# install.packages("reldist")
# install.packages("lawstat")
# install.packages("REAT")
# install.packages("stringr")
# install.packages("openxlsx")
# install.packages("installr")
######## 최초 실행시 한번만 수행 ###########

######## 라이브러리 로드 ##############
library(installr)
installr::install.rtools()
Sys.setenv(R_ZIPCMD= "C:/Rtools/bin/zip")
# library(haven)
library(dplyr)
library(tidyr)
library(ggplot2)
library(foreign)
library(labelled)
library(reldist)
library(lawstat)
library(REAT)
library(stringr)
library(openxlsx)
######## 라이브러리 로드 ##############


# sav 데이터 읽기
#gini_sav_data <- read_spss("c:/workgini/gini_sav_data_Lnoin.sav") #경로는 임의로 변경 가능함
gini_sav_data <- read_spss("./gini_sav_data_Lnoin.sav") #경로는 임의로 변경 가능함


#컬럼이름 정하기
colnames(gini_sav_data) <- c("year","h_reg5","h_reg7","h_din","h_cin","h01_1","h01_1s","h01_4","h01_11","h01_110",
                             "ex_1","ex_2","ex_3","ex_4","ex_5","ex_6","ex_7","ex_8","ex_9",
                             "total_living_expenses","total_expen_notax","earned_income","business_income","proberty_income",
                             "etc_income","pub_pension","pri_pension","total_income","equalized_produc_income","equalized_income",
                             "rank10","old","gen")  #컬럼은 임의로 변경 가능함

# 추가 
# 해당 명령어를 이용하여 데이터의 형태를 간략하게 조회해볼 수 있음 
glimpse(gini_sav_data)


# factor to character
## dbl+l 데이터가 존재하는 경우 이를 character로 만들기 위한 작업임 ##
# dplyr::mutate(as.character(to_factor(컬림이름))...) 을 통해 변경 가능

gini_sav_data_Lnoin <- gini_sav_data %>%
  dplyr::mutate(h_reg5 = as.character(as_factor(h_reg5)),
                h_reg7 = as.character(as_factor(h_reg7)),
                h01_4 = as.character(as_factor(h01_4)),
                h01_110 = as.character(as_factor(h01_110)),
                h01_11 = as.character(as_factor(h01_11)),
                rank10 = as.numeric(rank10),
                gen = as.numeric(gen))   
                
                
                
g_full <- c()

for (y in 2007:2018) {
  
  g <- giniValue(data = gini_sav_data_Lnoin,income_col = "total_income", # 총 소득
                 param1 = c("year","numeric",y), # y를 대입 (2007 ~ 2018)
                 param2 = c("old","numeric",under65old), #65세 미만 (under65old)
                 plotFlag = TRUE,
                 plotPath = "C:\\tmpimg\\", # 이미지를 저장할 경로를 지정한다.
                 plotName = paste0("65세_미만_총_소득_",y,"_년도") ) # paste0 는 여러개의 문자열을 서로 붙여주는 역할
  
  print(g)
  
  g_full <- append(g_full,g) #여러개의 지니계수를 배열 형태로 저장할 수 있음.
  
}

print(g_full)
res = data.frame(year = 2007:2018,gini = g_full) # 결과
print(res) # 결과 출력

 

 

[소스 코드]

######## 최초 실행시 한번만 수행 ###########
# install.packages("haven")
# install.packages("dplyr")
# install.packages("tidyr")
# install.packages("ggplot2")
# install.packages("labelled")
# install.packages("reldist")
# install.packages("lawstat")
# install.packages("REAT")
# install.packages("stringr")
# install.packages("openxlsx")
# install.packages("installr")
######## 최초 실행시 한번만 수행 ###########

######## 라이브러리 로드 ##############
library(installr)
installr::install.rtools()
Sys.setenv(R_ZIPCMD= "C:/Rtools/bin/zip")
# library(haven)
library(dplyr)
library(tidyr)
library(ggplot2)
library(foreign)
library(labelled)
library(reldist)
library(lawstat)
library(REAT)
library(stringr)
library(openxlsx)
######## 라이브러리 로드 ##############


# sav 데이터 읽기
#gini_sav_data <- read_spss("c:/workgini/gini_sav_data_Lnoin.sav") #경로는 임의로 변경 가능함
gini_sav_data <- read_spss("./gini_sav_data_Lnoin.sav") #경로는 임의로 변경 가능함


#컬럼이름 정하기
colnames(gini_sav_data) <- c("year","h_reg5","h_reg7","h_din","h_cin","h01_1","h01_1s","h01_4","h01_11","h01_110",
                             "ex_1","ex_2","ex_3","ex_4","ex_5","ex_6","ex_7","ex_8","ex_9",
                             "total_living_expenses","total_expen_notax","earned_income","business_income","proberty_income",
                             "etc_income","pub_pension","pri_pension","total_income","equalized_produc_income","equalized_income",
                             "rank10","old","gen")  #컬럼은 임의로 변경 가능함

# 추가 
# 해당 명령어를 이용하여 데이터의 형태를 간략하게 조회해볼 수 있음 
glimpse(gini_sav_data)


# factor to character
## dbl+l 데이터가 존재하는 경우 이를 character로 만들기 위한 작업임 ##
# dplyr::mutate(as.character(to_factor(컬림이름))...) 을 통해 변경 가능


## 수정사항 ##
# to_factor 대신 as_factor 을 사용하여도 동일한 효과를 쓸 수 있음
gini_sav_data_Lnoin <- gini_sav_data %>%
  dplyr::mutate(h_reg5 = as.character(as_factor(h_reg5)),
                h_reg7 = as.character(as_factor(h_reg7)),
                h01_4 = as.character(as_factor(h01_4)),
                h01_110 = as.character(as_factor(h01_110)),
                h01_11 = as.character(as_factor(h01_11)),
                rank10 = as.numeric(rank10),
                gen = as.numeric(gen))     

##### 컬럼 종류 #####

# year : 조사년도
# h_reg5 : 5개 권역별 지역구분
# h_reg7 : 7개 전국지역 구분
# h01_4 : 성별구분(1 남, 2 여) 
# h01_110 : 가구형태
# total_living_expenses : 총 생활비
# old : 세대주 나이
# earned_income : 근로소득
# business_income : 사업소득
# total_income : 전체소득
# equalized_income : 균등화소득
# rank10 : 10분위 구분 
# gen : 일반, 노인세대 구분(2는 노인)

##### 컬럼 종류 #####

############## sav ##################



############ 데이터 읽기 및 가공종료 ##############


########################################### FUNCTION ###########################################
giniValue <- function(data = data,income_col = income_col,plotFlag = FALSE,plotName = "result",plotPath = "./",...){
    
  addArg <- list(...)

  if(length(addArg) >= 1) {
    
   for (j in 1:length(addArg)) {
  
     nn <- as.data.frame(addArg[j])
     params <- data.frame(dum = c())

     if(nn[2,1] == "numeric"){ 
       
     for (k in 3:nrow(nn)) {
       param <- as.numeric(nn[k,1])
       params <- rbind(params,param)
     }
       
     } else if (nn[2,1] == "character") {
       
       for (k in 3:nrow(nn)) {
         param <- as.character(nn[k,1])
         params <- rbind(params,param)
       }
       
     } else {
       return("ERROR : 자료형의 입력이 잘못 되었습니다.")
     }
     
     # col name 재설정
     newColName = nn[1,1]
     colnames(params)[1] <- newColName
     # col name 재설정
     
     ### 데이터 필터링 ###
     data <- dplyr::inner_join(data,params,by=newColName)
     ### 데이터 필터링 ###
     
     # print(data)
     
     if(nrow(data) == 0) {
       return("조건에 맞는 데이터가 없습니다.")
     }

   }
   
  }

  # print(head(data))
  
  # 필터링된 데이터 선택 #
  x <- data[[income_col]]
  x_L1 <- sort(x)
  x_sum <- sum(x_L1)
  x_len <- length(x_L1)
  # 필터링된 데이터 선택 #
  
  # 지니계수 계산
  g = REAT::gini(x, lc = FALSE, lcx = "% of objects", 
                 lcy = "% of regarded variable", lctitle = "Lorenz curve", lcg = FALSE, 
                 lcgn = FALSE)
  
  
  # plotflag 여부에 따라 그림을 출력할지 말지 결정 #
  if(isTRUE(plotFlag)) {
    
  xx <- c()
  xy <- c()
  xw <- c()
  
  # x축 : 인구누적비율
  # y축 : 소득누적비율
  
  # 로렌츠 커브의 넓이를 계산하기 위한 코드
  tt <- 0
  
  for (i in 1:x_len) {
    
    isx <- i/x_len
    
    tt = tt+x_L1[i]
    
    isy <- tt/x_sum
    
    w <- (1/x_len) * isy
    
    xx <- append(xx,isx)
    xy <- append(xy,isy)
    xw <- append(xw,w)
    
  }
  # 로렌츠 커브의 넓이를 계산하기 위한 코드
  
  ############## w1 플로팅을 위한 위치좌표 설정  ############## 
  nis = seq(1,x_len,1)/x_len
  which_w1 = which.max(nis-xy)
  w1 = c(xx[which_w1],(nis[which_w1]+xy[which_w1])/2,0.5-round(sum(xw),4))
  ############## w1 플로팅을 위한 위치좌표 설정  ############## 
  
  # 결과를 데이터프레임에 저장
  res <- data.frame(x = xx, y = xy)
  
  ggplot(data = res) +
    theme_bw() +
    theme(text = element_text(size=14,face = "bold"),
          plot.title= element_text(size=18,face = "bold")) +
    geom_line(aes(x=xx,y=xy),size = 1.5,color = "blue") + 
    geom_text(aes(0.07, 0.97, label = paste0("G : ",round(g,4))),size = 5) +
    geom_text(aes(0.07, 0.92, label = paste0("N : ",length(x))),size = 5) +
    geom_text(aes(w1[1], w1[2], label = paste0("W1 = ",w1[3])),size = 5) +
    geom_abline(size = 1.5, color = "black") +
    # xlim(0,1) +
    # ylim(0,1) +
    xlab("인구 누적 비율") +
    ylab("소득 누적 비율") +
    labs(title = "gini value") +
    scale_x_continuous(expand = expansion(mult = c(0.01, 0.01)),breaks = c(0,0.25,0.5,0.75,1.0),limits = c(0,1)) +
    scale_y_continuous(expand = expansion(mult = c(0, 0)),breaks = c(0,0.25,0.5,0.75,1.0)) +
    ggsave(paste0(plotPath,plotName,".png"),dpi = 300,width = 8, height = 8)
  
  }
  
    # 지니계수 리턴
    return(g)

}



###### 실제로 적용하는 예시 1 #########

## 조건 1 
# 65세 이상 총 소득
# 2007년부터 2018년까지

# 비어있는 배열 생성
g_full <- c()

for (y in 2007:2018) {
  
  g <- giniValue(data = gini_sav_data_Lnoin,income_col = "total_income", # 총 소득
                 #param1 = c("year","numeric",y), # y를 대입 (2007 ~ 2018)
                 #param2 = c("old","numeric",over65old),
                 param3 = c("old","numeric",40), #65세 이상 (over65old)
                 plotFlag = TRUE,
                 plotPath = "C:\\tmpimg\\", # 이미지를 저장할 경로를 지정한다.
                 plotName = paste0("65세_이상_총_소득_",y,"_년도") ) # paste0 는 여러개의 문자열을 서로 붙여주는 역할
  
  print(g)
  
  g_full <- append(g_full,g) #여러개의 지니계수를 배열 형태로 저장할 수 있음.
  
  
}

print(g_full)
res = data.frame(year = 2007:2018,gini = g_full) # 결과
print(res) # 결과 출력

###### 실제로 적용하는 예시 1 #########


###### 실제로 적용하는 예시 2 #########

## 조건 2
# 65세 미만 총 소득
# 2007년부터 2018년까지

# 비어있는 배열 생성
g_full <- c()

for (y in 2007:2018) {
  
  g <- giniValue(data = gini_sav_data_Lnoin,income_col = "total_income", # 총 소득
                 param1 = c("year","numeric",y), # y를 대입 (2007 ~ 2018)
                 param2 = c("old","numeric",under65old), #65세 미만 (under65old)
                 plotFlag = TRUE,
                 plotPath = "C:\\tmpimg\\", # 이미지를 저장할 경로를 지정한다.
                 plotName = paste0("65세_미만_총_소득_",y,"_년도") ) # paste0 는 여러개의 문자열을 서로 붙여주는 역할
  
  print(g)
  
  g_full <- append(g_full,g) #여러개의 지니계수를 배열 형태로 저장할 수 있음.
  
}

print(g_full)
res = data.frame(year = 2007:2018,gini = g_full) # 결과
print(res) # 결과 출력


###### 실제로 적용하는 예시 2 #########

 

[결과물]

• 지니계수 결과 이미지

 

 참고 문헌

[논문]

• 없음

[보고서]

• 없음

[URL]

• 없음

 

 문의사항

[기상학/프로그래밍 언어]

• sangho.lee.1990@gmail.com

[해양학/천문학/빅데이터]

• saimang0804@gmail.com