공공데이터 with 샤이니 10_2 - 확률 밀도 함수, 회귀 분석, 주성분 분석 추가하기

1. 데이터 APP 개발하기 : 확률 밀도 함수, 회귀 분석, 주성분 분석 추가까지

 

1) 라이브러리 준비, 마커 클러스터링 까지

 

# 배포판 작성
# 10_1, 10_2, 10_3, 10_4, 10_5 종합적으로 작성

# 통계분석 8장 - 확률밀도함수, 회귀선, 주성분석 결과를 하나의 탭으로 3개의 열을 나눠서 표현하는 것을 추가!

# 1. 라이브러리 준비
library(sp);library(sf);library(leaflet);library(tidyverse);library(raster);library(rgdal);library(tmap);library(mapedit);library(leafem);library(ggfortify);library(grid);library(shiny);library(DT)
# 2. 데이터 불러오기
# 아파트 가격
load("./data/apt_price.rdata")
# 서울시 그리드 데이터
grid <- st_read("./data/sigun_grid/seoul.shp")
# 서울시 경계선 데이터
bnd <- st_read("./data/sigun_bnd/seoul.shp")
# 최고가 레스터 데이터
load("./data/raster_high.rdata")
# 급등 지역 레스터 데이터
load("./data/raster_hot.rdata")
# 3. 그리드 데이터 전처리
# 그리드 데이터 sf형을 sp형으로 변환
grid <- as(grid, "Spatial")
# sfg: simple feature geometry
# sfc: simple feature columns (the coordinate)
# sf: simple feature columns + attribute (data.frame)
# 그리드 데이터를 sfc로 변환
grid <- as(grid, "sfc")
# 아파트 가격 데이터와 그리드 데이터 합치기
# 데이터가 있는 그리드만 추출
grid <- grid[which(sapply(st_contains(st_sf(grid), apt_price),length) > 0)]
# 4. 마커 클러스터
# 마커 클러스터링 옵션 설정
# 상위 10%, 하위 10%, 나머지
# q10 <- quantile(apt_price$py, probs = seq(.1, .9, by = .1))
# 상위 10% = 90%
pct_90 <- as.numeric(quantile(apt_price$py, probs = seq(.1, .9, by = .1))[9])
# 하위 10% = 10%
pct_10 <- as.numeric(quantile(apt_price$py, probs = seq(.1, .9, by = .1))[1])
# 마커 클러스터링 함수 불러오기
load("./data/circle_marker/circle_marker.rdata")
# 색상 설정
circle.colors <- c()
#

 

2) 반응형 지도 만들기

 

# 5. 반응형 지도 만들기 => shiny application  내부
# 6. shiny application
# ui
ui <- fluidPage(
  # 상단: 좌 (반응형 지도)
  column(width = 9,
         selectModUI("selectmap"),
         div(style = "height:40px;")),
  # 상단: 우 (슬라이더-면적, 건축연도)
  column(width = 3,
         sliderInput(inputId = "range_area", # ID
                     "Area", # 제목
                     min = 0, max = 350, # 범위: 최소값, 최대값
                     value = c(0, 200)), # default 값
         sliderInput(inputId = "range_con", # ID
                     "Construction year", # 제목
                     min = 1960, max = 2021, # 범위: 최소값, 최대값
                     value = c(2000, 2021))), # default 값
  
  # 하단 : 탭 2개
  # 첫번째 탭 : 시각화 3개
  # 두번째 탭 : 입력값에 따른 데이터 테이블 형태로 출력
  
  
  tabsetPanel(
    tabPanel("Visualizations", # 탭에 노출되어 있는 이름
             column(width = 4, h4("Price", align = "center"), 
                    plotOutput("density", height = 250)), # 확률밀도함수
             column(width = 4, h4("Price Trends", align = "center"), 
                    plotOutput("regression", height = 250)), # 회귀선
             column(width = 4,h4("PCA", align = "center"),
                    plotOutput("pca", height = 250)) # 주성분분석
    ),
    tabPanel("Table", # 탭에 노출되어 있는 이름 
             dataTableOutput("table")) # output의 결과
  )
)

 

지난번과 다르게 하단 탭이 추가 됨.

  # 첫번째 탭 : 시각화 3개 (Price, Price Trends, PCA)
  # 두번째 탭 : 입력값에 따른 데이터 테이블 형태로 출력 (기존 데이터)

 

3) 반응식 설정하기 (서버)

# server
server <- function(input, output, session){
  # 반응식
  # A: 면적과 건축연도 범위에 해당하는 데이터 추출
  all <- reactive({
    all <- subset(apt_price, area >= input$range_area[1] & area <= input$range_area[2] & con_year >= input$range_con[1] & con_year <= input$range_con[2])
    return(all)})

 

4) 지도 입출력 모듈 설정하기 (서버와 이어짐)

 

  # B: 선택된 그리드 저장
  sel <- callModule(selectMod, "selectmap", 
                    leaflet() %>% 
                      addTiles(options = providerTileOptions(minZoom = 9, maxZoom = 18)) %>% 
                      # 서울시 경계선
                      addPolygons(data = bnd,
                                  color = "gold",
                                  fill = NA,
                                  weight = 5) %>% 
                      # 최고가 레스터 데이터
                      addRasterImage(raster_high,
                                     colors = colorNumeric(palette = c("blue","green","yellow","red"),
                                                           domain = values(raster_high),
                                                           na.color = "transparent"),
                                     opacity = 0.4) %>% 
                      # 급등 지역 레스터 데이터
                      addRasterImage(raster_hot,
                                     colors = colorNumeric(palette = c("blue", "green", "yellow", "red"),
                                                           domain = values(raster_hot),
                                                           na.color = "transparent"),
                                     opacity = 0.4) %>% 
                      # 체크 박스
                      addLayersControl(baseGroups = c("High Price in 2021", "Hot Price in 2021"),
                                       options = layersControlOptions(collapsed = F)) %>% 
                      # 마커 클러스터링
                      addCircleMarkers(data = apt_price,
                                       # 경도
                                       lng = unlist(map(apt_price$geometry,1)),
                                       # 위도
                                       lat = unlist(map(apt_price$geometry,2)),
                                       fillColor = circle.colors,
                                       weight = apt_price$py,
                                       clusterOptions = markerClusterOptions(iconCreateFunction=JS(avg.formula))) %>% 
                      # 그리드 추가
                      addFeatures(st_sf(grid), layerId = ~seq_len(length(grid)), color = "gray")
  )

 

5) 선택에 따른 반응 결과 저장하기 

  # A 와 B 교집합
  # 초기 반응값 저장
  rv <- reactiveValues(intersect = NULL, selectgrid = NULL)
  # 반응 결과 저장
  observe({
    # grid selected => 선택된 그리드 저장
    gs <- sel()
    # 선택(selected == TRUE)된 그리드 ID를 숫자형으로 변환하여 그리드 데이터에서 가져오기(sf형으로 변환)
    rv$selectgrid <- st_sf(grid[as.numeric(gs[which(gs$selected == TRUE),"id"])])
    # 데이터가 있는 경우
    if(length(rv$selectgrid) > 0){
      # A & B 교집합
      rv$intersect <- st_intersects(rv$selectgrid, all())
      # geometry 삭제
      # rv$intersect 결과는 리스트
      rv$sel <- st_drop_geometry(apt_price[unlist(rv$intersect), ])
      
      
    } else{ # 데이터가 없는 경우
      rv$intersect <- NULL
    }
  })

 

6) 확률밀도함수 추가

 

  # 확률밀도함수
  output$density <- renderPlot({
    # 데이터가 없는 경우
    if(nrow(rv$intersect) == 0)
      return(NULL)
    
    
    # 데이터가 있는 경우
    max_all <- density(all()$py)
    max_sel <- density(rv$sel$py)
    # x: 평당 가격
    # y: 확률 = 상대도수
    max_all <- max(max_all$y)
    max_sel <- max(max_sel$y)
    # y축 최대값 // y축에는 max_all 들어가야 함
    plot_high <- max(max_all,max_sel)
    # plot_high
    # 평균값
    avg_all <- mean(all()$py)
    avg_sel <- mean(rv$sel$py)
    # avg_all;avg_sel
    
    #
    #
    
    # 서울시 전체 지역 확률밀도함수 그리기
    plot(density(all()$py),
         ylim = c(0,plot_high),
         col = "blue",
         lwd = 2,
         main = NA
    )
    # 서울시 전체 지역 평균값 선 그리기
    abline(v = avg_all, 
           col = "blue",
           lwd = 2,
           lty = 2)
    # 평균값 텍스트 입력 
    text(avg_all + avg_all*0.2, # x축 위치
         plot_high * 0.1, # y축 위치
         round(avg_all))
    
    #
    #
    # 선택 지역 확률밀도함수 그리기
    lines(density(rv$sel$py),
          ylim = c(0,plot_high),
          col = "red",
          lwd = 2
    )
    # 서울시 선택 지역 평균값 선 그리기
    abline(v = avg_sel, 
           col = "red",
           lwd = 2,
           lty = 2)
    # 평균값 텍스트 입력 
    text(avg_sel + avg_sel*0.2, # x축 위치
         plot_high * 0.1, # y축 위치
         round(avg_sel))
    
    
    
  })

 

 

7) 회귀선 추가

  # 회귀선
  output$regression <-  renderPlot({
    # 데이터가 없는 경우
    if(nrow(rv$intersect) == 0)
      return(NULL)
    # 거래 연월일 변수의 타입을 날짜로 변환
    # str(all)
    # all$ymd <- as.Date(all$ymd, "%Y-%m-%d")
    # str(sel)
    # sel$ymd <- as.Date(sel$ymd, "%Y-%m-%d")
    
    # 월별 평균 가격
    # 서울시 전체 지역
    all <- aggregate(all()$py ~ lubridate::floor_date(all()$ymd, "month"), FUN = "mean")
    # 선택 지역
    sel <- aggregate(rv$sel$py ~ lubridate::floor_date(rv$sel$ymd, "month"), FUN = "mean")
    # 변수 이름 변경
    # head(all);head(sel)
    colnames(all) <- c("month","all_py")
    colnames(sel) <- c("month","sel_py") 
    # head(all);head(sel)
    
    
    # 회귀 모델
    # 서울시 전체 지역
    fit_all <- lm(all$all_py~ all$month)
    # 선택 지역
    fit_sel <- lm(sel$sel_py~ sel$month)
    # 회귀계수 = 기울기
    coef_all <- summary(fit_all)$coefficients[2]
    coef_sel <- summary(fit_sel)$coefficients[2]
    # 1년의 변화량
    coef_all <- round(coef_all,1)*365
    coef_sel <- round(coef_sel,1)*365
    # coef_all;coef_sel
    
    # 서울시 선택 지역 주석 추가
    # hjust(왼쪽, 오른쪽 0~1에 따라)
    # vjust(위, 아래 0~1에 따라)
    sel_grob <- grobTree(textGrob(paste0("Selected:", coef_sel),
                                  x = 0.1,
                                  y = 0.9,
                                  hjust = 0,
                                  gp = gpar(col = "red",
                                            fontsize = 14,
                                            fontface = "bold")))
    
    
    
    # 서울시 전체 지역 주석 추가
    all_grob <- grobTree(textGrob(paste0("All:", coef_all),
                                  x = 0.1,
                                  y = 0.85,
                                  hjust = 0,
                                  gp = gpar(col = "blue",
                                            fontsize = 14,
                                            fontface = "italic")))
    
    
    
    # 선택 지역 
    sel_p <- ggplot(sel, aes(x = month, y = sel_py))+
      geom_line(color = "red", size = 1) +
      xlab("Month")+
      ylab("Price")+
      stat_smooth(method = "lm",
                  color = "gray",
                  linetype = "dashed" ) +
      theme_bw()
    
    
    # 주석 추가하여 다시 시각화 # 서울시 전체 지역
    sel_p + geom_line(data = all, aes(x = month, y =all_py), color="blue", size = 1)+
      annotation_custom(sel_grob)+
      annotation_custom(all_grob)
    
    
  })

 

8) PCA 추가

 

  # PCA
  output$pca <-  renderPlot({
    # 데이터가 없는 경우
    if(nrow(rv$intersect) == 0)
      return(NULL)
    # 데이터가 있는 경우
    # 아파트 별 평균값
    pca_data <- aggregate(list(rv$sel$con_year, rv$sel$floor, rv$sel$py,         rv$sel$area), by = list(rv$sel$apt_nm),mean)
    # 변수 이름 변경
    colnames(pca_data) <- c("apt_nm","con_year","floor","py","area")
    # head(pca_data)
    # 주성분 분석
    pca <- prcomp(~ con_year + floor + py +area,
                  data = pca_data,
                  scale = T) # 표준화 = z점수 = 변동성 맞춰주는 작업
    
    # 아파트 이름 추가
    autoplot(pca,
             loadings.label = T,
             loadings.label.size = 5) +
      geom_label(aes(label=pca_data$apt_nm),
                 size=5)
    
    
  })

 

9) 반응 결과 렌더링 후 실행까지

 

  # rv$sel => 데이터테이블로 저장(렌더링 함수 사용)
  output$table <- DT::renderDataTable({
    dplyr::select(rv$sel, 
                  # 열 선택
                  ymd, addr_1, apt_nm, price, area, floor, py) %>% 
      arrange(desc(py))},
    extensions = 'Buttons',
    options = list(dom = 'Bfrtrip',
                   scrollY = 300,
                   scrollCollapse = T,
                   paging = T,
                   buttons = c('excel'))
  )
}
# shinyApp
shinyApp(ui, server)

 

Tap 1 (확률밀도함수, 회귀선, PCA)

 

Tab 2 (데이터 차트 뽑기)

 

 

 

* 핵심

- 샤이니 앱을 만드는 구조는 똑같다.

라이브러리 준비, 마커클러스터링 - 반응형 지도 - 반응식 설정 (서버) - 지도 입출력 모듈 설정
- 선택에 따른 반응 결과 저장 - (확률밀도함수, 회귀선, PCA) - 반응 결과 렌더링 후 실행까지

- 다만 원하는 분석 결과 값을 어디에 출력할 것인지에 따라 분석하는 코드의 위치만 바뀐다.