분석 주제 (연구문제 설정)

1. 서울시 골목상권 별 카페 총 매출액에 관련된 주제를 선정한 이유 & 목표

많은 자영업 분야들 중에서 자영업자들이 창업하는 비율을 봤을 때, 특히 카페 업종을 창업하는 비율이 급증하였다는 이야기를 많이 접하였습니다. 특히나 젊은 층은 공부를 하기 위해서나, 지인들 과의 만남 등을 위해 더더욱 카페를 많이 방문하기도 합니다. 그래서 저희 조는 서울시 골목상권 별 카페의 총 매출액의 주제가 요즘 시기와 잘 맞는 주제이기 때문에 조금 더 재밌게 분석할 수 있을 것 같다고 생각했습니다. 저희 15조가 이번 분석 주제에서 상권 데이터를 사용한 이유는, 카페를 창업할 때 상권이 매우 중요한 요소이기 때문이었습니다. 그래서 상권에 관한 데이터를 수집하여 분석 데이터로 사용하였습니다. 그래서 수집한 서울시 골목상권 별 상주인구, 직장인 수, 집객 시설 등의 상권 데이터를 바탕으로 어떤 변수들이 총 매출에 영향을 주는 지 알아보고자 하였습니다. 이러한 이유로 저희 팀은 서울시 우리마을가게 상권 데이터를 이용하여 분석을 진행하기로 하였습니다.

2. 데이터

library(tidyverse)
library(visdat)
library(caTools)
library(readr)
library(rlist)
library(rgdal)
library(sp)
library(plotly)
library(ggplot2)
library(ggthemes)
library(RColorBrewer)
library(raster)
library(cowplot)
library(reshape)
library(corrplot)
library(car)
library(Metrics)
library(patchwork)
library(gridExtra)
library(reshape2)
library(rcompanion)
library(sjPlot)
library(GGally)
library(factoextra)

# data 폴더안에 있는 모든 데이터 리스트
files <- list.files(path="data/main/")
files
## [1] "earningsSpendings.csv"           "estimatedFloatingPopulation.csv"
## [3] "estimatedSales.csv"              "infrastructure.csv"             
## [5] "numOfEmployee.csv"               "settledPopulation.csv"
# 데이터를 한번에 다 불러오기
csvFiles <- list()
for(i in 1:length(files)){
  code <- paste0("read_csv(paste0('data/main/', '", files[i],"'), locale = locale('ko', encoding = 'utf-8'))") # 함수식 설정
  result <- list(eval(parse(text=code))) # 함수식 실행
  result <- setNames(result, gsub("*.csv","",files[i])) # 이름을 파일명과 같이 설정
  csvFiles <- append(csvFiles, result) # 데이터를 리스트안에 몰아 넣기

  rm(result, code)
}

분석을 위해 서울 열린데이터 광장에서 제공하는 “서울시 우리마을가게 상권분석서비스”의 상권배후지 매출액 데이터와 다른 5개 데이터를 사용하였다.

데이터 파일명
추정매출 estimatedSales.csv
추정유동인구 estimatedFloatingPopulation.csv
상주인구 settledPopulation.csv
직장인구 numOfEmployee.csv
소득소비 earningsSpendings.csv
집객시설 infrastructure.csv

출처 : https://data.seoul.go.kr/

2.1 데이터 전처리

소득소비 데이터의 소득_구간_코드가 문자열로 되어있어서 연단위로 합칠 수 없으므로 제거해 주었다.

명목변수인 상권_변화_지표도 제거해 주고 카페에 해당되는 업종 코드명 “커피-음료”만 추출하였다. 다른 7개의 데이터에는 없는데 추정매출에만 있는 서비스_업종_코드, 서비스_업종_코드_명도 삭제해주었다.

# earningsSpendings에서 소득_구간_코드 열 제거
csvFiles$earningsSpendings <- csvFiles$earningsSpendings %>% within(rm(소득_구간_코드))

csvFiles$estimatedSales$서비스_업종_코드_명 %>% unique
##  [1] "전자상거래업"           "조명용품"               "인테리어"              
##  [4] "철물점"                 "가전제품"               "가구"                  
##  [7] "애완동물"               "화초"                   "섬유제품"              
## [10] "운동/경기용품"          "화장품"                 "문구"                  
## [13] "서적"                   "의료기기"               "의약품"                
## [16] "시계및귀금속"           "안경"                   "신발"                  
## [19] "일반의류"               "반찬가게"               "청과상"                
## [22] "수산물판매"             "육류판매"               "미곡판매"              
## [25] "핸드폰"                 "컴퓨터및주변장치판매"   "편의점"                
## [28] "슈퍼마켓"               "노래방"                 "고시원"                
## [31] "세탁소"                 "피부관리실"             "네일숍"                
## [34] "미용실"                 "자동차미용"             "자동차수리"            
## [37] "스포츠클럽"             "PC방"                   "당구장"                
## [40] "한의원"                 "치과의원"               "일반의원"              
## [43] "스포츠 강습"            "예술학원"               "외국어학원"            
## [46] "일반교습학원"           "커피-음료"              "호프-간이주점"         
## [49] "분식전문점"             "치킨전문점"             "패스트푸드점"          
## [52] "제과점"                 "양식음식점"             "일식음식점"            
## [55] "중식음식점"             "한식음식점"             "완구"                  
## [58] "여관"                   "자전거 및 기타운송장비" "골프연습장"            
## [61] "가방"                   "부동산중개업"           "가전제품수리"
csvFiles$estimatedSales <- csvFiles$estimatedSales %>% filter(`서비스_업종_코드_명` == "커피-음료")
csvFiles$estimatedSales <- csvFiles$estimatedSales %>% within(rm("서비스_업종_코드", "서비스_업종_코드_명"))

데이터가 2020년 4분기까지 밖에 없어서 비교를 위해 동일한 기간을 설정해주었다.

# 데이터 기간 설정
year <- 2020
exclude_quarter <- 1

# 데이터들의 1열~6열까지 변수명을 통일
for(i in 1:length(csvFiles)){
  colnames(csvFiles[[i]])[1:6] <- c("기준_년_코드", "기준_분기_코드", "상권_구분_코드", "상권_구분_코드_명", "상권_코드", "상권_코드_명")
}

# 분기별 데이터를 연간 데이터로 합쳐주는 함수
aggregateYearly <- function(data){
  data <- data %>% within(rm("기준_분기_코드", "상권_구분_코드", "상권_구분_코드_명", "상권_코드_명"))
  # 분기별로 나눠져있는 변수들을 연단위로 합치기
  data <- data %>% group_by(`상권_코드`, `기준_년_코드`) %>% dplyr::summarise(across(everything(), sum) %>% round(0), .groups="keep")
  return(data)
}

# 2,3,4분기만 남기는 함수
cutTime <- function(data){
  data <- data %>% filter(`기준_분기_코드` != exclude_quarter)
  data <- data %>% filter(`기준_년_코드` == year)
  data2020 <- aggregateYearly(data)
  return(list(data2020))
}

dataList2020 <- list()

for(i in 1:length(csvFiles)){
  cutData <- cutTime(csvFiles[[i]])
  result <- list(cutData[[1]])
  result <- setNames(result, names(csvFiles[i]))
  dataList2020 <- append(dataList2020, result)
}
rm(result) # 더이상 필요없는 변수 삭제

데이터가 제대로 추출 되었는지 확인해주었다.

checkLength <- function(){
  recordLength20 <- list()

  for(i in 1:length(csvFiles)){
    result <- nrow(dataList2020[[i]])
    result <- setNames(result, names(dataList2020[i]))
    recordLength20 <- append(recordLength20, result)
  }
  as.data.frame(recordLength20) %>% t
}

checkLength()
##                             [,1]
## earningsSpendings           1010
## estimatedFloatingPopulation 1010
## estimatedSales              1009
## infrastructure              1010
## numOfEmployee               1010
## settledPopulation           1010

추정매출 데이터에 관측값이 하나 모자르기 때문에 확인해보았다.

diff <- setdiff(dataList2020$earningsSpendings$상권_코드, dataList2020$estimatedSales$상권_코드)
diff
## [1] 1000750
csvFiles$estimatedSales[csvFiles$estimatedSales$상권_코드 == diff,] %>%
  group_by(`상권_코드_명`) %>%
  summarise(mean(`분기당_매출_금액`))

상권코드가 1000750인 사당로23길 상권에 매출 데이터가 없는것 같다. 따라서 다른 데이터들에도 1000750 상권의 데이터를 삭제해 주었다.

for(i in 1:length(dataList2020)){
  dataList2020[[i]] <- subset(dataList2020[[i]], `상권_코드` != diff)
}
checkLength()
##                             [,1]
## earningsSpendings           1009
## estimatedFloatingPopulation 1009
## estimatedSales              1009
## infrastructure              1009
## numOfEmployee               1009
## settledPopulation           1009

이제 모든 데이터의 관측값이 1009개로 동일해졌다.

8가지 데이터들 중에서 사용할 변수만 추출하여 주었다.

# 유동인구
EstimatedFloatingPopulation2020 <- dataList2020$estimatedFloatingPopulation[1:24]

# 직장인수
EmployeeNum2020 <- dataList2020$numOfEmployee[1:11]

# 상주인구
SettledPopulation2020 <- dataList2020$settledPopulation[c(1:11, 24)]

# 소득소비
earningsSpendings2020 <- dataList2020$earningsSpendings

# 집객시설
infrastructure2020 <- dataList2020$infrastructure[c(1:2, 4, 9:13, 16, 21:22)]

# 추정매출
sales2020 <- dataList2020$estimatedSales[c(1:3, 74)]

# 분기당 매출금액 변수명을 '총매출금액'으로 통일
colnames(sales2020)[3] <- "총매출금액"

분석을 위해 하나의 테이블로 만들어주었다.

data2020 <- list(sales2020,
                 EstimatedFloatingPopulation2020, 
                 EmployeeNum2020, 
                 SettledPopulation2020, 
                 earningsSpendings2020, 
                 infrastructure2020) %>% reduce(left_join, by = c("상권_코드", "기준_년_코드"))

NA값을 확인하였고 총 3625개의 NA값이 있다. 데이터를 확인해보았을 때 집객시설 데이터에 NA값이 매우 많았고 이는 해당 집객시설이 없다는 것으로 판단하고 모두 0으로 채워 넣어주었다.

vis_dat(data2020)

sum(is.na(data2020))
## [1] 3625
data2020[is.na(data2020)] <- 0

3. EDA

3.1 종속변수 매출총금액 EDA

우선 종속변수인 매출 총금액을 시각화해보았다.

3.1.1 총 매출 금액 히스토그램

다음으로 상권별 총 매출금액에 따른 상권의 빈도 분포를 히스토그램으로 그려보았다.

ggplot(sales2020, aes(x=총매출금액)) +
  geom_histogram(fill='skyblue', colour='black')+
  ggtitle('2020년도 총 매출 금액 히스토그램')

총 매출 금액이 왼쪽에 몰려있는 positive skewed data임을 알 수 있다. 해당 히스토그램에 따르면 종속변수인 총 매출금액은 정규분포를 따르지 않을 것으로 보인다.

3.1.2 분기별 매출금액 변화

분기별 총 매출금액의 변화를 라인그래프를 통해 알아보았다.

temp <- csvFiles$estimatedSales %>%
                group_by(기준_분기_코드) %>%
                dplyr::summarise(money= mean(분기당_매출_금액))


ggplot(temp, aes(x=기준_분기_코드, y=money))+
  geom_line()

코로나가 발생한 2분기를 기점으로 총 매출액이 2분기에서 3분기까지 줄어들다가 3분기 이후 급격하게 감소하고 있다. 이는 카페 매출이 겨울에 더 줄어드는 경향이 원인인지 사회적 거리두기 때문인지는 추가적인 분석이 필요해 보인다.

3.1.3 변수별 매출금액 평균 - 주중, 주말

다음으로 주중 주말과 요일, 시간대, 성별, 연령대 별 변수에 따른 매출액 양상을 막대그래프를 통해 알아보았다. 매출금액은 각 변수에 다른 상권들의 매출액 평균으로 집계하였다.

다음은 주중 주말의 매출금액에 대한 그래프이다.

draw.money <- function(df, var.names){
  df <- df %>% dplyr::select(var.names)
  df.t <- df %>% t
  plot_ly(x = rownames(df.t), y = df.t[, 1], type = "bar")
}

# 평균으로 매출금액 집계
avgSales2020 <- dataList2020$estimatedSales[28:50] %>% colMeans() %>% t %>% as.data.frame
var.names <- c("주중_매출_금액","주말_매출_금액")
draw.money(avgSales2020, var.names)

주말보다는 주중에 매출이 2배 넘게 더 많다. 변수에서 주중과 주말을 나누어 볼 필요가 있어보인다.

3.1.4 2020년 변수 별 매출금액 평균 - 요일

각 요일 별 매출금액에 대한 그래프이다.

draw.money(avgSales2020, c("월요일_매출_금액","화요일_매출_금액", "수요일_매출_금액", "목요일_매출_금액", "금요일_매출_금액", "토요일_매출_금액", "일요일_매출_금액"))

요일별 매출금액의 경우 일요일의 매출이 가장 낮다.

3.1.5 2020년 변수 별 매출금액 평균 - 시간대

시간대 별 매출금액에 대한 그래프이다.

draw.money(avgSales2020, c("시간대_00~06_매출_금액","시간대_06~11_매출_금액","시간대_11~14_매출_금액","시간대_14~17_매출_금액","시간대_17~21_매출_금액","시간대_21~24_매출_금액"))

시간대별 매출금액의 경우 점심시간 전후인 11시에서 14시 사이의 매출금액이 가장 높다

3.1.6 2020년 변수 별 매출금액 평균 - 성별

draw.money(avgSales2020, c("남성_매출_금액","여성_매출_금액" ))

성별 매출금액을 보았을 때, 여성이 남성보다 카페에 더 많은 돈을 사용한다는 것을 알 수 있다. 성별에 따라 변수를 나누는 것이 좋을 것으로 보인다.

3.1.7 2020년 변수 별 매출금액 평균 - 연령대

draw.money(avgSales2020, c("연령대_10_매출_금액",
                           "연령대_20_매출_금액",
                           "연령대_30_매출_금액",
                           "연령대_40_매출_금액",
                           "연령대_50_매출_금액",
                           "연령대_60_이상_매출_금액"))

연령대 별 매출금액의 경우, 10대의 카페 매출액이 가장 낮고 20대가 가장 매출액이 높다. 그리고 60대로 갈 수록 매출액이 점차 낮아지는 양상을 보인다. 가장 낮은 매출액을 보이는 10대와 가장 높은 매출을 보이는 20대, 그리고 가장 낮은 매출을 보이는 60대를 각각 나누어 봐야 할 필요성이 있다.

매출금액을 대상으로 EDA를 그려봄으로써 수많은 독립변수들을 어떻게 추려야 할지 간략하게 알아보았다.

3.2 점포수 대비 매출액이 높은 상권과 낮은 상권 비교

다음으로 점포수 대비 매출액이 높은 상권의 변수들과 낮은 상권의 변수들을 비교해봄으로써 독립변수에 대한 EDA를 더욱 세부적으로 해보고자 한다.

점포수 대비가 아닌 단순 매출금액을 대상으로 상권들을 지도로 시각화해보면 다음과 같다.

지도 데이터 출처 : http://www.gisdeveloper.co.kr/?p=2332_

골목상권 좌표 출처 : https://data.seoul.go.kr/dataList/OA-15560/S/1/datasetView.do

# 상권코드 추출
districtCode <- csvFiles$numOfEmployee[5:6] %>% unique
districtCode$상권_코드 <- as.factor(districtCode$상권_코드)

# 골목상권 좌표 추출
location <- read_csv("data/eda/seoulGolmok.csv", locale = locale('ko', encoding = 'utf-8'))
oldCoord <- location[6:7]
colnames(oldCoord) <- c("long","lat")
coordinates(oldCoord) <- c("long","lat")

proj4string(oldCoord) <- CRS("+init=epsg:5181")
newCoord <- spTransform(oldCoord, CRS("+init=epsg:4326"))
newCoord <- newCoord %>% as.data.frame
location <- bind_cols(location, newCoord)
location <- location[c(4, 11:12)]
location$상권_코드 <- as.factor(location$상권_코드)
location <- location[!(location$상권_코드==1000750),]
mapShp <- shapefile("SIG/TL_SCCO_SIG.shp")
map <- fortify(mapShp, region="SIG_CD")
map$id <- as.numeric(map$id)
seoulMap <- map[map$id<=11740,]

# 좌표계 변경
convert <- seoulMap[1:2]
coordinates(convert) <- c("long","lat")
proj4string(convert) <- CRS("+init=epsg:5179")
converted <- spTransform(convert, CRS("+init=epsg:4326"))
converted <- converted %>% as.data.frame
seoulMap[1:2] <- converted

# 서울시 구 좌표 추출
sigCode <- read_csv("data/eda/sig_code.csv", locale = locale('ko', encoding = 'utf-8'))
sigCode <- sigCode[c(3, 6:7)]
colnames(sigCode) <- c("gu","lat","long")

g <- seoulMap %>% ggplot() + 
  geom_polygon(aes(x=long, y=lat, group=group), fill="grey", color="white") + 
  geom_text(data=sigCode, aes(x=long, y=lat, label=gu), size=2.5)
sales2020$상권_코드 <- as.factor(sales2020$상권_코드)
mapping <- left_join(location, sales2020[,c(1,3)], by="상권_코드")
mapping <- mapping %>% arrange(`총매출금액`)

g + geom_point(data=mapping,
               aes(x=long, y=lat, color=`총매출금액`, size=`총매출금액`, alpha=`총매출금액`)) +
  guides(alpha = FALSE, size = FALSE) +
  scale_alpha(range=c(0.3,0.5)) +
  scale_size_continuous(range=c(1,15)) +
  scale_colour_gradient(low="#ffcccc", high="#ff0000") +
  theme_void()

지도를 보면 명동에 위치한 상권의 매출금액이 가장 높은 것을 볼 수 있다.

하지만 명동에 다른 상권보다 점포수가 많다면 매출금액이 높은 것은 당연한 결과일 수 있다. 그렇기 때문에 점포수를 고려한 매출금액을 통해 고매출 상권과 저매출 상권이 어느 곳일지 더 정확히 알아볼 필요가 있다.

3.2.1 점포수와 매출금액의 산점도

점포수와 매출금액에 선형관계가 있는지 보기 위한 산점도 그래프를 그려보았다.

ggplot(data2020, aes(x=점포수, y=총매출금액)) + geom_point(alpha= 0.3)

상권코드 별로 점포수가 많아질수록 총매출금액이 증가하는 것을 볼 수 있다.

따라서 점포수별 매출금액의 분포를 확인하기 위해 점포수 대비 매출금액을 계산하여 해당 그래프의 박스 플롯을 그려보았다.

3.2.2 점포수 대비 매출금액

sales.ratio <- data2020
# 점포수 대비 총매출금액으로 총매출금액 바꿔주기
sales.ratio$점포수_대비_매출금액<- sales.ratio$총매출금액/sales.ratio$점포수

# 지도
sales.ratio$상권_코드 <- as.factor(sales.ratio$상권_코드)
mapping <- left_join(location, sales.ratio[,c(1, 2, 66)], by="상권_코드")
mapping <- mapping %>% arrange(`점포수_대비_매출금액`)

g + geom_point(data=mapping,
               aes(x=long, y=lat, color=`점포수_대비_매출금액`, size=`점포수_대비_매출금액`, alpha=`점포수_대비_매출금액`)) +
  guides(alpha = FALSE, size = FALSE) +
  scale_alpha(range=c(0.3,0.5)) +
  scale_size_continuous(range=c(1,15)) +
  scale_colour_gradient(low="#ffcccc", high="#ff0000") +
  theme_void()

점포수 대비 매출금액은 영등포구가 높은것을 알 수 있다.

점포수 대비 매출금액의 분포를 확인해 보자.

# 점포수 대비 총 매출금액에 대한 박스플롯 그리기
plot_ly(sales.ratio, x=~점포수_대비_매출금액, type='box')
# 1분위, 3사분위 수 출력
(quantile(sales.ratio$점포수_대비_매출금액, prob=c(0.25, 0.75)))
##      25%      75% 
## 12045866 32406901

결과를 보면 12,045,866원이 제1 사분위수 값이고 32,406,901원이 제3 사분위수 값임을 알 수 있다. 해당 결과를 바탕으로 저매출 상권은 점포수 대비 매출금액이 12,045,866원 이하인 상권으로, 고매출 상권은 점포수 대비 매출금액이 32,406,901원 이상인 상권으로 설정하였다.

sales.ratio <- sales.ratio %>%
  mutate(상권매출구분 = ifelse(
    점포수_대비_매출금액 <= 12045866,
    '저매출상권',
    ifelse(점포수_대비_매출금액 >= 32406901, '고매출상권',' ')))

cpr <- sales.ratio %>%
  group_by(상권매출구분) %>%
  summarize_all(mean)

이제 두 상권에 대한 변수들의 비교를 간단히 해보고자 한다.

3.2.3 두 상권 별 유동인구 변수 비교

성별, 시간대별, 연령대별, 요일별 그래프를 그려보았지만 네개의 그래프 중 저매출상권과 고매출 상권사이의 차이가 보이는 그래프만 나타내보았다.

cpr <- cpr[c(which(cpr$상권매출구분== '저매출상권'), which(cpr$상권매출구분== '고매출상권')),]

# 데이터 조작과 그래프 그리는 함수선언
make_df <-function(df, varNames){
  temp <- df %>% dplyr::select(varNames)
  temp <-melt(as.data.frame(temp), id.vars= c("상권매출구분"))
  return(temp)
}

draw_df1 <- function(df, varNames){
  ggplot(df, aes(x= variable, y= value)) +
    geom_bar(stat="identity") + facet_wrap(~상권매출구분, nrow=1) +
    theme(axis.text.x=element_text(angle=25, hjust=1, size=7)) +
    scale_fill_brewer(palette="Spectral")
}

draw_df2 <- function(df){
  ggplot(df, aes(x=factor(상권매출구분), y=value, fill=variable)) +
    geom_bar(stat="identity", position="fill", width=0.4, colour="black") +
    scale_y_continuous(labels = scales::percent_format()) +
    theme(legend.position="right") + 
    theme(legend.title = element_text(size = 10)) + # legend title
    theme(legend.text = element_text(size = 8)) +
    scale_fill_brewer(palette="Spectral") +
    theme(axis.text.x=element_text(angle=25, hjust=1, size=7))
}

# 연령대별 유동인구
graph_df <- make_df(cpr, c("상권매출구분",
                           "연령대_10_유동인구_수",
                           "연령대_20_유동인구_수",
                           "연령대_30_유동인구_수",
                           "연령대_40_유동인구_수",
                           "연령대_50_유동인구_수",
                           "연령대_60_이상_유동인구_수"))
draw_df2(graph_df)+ ggtitle("연령대 별 유동인구")

연령대 별 유동인구에서 고매출상권엔 연령대 20,30대 유동인구가 저매출상권에 비해 많은 반면, 저매출 상권의 경우엔 연령대 60 이상의 유동인구가 많다.

3.2.4 두 상권 별 직장인구 변수 비교

# 직장인구수 
## 성별 직장인구
graph_df <- make_df(cpr, c("여성_직장_인구_수",
                          "남성_직장_인구_수",
                          "상권매출구분"))
emp1 <- draw_df2(graph_df) + ggtitle("성별 직장인구")
## 연령대별 유동인구
graph_df <- make_df(cpr, c("상권매출구분",
                          "연령대_10_직장_인구_수",
                          "연령대_20_직장_인구_수",
                          "연령대_30_직장_인구_수",
                          "연령대_40_직장_인구_수", 
                          "연령대_50_직장_인구_수",
                          "연령대_60_이상_직장_인구_수"))
emp2 <- draw_df2(graph_df) + ggtitle("연령대 별 직장인구")
grid.arrange(emp1, emp2, nrow=1, ncol=2)

성별 직장인구에선 매출이 높은 고매출상권에서 남성 직장 인구수가 좀 더 높은 것으로 보인다.

연령대 별 직장인구에선 고매출 상권엔 비교적으로 20, 30대 직장인구가 많고 저매출 상권엔 50, 60대 직장인구 수가 많은 것을 볼 수 있다.

3.2.5 두 상권 별 상주인구 변수 비교

연령대, 성별 직장인구 중 차이가 조금 보이는 성별 연령대 별 직장인구 만을 나타내었다.

## 연령대별 유동인구
graph_df <- make_df(cpr, c("상권매출구분",
                          "연령대_10_상주인구_수",
                          "연령대_20_상주인구_수",
                          "연령대_30_상주인구_수",
                          "연령대_40_상주인구_수", 
                          "연령대_50_상주인구_수",
                          "연령대_60_이상_상주인구_수"))
draw_df2(graph_df) + ggtitle("연령대 별 상주인구_수")

저매출 상권의 경우 연령대 60이상인 상주인구수가 가장 많은 비율을 차지하고 있다. 그리고 20,30 대 상주인구의 비율이 고매출상권보단 적은 비율을 기록한다. 이로 미뤄보아 저매출상권은 젊은 층이 많이 상주하지 않는 지역임을 알 수 있다.

직장인구, 상주인구, 유동인구 변수들을 보았을 때 연령대별 인구가 매출액에 영향을 끼칠 수 있는 요소가 될 수 있을 것으로 추정된다.

3.2.6 두 상권 별 소득금액, 지출금액 비교

# 소득금액, 지출금액
graph_df <- make_df(cpr, c("월_평균_소득_금액",
                          "상권매출구분"))
income <- draw_df1(graph_df)
graph_df <- make_df(cpr, c("문화_지출_총금액",
                          "교육_지출_총금액",
                          "상권매출구분"))
expense <- draw_df1(graph_df)
grid.arrange(income, expense, nrow=1, ncol=2)

확실히 고매출 상궈의 월평균 소득금액이 높다. 그러나 문화지출, 교육 지출 총금액의 경우 저매출상권의 금액이 좀 더 높은 것을 볼 수 있다.

3.2.7 두 상권 별 집객시설 비교

# 집객시설
graph_df<- make_df(cpr, c("상권매출구분",
                          "관공서_수",
                          "유치원_수",
                          "초등학교_수",
                          "중학교_수",
                          "고등학교_수",
                          "극장_수",
                          "지하철_역_수",
                          "버스_정거장_수"))
ggplot(graph_df, aes(x="", y=value, fill=variable)) +
  geom_bar(stat="identity", position="fill", width=0.4, colour="black") +
  coord_polar("y") +
  scale_y_continuous(labels = scales::percent_format()) +
  theme(legend.position="right") +
  theme(legend.title = element_text(size = 10)) + # legend title
  theme(legend.text = element_text(size = 8)) +
  scale_fill_brewer(palette="Spectral") +
  facet_wrap(.~상권매출구분)

총 집객시설의 경우 고매출 상권엔 저매출상권에 비해 관공서와 극장 수의 비율이 높은 것을 볼 수 있다. 저매출 상권엔 상대적으로 유치원, 초등학교, 중학교와 같은 교육시설의 비율이 높은 것을 볼 수 있다.

3.3 총직장인구수와 총매출금액의 산점도

마지막으로 직장인구, 상주인구, 유동인구와 총매출금액에 대한 산점도를 작성해보았다. 그 중 약간의 상관관계를 보였던 직장인구수와 매출금액의 산점도만을 나타내보았다.

# 2020년 총 직장인구수와 총 매출금액에 대한 상관계수 구하기
cor(data2020$총_직장_인구_수, data2020$총매출금액)
## [1] 0.4829954
# (편의상)영어로 변수 이름 바꿔주기
total.work.population2020 <- data2020$총_직장_인구_수
total.sales2020 <- data2020$총매출금액

# 총직장인구수와 총매출금액의 산점도 그리기
ggplot(data2020, aes(x=total.work.population2020, y=total.sales2020, colour=total.work.population2020)) + geom_point() + ggtitle(label = "2020")

총 직장인구수와 총매출금액 간 상관계수가 양의 관계로 나타났고, 그 정도가 0.48 정도로 나타난 것을 보았을 때, 즉 둘의 상관관계는 약하지 않은 것을 알 수 있다.

4. 주성분 분석과 다중 회귀 분석

분석에 사용할 데이터가 준비되었으니 데이터로 부터 주성분을 추출해낸 다음 추출한 주성분을 가지고 다중 회귀 분석을 수행하였다.

그전에 먼저 변수들 중 상관성이 높은 변수들이 있을 것으로 예상되어 변수간에 상관성을 한번 확인해보았다.

4.1 상관계수행렬

names(data2020)
##  [1] "상권_코드"                   "기준_년_코드"               
##  [3] "총매출금액"                  "점포수"                     
##  [5] "총_유동인구_수"              "남성_유동인구_수"           
##  [7] "여성_유동인구_수"            "연령대_10_유동인구_수"      
##  [9] "연령대_20_유동인구_수"       "연령대_30_유동인구_수"      
## [11] "연령대_40_유동인구_수"       "연령대_50_유동인구_수"      
## [13] "연령대_60_이상_유동인구_수"  "시간대_1_유동인구_수"       
## [15] "시간대_2_유동인구_수"        "시간대_3_유동인구_수"       
## [17] "시간대_4_유동인구_수"        "시간대_5_유동인구_수"       
## [19] "시간대_6_유동인구_수"        "월요일_유동인구_수"         
## [21] "화요일_유동인구_수"          "수요일_유동인구_수"         
## [23] "목요일_유동인구_수"          "금요일_유동인구_수"         
## [25] "토요일_유동인구_수"          "일요일_유동인구_수"         
## [27] "총_직장_인구_수"             "남성_직장_인구_수"          
## [29] "여성_직장_인구_수"           "연령대_10_직장_인구_수"     
## [31] "연령대_20_직장_인구_수"      "연령대_30_직장_인구_수"     
## [33] "연령대_40_직장_인구_수"      "연령대_50_직장_인구_수"     
## [35] "연령대_60_이상_직장_인구_수" "총_상주인구_수"             
## [37] "남성_상주인구_수"            "여성_상주인구_수"           
## [39] "연령대_10_상주인구_수"       "연령대_20_상주인구_수"      
## [41] "연령대_30_상주인구_수"       "연령대_40_상주인구_수"      
## [43] "연령대_50_상주인구_수"       "연령대_60_이상_상주인구_수" 
## [45] "총_가구_수"                  "월_평균_소득_금액"          
## [47] "지출_총금액"                 "식료품_지출_총금액"         
## [49] "의류_신발_지출_총금액"       "생활용품_지출_총금액"       
## [51] "의료비_지출_총금액"          "교통_지출_총금액"           
## [53] "여가_지출_총금액"            "문화_지출_총금액"           
## [55] "교육_지출_총금액"            "유흥_지출_총금액"           
## [57] "관공서_수"                   "유치원_수"                  
## [59] "초등학교_수"                 "중학교_수"                  
## [61] "고등학교_수"                 "대학교_수"                  
## [63] "극장_수"                     "지하철_역_수"               
## [65] "버스_정거장_수"
col.length <- length(data2020)
cor <- cor(data2020[, 3:col.length], method = c("spearman"))
corrplot(cor, tl.pos='n', mar=c(0,0,1,0))

예상대로 상관성이 높은 변수들이 많다. 이대로 주성분 회귀 분석(Principal Component Regression)을 돌려보았을 때 성능이 좋지 않았기에 변수를 좀 줄여주었다.

4.2 변수 선택

변수들 중 총_xxx 와 같은 변수와 연령대별, 시간대별로 나누어져 있는 변수들을 제거해 주었다.

# 총계 변수 제거
drops <- c("총_유동인구_수",
           "총_직장_인구_수",
           "총_상주인구_수",
           "지출_총금액")
data2020 <- data2020[ ,!(names(data2020) %in% drops)]

# 연령대와 시간대 등으로 나눈 변수 제거

data2020 <- subset(data2020, select=-c(7:18, 28:33, 36:41))

주성분 회귀분석은 outlier에 취약하기에 매출액이 앞도적으로 높은 명동길을 제거하였다.(참조)

data2020 <- subset(data2020, 상권_코드!=1000039)

이대로 다시 주성분 회귀 분석을 돌려보았만 수정결정계수가 0.5256으로 낮게 나와 변수들을 좀더 줄여주었다. 이번엔 여러 변수들을 더 큰 개념으로 묶어주었다.

data2020 <- data2020 %>% mutate(학교_수 = 유치원_수 +
                                          초등학교_수 +
                                          중학교_수 +
                                          고등학교_수)
drops <- c("유치원_수", 
           "초등학교_수", 
           "중학교_수", 
           "고등학교_수")
data2020 <- data2020[ ,!(names(data2020) %in% drops)]

data2020 <- data2020 %>% mutate(대중교통_수 = 지하철_역_수 +
                                              버스_정거장_수)
drops <- c("지하철_역_수", 
           "버스_정거장_수")
data2020 <- data2020[ ,!(names(data2020) %in% drops)]

data2020 <- data2020 %>% mutate(주중_유동인구_수 = 월요일_유동인구_수 +
                                                   화요일_유동인구_수 +
                                                   수요일_유동인구_수 +
                                                   목요일_유동인구_수 +
                                                   금요일_유동인구_수)
drops <- c("월요일_유동인구_수",
           "화요일_유동인구_수",
           "수요일_유동인구_수",
           "목요일_유동인구_수",
           "금요일_유동인구_수")
data2020 <- data2020[ ,!(names(data2020) %in% drops)]

data2020 <- data2020 %>% mutate(주말_유동인구_수 = 토요일_유동인구_수 +
                                                   일요일_유동인구_수)
drops <- c("토요일_유동인구_수",
           "일요일_유동인구_수")
data2020 <- data2020[ ,!(names(data2020) %in% drops)]

4.3 모델 생성

종속변수가 정규분포를 따르게 변화시켜 분석을 실행하여 보았다. 변환에 사용한 방법은 Tuckey’s Ladder of Powers 방법으로, 이는 데이터를 최대한 정규분포로 만들어주는 lambda 값을 찾고 이를 변수에 n제곱 시켜 값을 변환시켜준다.

왼쪽으로 치우쳐진 종속변수를 Tukey’s Ladder of Powers 방법으로 데이터를 최대한 정규분포에 가깝게 변환시켜 주었다. lambda 값으로는 0 보다 큰 0.15가 나왔으므로 변수들에 0.15 제곱을 해주었다. 또한 독립변수도 종속변수가 최대한 선형관계가 되도록 종속변수에 적용된 lambda값으로 제곱하여 변환시켜주었다.

모델 생성하기 용이하게 주성분 분석과 다중 회귀 분석을 함수로 선언하여 주었다.

pca <- function (data) {
  res = list()
  dep.var <- transformTukey(data$총매출금액)
  lambda <- transformTukey(data$총매출금액, plotit = FALSE,, quiet = TRUE, returnLambda = TRUE)
  indep.var <- data[, -c(1:3)]
  indep.var <- indep.var^lambda
  pca <- prcomp(indep.var, center=TRUE, scale.=TRUE)
  print(summary(pca))
  return(list(pca, dep.var, indep.var))
}

pcr <- function (x) {
  model <- lm(총매출금액 ~ ., data = x)
  print(summary(model))
  return(model)
}

4.3.1 주성분 분석

par(mfrow=c(2,2))
res <- pca(data2020)
## 
##     lambda      W Shapiro.p.value
## 407   0.15 0.9965         0.02412
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     3.7376 1.7932 1.3029 1.08357 0.96929 0.93013 0.91206
## Proportion of Variance 0.5588 0.1286 0.0679 0.04697 0.03758 0.03461 0.03327
## Cumulative Proportion  0.5588 0.6874 0.7553 0.80227 0.83985 0.87446 0.90773
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85888 0.74357 0.65401 0.55801 0.27620 0.24551 0.22981
## Proportion of Variance 0.02951 0.02212 0.01711 0.01245 0.00305 0.00241 0.00211
## Cumulative Proportion  0.93724 0.95935 0.97646 0.98892 0.99197 0.99438 0.99649
##                           PC15   PC16   PC17    PC18    PC19    PC20    PC21
## Standard deviation     0.19119 0.1501 0.1003 0.08141 0.07417 0.05540 0.04901
## Proportion of Variance 0.00146 0.0009 0.0004 0.00027 0.00022 0.00012 0.00010
## Cumulative Proportion  0.99796 0.9989 0.9993 0.99952 0.99974 0.99987 0.99996
##                           PC22    PC23     PC24     PC25
## Standard deviation     0.02369 0.01842 0.005119 0.004086
## Proportion of Variance 0.00002 0.00001 0.000000 0.000000
## Cumulative Proportion  0.99998 1.00000 1.000000 1.000000
pca.res <- res[[1]]

fviz_pca_biplot(pca.res, repel = FALSE, alpha.ind = 0.5, col.var = "#cc3300", col.ind = "#696969", label = c("var"))

pc.cut <- 4
pca.loadings <- pca.res$rotation[,1:pc.cut]
pca.loadings
##                                 PC1          PC2         PC3           PC4
## 점포수                -0.0009814331  0.407124865 -0.33838370 -0.0184965707
## 남성_유동인구_수      -0.2232869912  0.009213452 -0.37149132 -0.1576310794
## 여성_유동인구_수      -0.2277796261 -0.013028984 -0.34937794 -0.1415889353
## 남성_직장_인구_수      0.0345530718  0.501966051  0.10021644 -0.0369369115
## 여성_직장_인구_수      0.0232056337  0.506016038  0.06361995 -0.0392896362
## 남성_상주인구_수      -0.2427870682 -0.107545847  0.11138376  0.1031397266
## 여성_상주인구_수      -0.2447276806 -0.078588658  0.10860716  0.1005930575
## 총_가구_수            -0.2464960123 -0.078137662 -0.01705515  0.0796142503
## 월_평균_소득_금액      0.0642664302  0.336597789  0.20904170 -0.0141368705
## 식료품_지출_총금액    -0.2622430133 -0.010328616  0.09003364  0.0353164684
## 의류_신발_지출_총금액 -0.2595157334  0.075876598  0.11290596  0.0172814848
## 생활용품_지출_총금액  -0.2582528015  0.090377755  0.10476343  0.0024537790
## 의료비_지출_총금액    -0.2615071962 -0.011584948  0.10221096  0.0423406134
## 교통_지출_총금액      -0.2534378551  0.096861905  0.16298622  0.0097043889
## 여가_지출_총금액      -0.2502109961  0.123931768  0.13546340 -0.0009416323
## 문화_지출_총금액      -0.2614350045  0.045157740  0.05786548  0.0179611027
## 교육_지출_총금액      -0.2462369666  0.118526186  0.15169598  0.0171347015
## 유흥_지출_총금액      -0.2595055383  0.070932536  0.10102852  0.0222220831
## 관공서_수             -0.0052744195  0.144269471 -0.18303669  0.4985229337
## 대학교_수              0.0545276939 -0.014274516 -0.31023070  0.3072338341
## 극장_수                0.0924660090  0.255448324 -0.16365433 -0.0011380308
## 학교_수               -0.0470617062 -0.169631361  0.04685296  0.3868187484
## 대중교통_수           -0.0199940204  0.103459987 -0.11928896  0.6223189775
## 주중_유동인구_수      -0.2230470864  0.017191396 -0.37308693 -0.1597998659
## 주말_유동인구_수      -0.2334995726 -0.052482385 -0.32747307 -0.1210182306

pca.loadings의 절댓값을 확인 해 보았을 때 PC1에 많은 영향을 미치는 변수는 유동인구, 상주인구, 지출금액, 총가구수와 같은 주거와 관련이 있는것들임을 알 수 있다. 따라서 PC1을 주거지역이라 할 수 있을것 같다. PC2에 많은 영향을 미치는 변수는 직장인구수로 PC2를 회사들이 모여있는 업무지구라고 할 수 있을것 같다. PC3에 많은 영향을 미치는 변수는 유동인구수로 그대로 유동인구라 할 수 있을것 같다. PC4에 많은 영향을 미치는 변수는 대중교통, 관공서, 학교로 공공기관이라 할 수 있을것 같다.

4.3.2 다중 회귀 분석

dep.var <- res[[2]]
pcs <- as.data.frame(pca.res$x)
pcs.cut <- pcs[1:pc.cut]
x <- bind_cols(dep.var, pcs.cut)
colnames(x) <- c("총매출금액", "주거지역", "업무지구", "유동인구", "공공기관")
model <- pcr(x = x)
## 
## Call:
## lm(formula = 총매출금액 ~ ., data = x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3061 -1.8423  0.0298  1.8426 10.6371 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 25.07956    0.08729 287.298   < 2e-16 ***
## 주거지역     0.11342    0.02337   4.854 0.0000014 ***
## 업무지구     1.74761    0.04871  35.881   < 2e-16 ***
## 유동인구    -1.13521    0.06703 -16.935   < 2e-16 ***
## 공공기관    -0.22060    0.08060  -2.737   0.00631 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.772 on 1003 degrees of freedom
## Multiple R-squared:  0.6155, Adjusted R-squared:  0.6139 
## F-statistic: 401.3 on 4 and 1003 DF,  p-value: < 2.2e-16

주거지역, 업무지구, 유동인구, 공공기관 변수의 p-value가 0.05이하로 종속변수인 총매출금액에 유의한 영향을 끼친다는 것을 알 수 있다.

tab_model(model,
          show.se = TRUE,
          show.ci = FALSE,
          show.stat = TRUE,
          pred.labels = c("(Intercept)", "주거지역", "업무지구", "유동인구", "공공기관"),
          dv.labels = "총매출금액")
  총매출금액
Predictors Estimates std. Error Statistic p
(Intercept) 25.08 0.09 287.30 <0.001
주거지역 0.11 0.02 4.85 <0.001
업무지구 1.75 0.05 35.88 <0.001
유동인구 -1.14 0.07 -16.94 <0.001
공공기관 -0.22 0.08 -2.74 0.006
Observations 1008
R2 / R2 adjusted 0.615 / 0.614

R-squared값은 0.615, Adjusted R-squared값은 0.614으로 종속변수의 분산 중 독립변수에 의해 설명되는 분산이 61.4%임을 알 수 있다.

shapiro.test(model$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model$residuals
## W = 0.99759, p-value = 0.1447

회귀분석의 가정중 하나인 잔차가 정규성을 지니는지 확인해보기 위해잔차의 정규성을 Shapiro-Wilk test로 알아보았다. p-value가 0.1447로 0.05보다 크게 나왔기 때문에 귀무가설을 기각하고 잔차가 정규성을 지니고 있다는 대립가설을 채택하였다.

par(mfrow=c(2,2))
plot(model)

모델의 plot을 보면 이상치들이 있는것을 알 수 있다. 이러한 이유는 아마 분석 대상이 서울 전체 상권으로 설정하였기에 상권별로 매출금액의 차이가 많을것으로 예상할 수 있다.

4.4 회귀 모형

선택한 주성분들을 독립변수로, 총매출금액을 종속변수로 설정하고 다중회귀분석을 진행하였을 때 아래와 같은 회귀식을 추정할 수 있었다.

\[Y^{0.15}=25.08+0.11*주거지역+1.75*업무지구-1.13*유동인구-0.22*공공기관\]

종속변수에 Tukey’s Ladder of Powers 변환을 해주었기 때문에 회귀식 해석 시 이러한 변환을 고려한 해석이 필요하다. 종속변수 총매출금액에 0.15를 제곱해주었기 때문에 독립변수인 주성분이 한단위 증가할 수록 총매출금액이 0.15승 증가한다고 해석할 수 있다.

5. 결론

본 분석은 서울시 골목상권별 카페업종에 대해 매출에 영향을 미치는 변수들이 무엇인지 주성분 분석을 통해 주성분을 추출해 내고 추출해낸 주성분을 독립변수로 하여 다중회귀분석을 진행하였다. 데이터들을 전처리하는 과정을 거치고 의미있는 변수들을 선택하여 62개의 변수 중 25개를 선택하였다. 이후 주성분 분석을 통해 25개의 독립변수를 주거지역, 업무지구, 유동인구, 공공기관 4개의 변수로 차원을 축소할 수 있었다. 추출한 4가지 변수중 주거지역이 전체 분산의 55.88%를 설명하며 가장 중요한 변수임을 확인하였다. 그 다음으로 업무지구(12.86%), 유동인구(6.79%), 공공기관(4.69%) 순 이었다. 이렇게 추출해낸 변수를 독립변수로 설정하고 다중회귀분석을 실시한 결과 수정된 결정계수(Adjusted R-squared)값이 0.614로 종속변수의 분산을 독립변수가 61.4% 설명할 수 있는 모델을 구축하였다. 분석 대상을 서울시 전체에서 상권으로 좁힌다면 더 좋은 모델을 구축할 수 있을것이라 예상된다.