[경제자료분석] 2019-2학기 경제자료분석

건국대학교 최고의 학과~ 응용통계학과 전공 과목인 유규상 교수님의 경제자료분석 팀프로젝트로 진행한 데이터 분석입니다.

경제자료분석 과목은 회귀분석2라고 할 수 있을 정도로 1학기 회귀분석 과목보다 확장되고 다양한 방법론, 그리고 특히 경제자료를 중점으로 어떻게 분석할 수있나에 대해서 배우는 수업이었습니다.

프로젝트는 4-5인 1조, 자유주제, 자유방법, 모든 언어 가능이었고 kick off 발표와 final 발표, 그리고 중간에 교수님과의 면담 2회 이상으로 진행되었습니다. 

프로젝트의 주요 채점 기준은 세련되거나 발전된 분석기법의 활용보다는 수업내용의 활용과 충실 등이라고 먼저 공지하셨습니다. 

제가 한 분석의 주요 idea는 '버거지수' 에서 영감 받은 '편의점지수' 였으며, improve해야하는 점은 마지막에 적어놓겠습니다~!

혹시 code를 검색하다가 유입되셨거나 관련 code를 찾고 계시면 하단부에 포스팅 되어 있으니 Ctrl + F로 찾아주세요

---

<Final PPT> 위주로 설명

이 예 진 ( 응통 17 )
Te am L e ade r ,
P r e p r o c e s s i n g , A n a l y s i s , PPT ma k i n g

정 재 윤 ( 응통 17 )
Da t a An a l y s i s ,
P r e p r o c e s s i n g , V i s u a l i z a t i o n

-> 데이터 출처를 보시면 크롤링과 공시 자료로 이루어져있음을 보실 수 있습니다.

 

 

 

 

 

 

 

 

---

PPT는 43슬라이드이며, 포스팅에 포함하지 않은 최종 보고서는 39페이지 입니다.

슬라이드와 보고서가 매우 자세하게 적혀있기 때문에 읽어보시고 궁금하신 점이나 조언 해주실 사항은 댓글 부탁드립니다!

---

<Code> R code로 진행 (R studio)

<Code>

1.1  참고 코드

# csv file to open

# make dataframe "total"

setwd=("C:/Users/jinji/Desktop/2019년/15.경제자료분석/데이터자료")

 

storedist = read.csv("C:/Users/jinji/Desktop/2019년/15.경제자료분석/데이터자료/storedist2.csv",head=T)

 

variable = read.csv("C:/Users/jinji/Desktop/2019년/15.경제자료분석/데이터자료/variable4.csv", head=T)

 

total = merge(storedist,variable,key=c("시","구"))

as.data.frame(total)

str(total)

 

storedist

variable

 

 

names(total) <- c("City","Section","CU","GS","SE","EM","MS","Latitude","Longitude","Avg_age","Avg_m","Avg_f","UnivNum","SchNum","TotalNum","MaleNum","FemaleNum","Nat","Foreign","GenRate","Area","Density","Solo","SubNum","CityRate","Theater","Hospital","SB","Smoking","Par",'Tou',"Bread")

total

str(total)

summary(total)

 

# Modified variable name

# Section

 

# make CIC Formula

## CIC = (SE + EM + MS)/(CU*0.5 + GS*0.5)

 

 

## transform()

 

 

total = transform(total,CIC = (SE + EM + MS)/(CU*0.5 + GS*0.5))

total = transform(total,CIC_v =(CU*0.5 + GS*0.5)/(SE + EM + MS))

total = transform(total,CIC_2 = (SE^2 + EM^2 + MS^2)/(CU^2*0.5 + GS^2*0.5))

total = transform(total,CIC_CR = CIC*CityRate)

total = transform(total,CIC_w = (CU+GS+SE+EM+MS)*CityRate)

total = transform(total,CIC_test = (EM + MS)/(CU*(1/3) + GS*(1/3) + SE*(1/3)))

 

str(total)

str(storedist)

 

# between variable (storedist), correlation, plot

 

plot(storedist[3:7])

cor(storedist[3:7])

 

## between store brand, high Multicollinearity

 

plot(total[20:25])

 

plot_test = data.frame(total$CIC,total$Avg_age,total$UnivNum,total$SchNum,total$Density,total$GenRate,total$Foreign,total$Solo,total$SubNum,total$Hospital,total$SB,total$Smoking,total$Bread)

plot(plot_test)

 

plot(total$CIC,total$CityRate)

cor(total$CIC,total$CityRate)

 

# plot (인구 수, 학교 수), (1인 가구, 학교 수), (1인 가구, 대학교 수), (학교 수, 대학교 수)

plot(total$TotalNum,total$SchNum)

plot(total$Solo,total$SchNum)

plot(total$Solo,total$UnivNum)

plot(total$SchNum,total$UnivNum)

 

par(mfrow = c(1,1))

 

 

cor(total$TotalNum,total$SchNum)

cor(total$Solo,total$SchNum)

cor(total$Solo,total$UnivNum)

cor(total$SchNum,total$UnivNum)

 

 

 

# between variable (total), correlation, plot

 

# Just try model, dont care about stepwise~,slect variable~

 

# model 1 : Use CIC Formula

summary(lm(CIC~Avg_age+UnivNum+SchNum+TotalNum+Foreign+GenRate+Area+Density+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

             ,data=total))

 

## Multiple R-squared:   0.1231

 

summary(lm(CityRate~Avg_age+UnivNum+SchNum+TotalNum+Foreign+GenRate+Area+Density+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total))

## Multiple R-squared:  0.7575

 

# model 2 : Use CIC_CR

 

model_CR = lm(CIC_CR~Avg_age+UnivNum+SchNum+TotalNum+Foreign+GenRate+Area+Density+Solo+SubNum            +Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total)

 

summary(model_CR)

 

library(car)

vif(model_CR)

 

step(model_CR,direction="both")

 

 

 

## Multiple R-squared:  0.3291

 

 

## After Professor Interview,,

 

# Number of convenience stores per unit area

# Number of convenience stores per population

 

 

 

# model 3 : Use Number of convenience stores per population

## (CU + GS + SE +EM + MS)/ Density to response

 

summary(lm((CU + GS + SE +EM + MS)/ Density~Avg_age+UnivNum+SchNum+Foreign+GenRate +Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total))

 

# model 4 : Use Number of convenience stores per unit area

## (CU + GS + SE +EM + MS)/ Area to response

 

summary(lm((CU + GS + SE +EM + MS)/ Area~Avg_age+UnivNum+SchNum+TotalNum+Foreign+GenRate+Density+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou,data=total))

 

density_test = data.frame(total$Density,total$SchNum,total$SubNum,total$Theater,total$SB,total$Tou)

cor(density_test)

 

 

# model4 is our final model

## stepwise for final model

 

model4 = lm((CU + GS + SE +EM + MS)/ Area~Avg_age+Avg_m+UnivNum+SchNum+TotalNum+Foreign+GenRate

            +Density+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

            ,data=total)

summary(model4)

 

model4_modified = step(model4, direction="both")

summary(model4_modified)

 

#########################################################################

 

# fimal model formula = (CU + GS + SE + EM + MS)/Area ~ Avg_age + Avg_m + SchNum + TotalNum + Foreign + Density + Solo + SubNum + Smoking)

 

# error test of model4

 

##0. 회귀계수들 간의 종속성 검정

library(car)

vif(model4_modified)

 

##1. 오차항의 독립성 검정

library(lmtest)

dwtest(model4_modified)

 

##2. 오차항의 정규성 검정

shapiro.test(resid(model4_modified))

 

##3. 오차항의 등분산성 검정

par(mfrow=c(2,2))

plot(model4_modified)

par(mfrow=c(1,1))

plot(resid(model4_modified))

 

##4. find outlier and leverage

### find outlier

outlierTest(model4_modified)

total[123,]

total[148,]

### leverage -> cooks distance

cooks.distance(model4_modified)

summary(cooks.distance(model4_modified))

plot(cooks.distance(model4_modified))

 

cooks.distance(model4_modified)[cooks.distance(model4_modified)>1]

total[148,]

 

leveragePlot(model4)

#########################################################################

# model appendix : number each convenience store to response

 

CU_lm = lm(CU~Avg_age+UnivNum+SchNum+Foreign+GenRate+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total)

summary(CU_lm)

step(CU_lm,direction = 'both')

plot(CU_lm)

 

 

GS_lm = lm(GS~Avg_age+UnivNum+SchNum+Foreign+GenRate+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

                   ,data=total)

summary(GS_lm)

step(GS_lm,direction = 'both')

plot(GS_lm)

 

 

SE_lm = lm(SE~Avg_age+UnivNum+SchNum+Foreign+GenRate+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total)

summary(SE_lm)

plot(SE_lm)

 

 

EM_lm = lm(EM~Avg_age+UnivNum+SchNum+Foreign+GenRate+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total)

summary(EM_lm)

step(EM_lm,direction = 'both')

plot(EM_lm)

 

 

MS_lm = lm(MS~Avg_age+UnivNum+SchNum+Foreign+GenRate+Solo+SubNum+Theater+Hospital+SB+Smoking+Par+Tou

           ,data=total)

summary(MS_lm)

step(MS_lm,direction = 'both')

plot(MS_lm)

 

 

#################################################################

#Appendix

 

test_lm = lm(CU~GS + SE,data=total)

summary(test_lm)

step(test_lm,direction = 'both')

plot(test_lm)

library(car)

library(lmtest)

vif(test_lm)

shapiro.test(resid(test_lm))

dwtest(test_lm)

outlierTest(test_lm)

cooks.distance(test_lm)[cooks.distance(test_lm)>1]

 

 

summary(lm(CU~GS + SE + EM + MS,data=total))

summary(lm(GS~CU + SE + EM + MS,data=total))

summary(lm(SE~CU + GS + EM + MS,data=total))

summary(lm(EM~CU + GS + SE + MS,data=total))

summary(lm(MS~CU + GS + SE + EM,data=total))

 

vif(lm(CU~GS + SE + EM + MS,data=total))

vif(lm(GS~CU + SE + EM + MS,data=total))

vif(lm(SE~CU + GS + EM + MS,data=total))

vif(lm(EM~CU + GS + SE + MS,data=total))

vif(lm(MS~CU + GS + SE + EM,data=total))

 

Appendix 참고 코드

####시각화

variable=read.csv("C:/Users/정재윤/Desktop/재윤/건국대/경제자료분석/자료종합1(지역코드추가).csv",header=T)

 

City=variable["시"]

Section=variable["구"]

CU=variable["CU"]

GS=variable["GS25"]

SE=variable["세븐일레븐"]

EM=variable["이마트24"]

MS=variable["미니스톱"]

Avg_age=variable["총평균연령"]

Avg_m=variable["남평균연령"]

Avg_f=variable["여평균연령"]

UnivNum=variable["대학"]

SchNum=variable["초중고"]

TotalNum=variable["총인구"]

MaleNum=variable["남자"]

FemaleNum=variable["여자"]

Nat=variable["내국인"]

Foreign=variable["외국인"]

GenRate=variable["성비"]

Area=variable["면적"]

Density=variable["인구밀도"]

Solo=variable["X1인가구"]

SubNum=variable["전철역수"]

Urban_level=variable["도시화율"]

lat=variable["위도"]

long=variable["경도"]

Hospital=variable["X3차.병원"]

Theater=variable["X3대.영화관"]

ST=variable["스타벅스"]

Smoking=variable["흡연율"]

PAR=variable["파리바게트"]

TOU=variable["뚜레쥬르"]

Bread=variable["빵집"]

id=variable["id"]

code=variable["code"]

 

CIC=(SE+EM+MS)/(0.5*CU+0.5*GS)

CIC_CR=CIC*Urban_level

CIC_Area=(GS+CU+SE+EM+MS)/Area

CIC_Density=(GS+CU+SE+EM+MS)/Density

variable2=data.frame(City,Section,id,CU,GS,SE,EM,MS,lat,long,Theater,Hospital,ST,Smoking,PAR,TOU,Bread,code,Avg_age,Avg_m,Avg_f,UnivNum,SchNum,TotalNum,MaleNum,FemaleNum,Nat,Foreign,GenRate,Area,Density,Solo,SubNum,Urban_level,CIC,CIC_CR,CIC_Area,CIC_Density)

attach(variable2)

names(variable2)=c("City","Section","id","CU","GS","SE","EM","MS","lat","long","Theater","Hospital","ST","Smoking","PAR","TOU","Bread","code","Avg_age","Avg_m","Avg_f","UnivNum","SchNum","TotalNum","MaleNum","FemaleNum","Nat","Foreign","GenRate","Area","Density","Solo","SubNum","Urban_level","CIC","CIC_CR","CIC_Area","CIC_Density")

 

install.packages("devtools")

devtools::install_github("cardiomoon/Kormaps2014")

devtools::install_github("cardiomoon/moonBook2")

library(kormaps2014)

library(moonBook2)

 

variable3=merge(kormap2,variable2,by="id",all.x=T)

 

areacode

areacode1=changeCode(areacode)

areacode1

 

install.packages("ggplot2")

library(ggplot2)

par(mfrow=c(1,1))

ggplot(variable3,aes(map_id=code.x,fill=CU))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 CU분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=GS))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 GS25분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=SE))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 세븐일레븐 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=EM))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 Emart 24 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=MS))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 Ministop 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=CIC))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 CIC 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=Urban_level))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 도시화율 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=CIC_CR))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 CIC_CR 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=CIC_Area))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 CIC_Area분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=Foreign))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 외국인 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=Smoking))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 흡연율 분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=CIC_Density))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 CIC_Density분포도")+

  coord_map()

 

ggplot(variable3,aes(map_id=code.x,fill=ST))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 스타벅스 분포도")+

  coord_map()

 

 

ggplot(variable3,aes(map_id=code.x,fill=GenRate))+

  geom_map(map=kormap2,colour="black",size=0.1)+

  expand_limits(x=kormap2$long,y=kormap2$lat)+

  scale_fill_gradientn(colours=c('white','orange','red'))+

  ggtitle("시도별 성비 분포도")+

  coord_map()

 

---

<Improve>

분석 후 시간이 흐른 후 지금 한계점과 분석 방법에 대해서 다시 생각해보면, 놀랍게도 그 때는 안보였던 개선할 점이 많이 보입니다.

위에 언급했다시피 제한된 data와 상관성, 그리고 회귀분석을 이용한 점 등등의 관점으로 봤을 때는 나름 최선을 다한 회귀분석 과제였지만, 모델 선택에 대해서 R^2 가 아니라 AIC를 더 적극적으로 활용하고, R^2와 adj-R^2는 선형성 설명을 위주로 분석해야했다는 아쉬움이 강하게 남습니다.

지금 듣고 수강 중인 권성훈 교수님의 data mining 수업에서 배운 내용을 토대로 개선점 이유를 설명하면 다음과 같습니다. 아래와 같은 진행으로 variable selection과 model 비교를 진행했으면 10% 더 훌륭한 프로젝트가 되었을 듯 합니다 ㅎㅎ 

 

<Variable selection for multiple linear regression>

• 변수 선택의 과정은 1. 모형 만들기 2. 모형 비교하기 로 진행됩니다.
• 이 때 비교를 원하는 후보 모형 model군을 만들어서 비교하게 되는데 2^p개를 모두 비교하기에는 무리가 있기 때문에 전체를 해보는 방법과 패널티를 이용하는 방법으로 나뉩니다. 
• -> 여기서 전체를 해보는 방법(BBS)이 제가 사용한 forward, backward, stepwise selection이고 실제로 stepwise가 세 가지 중 가장 개선된 방법이기에 이전 회귀분석 과제와 경자분 프로젝트에서 사용했으나, 위의 프로젝트 만큼의 변수 갯수가 되어도 굉장히 많은 모델들을 비교해야했습니다. ( ex 변수 10개 -> 모델 2^10 개)
• 그리고 패널티를 이용하는 방법(Penalized estimation, PE)이 한번쯤 들어봤을법한 LASSO, SCAD, MCP를 이용한 모델 선택 방법입니다. 
• (덧붙이자면, best subset 방법은 true model S*을 포함하는지 모르지만, LASSO 등의 PE 방법은 true model S*을 포함하는 것이 잘 알려져 있습니다.) 
• 추가적으로 모델 평가는 분포 가정을 통한 비교 (통계학에서 많이 사용) 와 predictive approaches 가 주로 사용됩니다. (ML 에서 우리는 accuracy를 높이고 loss를 줄이려고 갖은 노력을 하죠...) 분포 가정을 통한 비교가 GIC 입니다. GIC에 널리 알려진 AIC, BIC 등이 있습니다.

 

그래서 결과적으로 더 나은 변수선택을 통해서 최종 model을 어떻게 선택했어야 하나?!

-> 큰 흐름은 1번 모형 만들기와 2번 모형 비교하기를 염두하고 진행해야합니다.

-> 따라서 model subset 을 만들 때 전체를 만들기로 했다. 그러면 stepwise 방법(모형 만들기)을 사용하면서 AIC(모형비교하기) 를 이용해서 모델들을 비교 후에 결정했어야합니다. ( 이 때 R^2은 단일 모델의 분석 설명을 위한 수치이기 때문에 언급하면 안되고, 여러 모델 비교시에 이용되는 AIC, BIC 등의 GIC를 근거로 설명해야합니다.) 

관련 논문, [ Yang 2005 that proves confiction between selection and prediction ] 

-> 다른 방법으로 model subset 을 만들 때 PE를 사용하기로 했다. 그러면 LASSO 등의 패널티(모형 만들기)를 이용해서 모델을 만들어서 AIC(모형 비교하기) 등을 이용해서 모델들을 비교했어야합니다. 

 

---

전체적으로 제 프로젝트 과정은 이 부분이 빈약합니다. 그 동안 code로 생각 없이 진행해왔던 분석들에 대해서 data mining 수업을 들으면서 통계학적 데이터 분석 방법에 대해서 논리적, 이론적, 학문적으로 정리할 시간이 필요하기 때문에 추후에 공부하면서 더 포스팅하겠습니다!

읽어보시고 궁금하신 점이나 조언 해주실 사항은 언제든지 댓글 부탁드립니다!