5 minute read

지리적 가중 포아송 회귀를 최종 모델로 선택했고, 2개의 사고유형, 6개의 연령대별을 합쳐서 총 12개의 그룹에 대해 각각 다른 모델을 만들었다. 또한, 이 12개의 모델의 회귀계수를 ‘coef_사고유형_연령대’ 형식의 csv로 저장했다. 이 회귀계수 데이터에는 회귀계수 뿐만 아니라 각 row별로 잔차가 포함되어있는데 이 잔차를 활용하여 추가적인 분석을 진행해보자.

library(sp)
library(dplyr)
library(stringr)
library(rgeos)
library(tmap)
library(raster)
library(spdep)
library(gstat)
library(spgwr)
library(GWmodel)
library(regclass)
library(ggplot2)
library(ggcorrplot)
library(lmtest)
library(leaflet)
library(extrafont)
acci_count_filter25 <- read.csv('accident_count_filter_25.csv')
coef_60_car <- read.csv("coef_car_over60.csv")
coef_50_car <- read.csv("coef_car_50.csv")
coef_40_car <- read.csv("coef_car_40.csv")
coef_30_car <- read.csv("coef_car_30.csv")
coef_20_car <- read.csv("coef_car_20.csv")
coef_under20_car <- read.csv("coef_car_under20.csv")

coef_60_person <- read.csv("coef_person_over60.csv")
coef_50_person <- read.csv("coef_person_50.csv")
coef_40_person <- read.csv("coef_person_40.csv")
coef_30_person <- read.csv("coef_person_30.csv")
coef_20_person <- read.csv("coef_person_20.csv")
coef_under20_person <- read.csv("coef_person_under20.csv")

Modeling Part.4(by using R)

잔차분석

모델이 위험한 지역을 위험하지 않다고 한 것과 위험하지 않은 지역을 위험하다고 한 것중 무엇이 더 중요할까? 만약 위험한 지역을 위험하지 않는다고 예측해서 교통 사고 예방에 대한 조치가 취해지지 않는다면 사람의 생명과 직접적인 연관이 있다. 반면, 위험하지 않은 지역을 위험하다고 예측해서 조치를 취한다면 경제적인 비용이 발생하긴 하지만, 그 지역은 더 위험하지 않은 지역이 될 수 있다.

결국 이와 같은 상황에서는 위험한 지역을 위험하지 않는다고 한 경우가 위험하지 않는 지역을 위험하다고 예측한 경우보다 더 중요하다. 이를 잔차에 비유하자면, 실제로 위험한 지역을 위험하지 않는다고 한 경우는 양의 잔차(예측한 사고건수보다 실제로 그 지역에서 발생한 사고건수가 큰 것)에 해당한다. 이 양의 잔차의 크기가 크다면 모델이 예측한 것보다 더 위험한 지역 즉, 대책이 필요한 지역일 가능성이 높다. 따라서 이 잔차분석을통해 양의 잔차가 큰 지역의 특성과 사고들의 유형 분석을 통해 적합한 교통 안전시설물 설치와 같은 대책을 제안해볼 수 있다.

1) GWPR로 구한 회귀계수 데이터에 고유 인덱스 부여 및 데이터 병합

도출된 회귀계수 데이터에 고유 인덱스 부여

coef_under20_car$index<-car_under20$X
coef_20_car$index<-car_20$X
coef_30_car$index<-car_30$X
coef_40_car$index<-car_40$X
coef_50_car$index<-car_50$X
coef_60_car$index<-car_over60$X
coef_under20_person$index<-person_under20$X
coef_20_person$index<-person_20$X
coef_30_person$index<-person_30$X
coef_40_person$index<-person_40$X
coef_50_person$index<-person_50$X
coef_60_person$index<-person_over60$X



각 그룹의 잔차 컬럼명 구분 작업(residual로 되어있는 변수이름에 그룹명을 추가)

colnames(coef_under20_car)[which(colnames(coef_under20_car)=="residuals")]<-"residuals_under20_car"
colnames(coef_20_car)[which(colnames(coef_20_car)=="residuals")]<-"residuals_20_car"
colnames(coef_30_car)[which(colnames(coef_30_car)=="residuals")]<-"residuals_30_car"
colnames(coef_40_car)[which(colnames(coef_40_car)=="residuals")]<-"residuals_40_car"
colnames(coef_50_car)[which(colnames(coef_50_car)=="residuals")]<-"residuals_50_car"
colnames(coef_60_car)[which(colnames(coef_60_car)=="residuals")]<-"residuals_60_car"
colnames(coef_under20_person)[which(colnames(coef_under20_person)=="residuals")]<-"residuals_under20_person"
colnames(coef_20_person)[which(colnames(coef_20_person)=="residuals")]<-"residuals_20_person"
colnames(coef_30_person)[which(colnames(coef_30_person)=="residuals")]<-"residuals_30_person"
colnames(coef_40_person)[which(colnames(coef_40_person)=="residuals")]<-"residuals_40_person"
colnames(coef_50_person)[which(colnames(coef_50_person)=="residuals")]<-"residuals_50_person"
colnames(coef_60_person)[which(colnames(coef_60_person)=="residuals")]<-"residuals_60_person"



회귀계수 데이터 병합

residual_join_dta<-coef_under20_car%>%
  dplyr::full_join(coef_20_car,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car)%>%
  dplyr::full_join(coef_30_car,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car)%>%
  dplyr::full_join(coef_40_car,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car)%>%
  dplyr::full_join(coef_50_car,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car)%>%
  dplyr::full_join(coef_60_car,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car)%>%
  dplyr::full_join(coef_under20_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person)%>%
  dplyr::full_join(coef_20_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person,residuals_20_person)%>%
  dplyr::full_join(coef_30_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person,residuals_20_person,residuals_30_person)%>%
  dplyr::full_join(coef_40_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person,residuals_20_person,residuals_30_person,residuals_40_person)%>%
  dplyr::full_join(coef_50_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person,residuals_20_person,residuals_30_person,residuals_40_person,residuals_50_person)%>%
  dplyr::full_join(coef_60_person,"index")%>%
  dplyr::select(index,residuals_under20_car,residuals_20_car,residuals_30_car,residuals_40_car,residuals_50_car,residuals_60_car,residuals_under20_person,residuals_20_person,residuals_30_person,residuals_40_person,residuals_50_person,residuals_60_person)



고유 index를 기준으로 full_join 해서 합친 데이터의 NA 부분을 0으로 바꿔준다.


residual_join_dta$residuals_under20_car<-ifelse(is.na(residual_join_dta$residuals_under20_car),0,residual_join_dta$residuals_under20_car)
residual_join_dta$residuals_20_car<-ifelse(is.na(residual_join_dta$residuals_20_car),0,residual_join_dta$residuals_20_car)
residual_join_dta$residuals_30_car<-ifelse(is.na(residual_join_dta$residuals_30_car),0,residual_join_dta$residuals_30_car)
residual_join_dta$residuals_40_car<-ifelse(is.na(residual_join_dta$residuals_40_car),0,residual_join_dta$residuals_40_car)
residual_join_dta$residuals_50_car<-ifelse(is.na(residual_join_dta$residuals_50_car),0,residual_join_dta$residuals_50_car)
residual_join_dta$residuals_60_car<-ifelse(is.na(residual_join_dta$residuals_60_car),0,residual_join_dta$residuals_60_car)
residual_join_dta$residuals_under20_person<-ifelse(is.na(residual_join_dta$residuals_under20_person),0,residual_join_dta$residuals_under20_person)
residual_join_dta$residuals_20_person<-ifelse(is.na(residual_join_dta$residuals_20_person),0,residual_join_dta$residuals_20_person)
residual_join_dta$residuals_30_person<-ifelse(is.na(residual_join_dta$residuals_30_person),0,residual_join_dta$residuals_30_person)
residual_join_dta$residuals_40_person<-ifelse(is.na(residual_join_dta$residuals_40_person),0,residual_join_dta$residuals_40_person)
residual_join_dta$residuals_50_person<-ifelse(is.na(residual_join_dta$residuals_50_person),0,residual_join_dta$residuals_50_person)
residual_join_dta$residuals_60_person<-ifelse(is.na(residual_join_dta$residuals_60_person),0,residual_join_dta$residuals_60_person)



DGI(Daejeon Gid Index = 대전격자지수)

image

위와 같은 수식처럼 각 그룹별로 도출된 가중치 데이터와 잔차를 활용하여 대전광역시 격자별 위험지수를 계산한다.

DGI지수 를 구하는 함수

# DGI(Daejeon Gid Index-대전격자지수)
DGI<-function(weightsvec=NULL,residualsmat=NULL,DIratiovec,threshold_quantiles=NULL,gid=NULL,gid_include=FALSE){
  #residualsmat=>(r1vec,r2vec,r3vec,...,rnvec)
  r_mat<-apply(residualsmat,2,function(x){ifelse(x<=0,0,x)})
  DGIj<-r_mat%*%weightsvec
  if(gid_include){
    DGIj<-as.data.frame(cbind(DGIj,DIratiovec,1:nrow(residualsmat)))
    DGIj$gid<-gid
    colnames(DGIj)<-c("DGI","DI_ratio","index","gid")
  }
  else{
    DGIj<-as.data.frame(cbind(DGIj,DIratiovec,1:nrow(residualsmat)))
    colnames(DGIj)<-c("DGI","DI_ratio","index")
  }
  DGIj_sorted<-DGIj[order(DGIj$DGI,decreasing = T),]
  
  DGIj_sorted$difference<-c(0,diff(DGIj_sorted$DGI))
  DGIj_sorted$false_rank<-1:nrow(residualsmat)
  threshold<-quantile(-diff(DGIj_sorted$DGI),threshold_quantiles)
  s_thres<-0
  switching<-0
  for(i in 1:(nrow(residualsmat)-1) ){
    if((DGIj_sorted[i,"DGI"]-DGIj_sorted[i+1,"DGI"])< threshold){
      if(DGIj_sorted$DI_ratio[i]<DGIj_sorted$DI_ratio[i+1]){
        cache<-DGIj_sorted[i,]
        DGIj_sorted[i,]<-DGIj_sorted[i+1,]
        DGIj_sorted[i+1,]<-cache
        switching<-switching+1
      }
      s_thres<-s_thres+1
    }
    
  }
  DGIj_sorted$true_rank<-1:nrow(residualsmat)  
  return(list(DGI=DGIj_sorted,iterations_of_switching=switching,s_thres=s_thres,weights=weightsvec,residuals_matrix=residualsmat))
}



격자 데이터(acci_count_filter25)의 index를 기준으로 각 그룹별로 잔차를 병합한 잔차 데이터(residual_join_dta)와 병합한 후 DGI 지수를 계산하기 위해 필요한 변수들만 추출한다.

finalie_last_2<-acci_count_filter25%>%
  mutate(중상자이상수=(사망자수+중상자수) )%>%
  right_join(residual_join_dta,by=c("X"="index"))%>%
  dplyr::select(gid,x,y,중상자이상수,residuals_under20_car:residuals_60_person)

finalie_last_2$index<-as.numeric(rownames(finalie_last_2))
finalie_last_real<-finalie_last_2%>%
  left_join(residual_join_dta,by="index")%>%
  dplyr::select(gid,x,y,중상자이상수,residuals_under20_car.x,residuals_20_car.x,residuals_30_car.x,residuals_40_car.x,residuals_50_car.x,residuals_60_car.x,residuals_under20_person.x,residuals_20_person.x,residuals_30_person.x,residuals_40_person.x,residuals_50_person.x,residuals_60_person.x)



DGI 지수 구하기 and 최종위험지역 100개 도출

#weightsvec->전처리 파트에서 구한 그룹별 가중치 벡터
#residualsmat->잔차 행렬(행:관측치,열:사고유형-연령대 그룹,값:잔차)
#DIratiovec->각 격자가 가지는 사상자 수에 대한 사망자 수 + 중상자 수의 비율
#threshold_quantiles->내림차순으로 정렬한 격자별 DGI 지수의 차이(절댓값)의 하위 n % 백분위수를 threshold로 설정하는데, 여기서 필요한 n 설정   
#함수 리턴 결과가 gid를 포함한 결과이고 싶을 때 gid_include=TRUE 하고 gid 값을 설정해 주면 된다.

dgi<-DGI(weightsvec=c(0.07748547,0.07605981,0.07829668,0.08305307,0.08640373,0.09488054,0.07294598,0.06206129,0.07528174,0.08830434,0.09012788,0.11509948),residualsmat=as.matrix(dplyr::select(finalie_last_2,residuals_under20_car:residuals_60_person),nrow=5138,ncol=12),DIratiovec=finalie_last_real$`중상자이상수`,threshold_quantiles=0.996,gid=finalie_last_2$gid,gid_include = TRUE)
top_100<-dgi$DGI[which(dgi$DGI$true_rank<=100),] # DGI 지수가 가장 높은 즉, 위험지역 100개지역 도출
write.csv(top_100_visu, "top_100_visu.csv", row.names = FALSE)

top_100_visu<-top_100%>% 
  left_join(finalie_last_2,by = c("gid","index"))%>%
  left_join(acci_count_filter25,by=c("gid","x","y"))%>%
  dplyr::select(DGI:residuals_60_person,acci_cnt:중상자수)

plot(finalie_last_real$x,finalie_last_real$y,col="green",pch=19)
lines(top_100_visu$x,top_100_visu$y,type="p",pch=19,col="red")

image



최종 위험지역 100개 데이터 전처리

위도, 경도 좌표를 활용하여 정확한 주소를 입력하고 반경 범위를 표시한다.

# gid 격자의 시군구 정보 추가
top_100_visu <- read.csv('top_100_visu.csv')
accident_list <- read.csv('1.대전광역시_교통사고내역(2017~2019).csv')
gid_df <- accident_list[, c('시군구', 'gid')]

final_df <- top_100_visu %>% semi_join(accident_list, by='gid') # 기존 데이터와 병합
dim(final_df)
final_df %>% head(3)

image



head(final_df[, c('y', 'x')])

image



각 격자의 위도, 경도 좌표를 바탕으로 세부적인 주소를 네이버 지도를 통해서 획득

gid_location = c('대전광역시 서구 갈마동 706', '대전광역시 서구 둔산동 1077', '대전광역시 동구 삼성동 458', 
                 '대전광역시 서구 갈마동 307-1', '대전광역시 유성구 봉명동 1063', '대전광역시 서구 관저동 616-10', 
                 '대전광역시 서구 탄방동 544', '대전광역시 유성구 원내동 360-27', '대전광역시 유성구 화암동 216-5', 
                 '대전광역시 서구 둔산동 959-1', '대전광역시 동구 대동 400-1', '대전광역시 대덕구 중리동 450', 
                 '대전광역시 서구 둔산동 911', '대전광역시 서구 둔산동 1204', '대전광역시 동구 판암동 467-13', 
                 '대전광역시 중구 태평동 520-1', '대전광역시 서구 갈마동 1459-1', '대전광역시 유성구 반석동 685', 
                 '대전광역시 동구 가양동 424-6', '대전광역시 유성구 도룡동 8-59', '대전광역시 중구 선화동 383', 
                 '대전광역시 중구 은행동 154', '대전광역시 중구 오류동 165-15', '대전광역시 유성구 봉산동 1002', 
                 '대전광역시 유성구 봉명동 1058', '대전광역시 유성구 장대동 117-13', '대전광역시 중구 대흥동 233-2', 
                 '대전광역시 서구 둔산동 2166', '대전광역시 유성구 구암동 527-61', '대전광역시 중구 대흥동 201', 
                 '대전광역시 서구 둔산동 1162', '대전광역시 동구 신흥동 6', '대전광역시 서구 갈마동 705', 
                 '대전광역시 서구 월평동 282-1', '대전광역시 중구 오류동 196-1', '대전광역시 유성구 지족동 1005', 
                 '대전광역시 유성구 구암동 641', '대전광역시 서구 변동 12-10', '대전광역시 동구 홍도동 838', 
                 '대전광역시 대덕구 중리동 469', '대전광역시 유성구 원내동 711', '대전광역시 서구 월평동 1636', 
                 '대전광역시 유성고 도룡동 465-27', '대전광역시 서구 갈마동 986', '대전광역시 유성구 봉명동 539-1', 
                 '대전광역시 대덕구 읍내동 505-23', '대전광역시 서구 월평동 184-2', '대전광역시 서구 갈마동 1459', 
                 '대전광역시 동구 중동 75-8', '대전광역시 중구 대흥동 508-88', '대전광역시 중구 문화동 1-29', 
                 '대전광역시 유성구 노은동 132', '대전광역시 서구 둔산동 948', '대전광역시 서구 갈마동 685', 
                 '대전광역시 서구 만년동 682', '대전광역시 동구 정동 45-1', '대전광역시 서구 도안동 1493-1', 
                 '대전광역시 서구 가장동 202', '대전광역시 서구 탄방동 544', '대전광역시 대덕구 오정동 43', 
                 '대전광역시 동구 용전동 220', '대전광역시 서구 변동 292', '대전광역시 유성구 노은동 61-13', 
                 '대전광역시 서구 둔산동 1546', '대전광역시 서구 월평동 1503-1', '대전광역시 서구 복수동 611', 
                 '대전광역시 유성구 봉명동 574', '대전광역시 유성구 송강동 171-1', '대전광역시 서구 둔산동 1100', 
                 '대전광역시 서구 둔산동 950', '대전광역시 서구 월평동 1528', '대전광역시 동구 용전동 227', 
                 '대전광역시 대덕구 중리동 150-3', '대전광역시 유성구 구암동 593-2', '대전광역시 유성구 봉명동 1058', 
                 '대전광역시 대덕구 오정동 64-18', '대전광역시 동구 가양동 630', '대전광역시 중구 선화동 862-3', 
                 '대전광역시 대덕구 중리동 424', '대전광역시 중구 용두동 115-2', '대전광역시 서구 도마동 227', 
                 '대전광역시 동구 자양동 67-17', '대전광역시 대덕구 비래동 103-4', '대전광역시 동구 인동 352', 
                 '대전광역시 유성구 구암동 95-5', '대전광역시 유성구 장대동 115-1', '대전광역시 동구 용전동 68-2', 
                 '대전광역시 유성구 반석동 685', '대전광역시 중구 용두동 112-10', '대전광역시 동구 원동 60-10', 
                 '대전광역시 중구 산성동 48-4', '대전광역시 유성구 지족동 1005', '대전광역시 중구 태평동 531', 
                 '대전광역시 동구 성남동 509', '대전광역시 유성구 봉명동 563-23', '대전광역시 대덕구 신탄진동 144-1',  
                 '대전광역시 대덕구 중리동 126-4', '대전광역시 중구 유천동 185-3', '대전광역시 유성구 궁동 220-9', '대전광역시 유성구 봉명동 449-4')



주소, 반경, 순위 입력

danger_df <- final_df[, c('x', 'y')]
danger_df$'시설명/주소지' <- gid_location
danger_df$'위험순위' <- c(1:100)
danger_df$'반경범위' <- '100M'
dim(danger_df)
danger_df %>% head()

image



danger_df <- danger_df[, c('위험순위', '시설명/주소지', 'x', 'y', '반경범위')]
dim(danger_df)
danger_df %>% head()

image



write.csv(danger_df, "danger_top_100.csv", row.names = F)



이로써 전처리, 모델링, 잔차분석이라는 긴 과정을 거쳐서 최종 위험지역 100개 지역을 도출했다. 다음 포스팅에서는 이 100개의 교통사고 위험지역에 어떠한 조치를 취하면 좋을지 시각화와 함께 대안책을 정리해보자.

image

Categories:

Updated: