다수의 회귀분석을 들어 갈 것인데, 다수의 회귀 분석에 들어가기 앞서서, 데이터 정리를 앞시간에 하였다. 앞 글에서 보았던 내용을 미리, 코딩 한 다음에 와야 아래의 내용을 이어 나갈 수 있다.
회귀분석을 다수로 할 수 있는 아파트 단지가 있는데, 2017.6~ 2021.3 까지 아파트 거래가 7건이상인 단지는 4,404이다. 이것을 전부 회귀 모형을 돌려서 R제곱 , Y 절편, 기울기를 한꺼번에 구할 것이다.
이것을 구할 수 있는 방법은 for 루프문, lapply, data.table등 여러가지가 있다. 물론 멀티 연산도 있다. 하지만, 이번에는 purrr을 이용해서 할 것이다.
다수의 회귀 분석도 전처리 분야 중에서 하나이다. 즉 어떤 동네에 있는 아파트가 가격을 추세적으로 잘 설명할 수 있는가 인데, 일반적으로 R2으로 설명을 한다.
R을 진입 하는 시기에는 항상 시계열 그래프 그리는데, 한참이나 헤메었다. 이유는 데이터 속성 때문인데, 데이터 타입이 Date Type 인가? character Type 인가? 에 따라서 그래프 그리는 방법이 달랐기 때문이다.
http://rdmkyg.iptime.org/index.php/s/Oex9srbPk7bQtGG/download
아래 그림과 같이 아파트 시군구가 묶여져 있다. 동별로 묶어서 데이터를 분석 하는데, 문제는 없으나, 시단위로 집계 할 수 없고, 구단위로 집게 할 수 없다. 그럼 아래의 데이터는 어떻게 하는 것이 좋은가?
시스템에서 일반적으로 데이터를 받으면 항상 겪는 문제가 있다. 일자, 월을 받는 것인데, 일자나 월은 1자리수도 있고 2자리수도 있다. 이것을 데이터로 묶어서, 년월일을 8자리로 바꾸거나 아니면 2021-05-10 처음 년월일 사이에 "-" 하이픈을 넣은 다음 Date 형태로 데이터 유형을 바꾸는 경우가 있다.
년월일 8자리로 묶었을 때의 장점은 데이터 소팅이 가능 하다는 것이다. 기본적으로 데이터 소팅은 자리수가 같아야 할 수 있다. 자리수가 다를 경우 아래와 같이 소팅이 잘못 될 수도 있다.
2015년 1월 11일과 2015년 2월 2일의 소팅 값이 바뀌었고, 눈으로 식별 하기가 어렵다.
> a <- c(2015111, 201522, 20151001)
> sort(a)
201522 2015111 20151001
지금 까지 친 소스코드를 파이프연산자 %>%를 활용 해서, 간단하게 정리 하였고, 계약일자의 자리수가 문제가 있다.
# 라이브러리 불러오기
library(rio)
library(dplyr)library(readxl)
library(httr)
# 파일을 불러오는 path를 url 변수로 정의
url = "https://drive.google.com/u/0/uc?id=1PT_0-xJvM9KGp5zzx1jjcPAaS25epx4l&export=download"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
apt <- read_excel(tf, 1L) # 1L은 sheet 1을 이야기 하는 것이다.
# 아파트 데이터 전처리
apt <- apt %>%
rename(전용면적_m3 = '전용면적(㎡)') %>%
rename(거래금액_만원 = '거래금액(만원)') %>%
mutate(건축년도 = as.integer(건축년도),
층 = as.integer(층)) %>%
mutate(전용면적_m3 = as.double(전용면적_m3)) %>%
mutate(거래금액_만원 = as.integer(gsub(",", "", 거래금액_만원)))
glimpse(apt %>% select(계약년월, 계약일)) ## Rows: 296,656 ## Columns: 2 ## $ 계약년월 <chr> "202103", "202103", "202103", "202103", "202103", "202103", "2021… ## $계약일 <chr> "31", "11", "15", "26", "31", "2", "9", "27", "9", "10", "2", "1…
# 문자열을 카운트 한 다음 앞에 0일 붙임
apt <- apt %>%
mutate(계약일 = ifelse(nchar(계약일) ==1,
paste( "0", 계약일, sep = ""),
계약일 ))
apt$계약일[1:10] ## [1] "31" "11" "15" "26" "31" "02" "09" "27" "09" "10"자리수가 바뀌었다. 자리수가 바뀌었다고 모두 끝나는 것이 아니다. 계약일자로 8자리 charactor로 만들어야 한다.
# 일자 자릿수를 8자리를 만든다.
apt <- apt %>%
mutate(계약일자 = paste0(계약년월, 계약일))
apt$계약일자[1:10] ## [1] "20210331" "20210311" "20210315" "20210326" "20210331" "20210302"
## [7] "20210309" "20210327" "20210309" "20210310"여기에서 데이터 관리의 노하우가 나오는데, 예를 들면, 아파트 거래시 같은 일자에 여러번 계약되는 경우가 있다. 이를 경우에는 데이터를 그룹화 하여 평당 가격을 평균으로 구한다.
# 강남구 논현동 계약일자별 평당 가격
apt %>%
filter(시군구 == "서울특별시 강남구 논현동") %>%
mutate(평당가격 = 거래금액_만원 / (전용면적_m3/3.3) ) %>%
group_by(시군구,계약일자) %>%
summarise(평당가격_일별= mean(평당가격))## # A tibble: 485 x 3
## # Groups: 시군구 [1]
## 시군구 계약일자 평당가격_일별
## <chr> <chr> <dbl>
## 1 서울특별시 강남구 논현동 20170701 3292.
## 2 서울특별시 강남구 논현동 20170703 3213.
## 3 서울특별시 강남구 논현동 20170704 1994.
## 4 서울특별시 강남구 논현동 20170707 2853.
## 5 서울특별시 강남구 논현동 20170710 3720.
## 6 서울특별시 강남구 논현동 20170711 2689.
## 7 서울특별시 강남구 논현동 20170712 4430.
## 8 서울특별시 강남구 논현동 20170713 5819.
## 9 서울특별시 강남구 논현동 20170714 3510.
## 10 서울특별시 강남구 논현동 20170715 3479.
## # … with 475 more rows그 동안 데이터의 유형을 배울 때, 그냥 아무런 의미 없이 추상적으로만 하였는데, 지금은 예제로 할 것이다.
필자가 근무하고 있는 부서는 보전(maintenance Dept') 부서이다. 생존분석이 매우 중요하다. 고장확률은 아래와 같다.
고장확률 = 1 - 생존확를
센서에서 고장에 대한 Event 알람 신호를 받아, 알람 신호기간 동안 고장 누적 분포 함수를 구할 수 있다.
아래 데이터는 Worcester Heart Attack Study 건으로 심장마비 후 생존 시간에 영향을 미칠 수 있는 연령, 성별 BMI 등 여러가지 요인 조사하였다. 모든 참가자는 심장마비 후 병원 입원시 시작되어 후속 조치에 따른 손실 또는 사망으로 끝난다.
아래 사용된 변수는 아래와 같다.
패키지를 불러 올 때 tidyvers 만 불러오면 이 모든 것을 불러 올 수 있는데, tidyvers 설치 할때 같이 설치 되는 패키지라 별도로 하는 것이 좋다. 왜냐하면 패키지끼지 같은 함수 이름을 쓰는데, 그런 경우에는 별도의 명령어를 써야 하는데, 매우 불편 하다.
그냥 아는 것만 꺼내다 쓰는 것이 좋다.
국토 교통부에서 데이터 셋을 받으면, 그 데이터가 좋지 않다. 데이터의 손질이 필요 하다. 손질을 하기 위해 가장 먼저 알아야 할 것은 필드 이름을 알아야 한다. 즉 필드 이름이 매우 잘되어야 데이터를 Manipulation(데이터 조작)을 잘 할 수 있다.
아래 예제와 같이 아파트 데이터 셋을 받는다. 데이터 셋을 받고 나서 colnames 함수를 이용해서 필드(행)이름을 아래와 같이 본다.
# 라이브러리 블러오기
library(readxl)
library(httr)
# 파일을 불러오는 path를 url 변수로 정의
url = "https://drive.google.com/u/0/uc?id=1PT_0-xJvM9KGp5zzx1jjcPAaS25epx4l&export=download"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
apt <- read_excel(tf, 1L) # 1L은 sheet 1을 이야기 하는 것이다. # 필드 이름 보기
colnames(apt)## [1] "시군구" "번지" "본번" "부번" ## [5] "단지명" "전용면적(㎡)" "계약년월" "계약일" ## [9] "거래금액(만원)" "층" "건축년도" "도로명" ## [13] "해제사유발생일"
# dplyr 라이브러리 불러오기
library(dplyr)# 아파트 데이터 프레임 필드(열) 이름 바꾸기
apt1 <- apt %>%
rename(전용면적_m3 = '전용면적(㎡)') %>%
rename(거래금액_만원 = '거래금액(만원)') %>%
rename(해제사유일 = 해제사유발생일)
colnames(apt1) ## [1] "시군구" "번지" "본번" "부번"
## [5] "단지명" "전용면적_m3" "계약년월" "계약일"
## [9] "거래금액_만원" "층" "건축년도" "도로명"
## [13] "해제사유일"
가장 얻기는 쉽지만, 활용하기 어려운 데이터 셋 주식 데이터 이다. 데이터의 시계열의 흐름을 보고 미래 사항을 예측 하지만, 아무리 큰 시장에도 심리적인 사항이 있어 정확하게 예측 하기는 어렵다. 하지만, 포토 폴리오를 예측을 할때, 해당 종목이 계절성을 띄고 있는지 선형성을 띄고 있는지 분석 하여, 전략을 수립 할 수 있다.
## 불필요 객체 지우기
rm(list = ls())
# library loading
library(tidyquant)
library(tqk)
# 주식 가져오는 함수
data_get_fun <- function(days, stock_code){
# Data get from tqk
stock <- tqk_get(stock_code, from= days)
stock$종목코드 <- stock_code
# R은 데이터를 가져와 전저리 할때 단순함
stock <- stock %>%
arrange(date) %>%
mutate(전일비 = c(NA, diff(adjusted))) %>%
filter(!is.na(전일비)) %>%
mutate(전일시세 = adjusted - 전일비) %>%
mutate(등락율 = round((adjusted - 전일시세)/전일시세 * 100, 2))
return(stock)
}
# 133690(나스닥), 319640 (골드선물)
# 오늘 보다 2019년 4월 부터 데이터를 가져온다.
days = as.Date("2019-04-01")
nasdaq <- data_get_fun(days, "133690")
gold <- data_get_fun(days, "319640")
# 데이터 프레임으로 취합 한다.
df <- rbind(nasdaq, gold)
# 월 단위로 묶는다.
df_sum <- df %>%
mutate(월 = paste0(substr(date, 1,4), substr(date, 6,7))) %>% group_by(종목코드, 월) %>%
summarise(가격 = mean(close))# 미국 나스닥
x <-ts( (df_sum %>% filter(종목코드 == "133690"))$가격 ,
start = c(2019, 04),
end = c(substr(today(),1,4),substr(today(),6,7)),
frequency= 12)
# 더블지수 평활
fd <- holt(x, h=5)
fd1 <- holt(x, damped = TRUE, phi = 0.9, h=5)
# 더블지수 평활
autoplot(x) +
autolayer(fd, series="홀트 기법", PI=FALSE) +
autolayer(fd1, series="감쇠 홀트 기법", PI=FALSE) +
ggtitle("더블지수 평활로 얻은 예측값") + xlab("연도") +
ylab("미국나스닥 가격") +
guides(colour=guide_legend(title="예측값"))
더블지수 평활
# 골드 선물
x <-ts( (df_sum %>% filter(종목코드 == "319640"))$가격 ,
start = c(2019, 04),
end = c(substr(today(),1,4),substr(today(),6,7)),
frequency= 12)
# 홀트 -윈터스의 덧셈
fit1 <- hw(x, seasonal = "additive" , h =5)
fit2 <- hw(x,seasonal="multiplicative", h =5)
autoplot(x) +
autolayer(fit1, series="HW 덧셈 예측", PI=FALSE) +
autolayer(fit2, series="HW 곱셈 예측",
PI=FALSE) +
ggtitle("홀트 윈터스로 얻은 예측값") + xlab("연도") +
ylab("골드 선물 가격") +
guides(colour=guide_legend(title="예측"))
홀트 윈터스 평활
통계 패키지 R에는 우리가 필요하다고 판단 되는 것들은 모두 있는 것 같다. 이번에는 주식 패키지 설치와 관련된 것인데, 주식패키지는 야후 웹페이지에서 크롤링 해서 가져오는 것이라, 속도가 매우 늦다. 차라리 직접 크롤링 하는것이 속도가 빠르다. 초 단타로 주식을 거래 하지 않는 이상은 실시간 데이터를 가져올 이유가 없다. 10분 단위도 문제 없기 때문에, 직접 크롤링 하는 것 보다. 아래 패키지를 사용하는 것이 오히려 편리 하다. 직접 크롤링 하는 것은 나중에 curl 하는 페이지에서 직접 다룰 예정이다.
앞에서는 html의 간단한 sheet를 소개 하였습니다. html은 주로 골격을 나타나는 것이라, 디자인을 하는데는 css로 하여야 합니다. 아래 코드와 같이 css 관련 하여 매우 간단하게 코딩 하겠습니다. body 부분의 css 코딩 ...
