rm(list=ls())
Sys.setenv(LANG = "en")
if(!require("tidyverse")) install.packages("tidyverse");library(tidyverse)
if(!require("lubridate")) install.packages("lubridate");library(lubridate)
if(!require("data.table")) install.packages("data.table");library(data.table)
if(!require("htmlTable")) install.packages("htmlTable");library(htmlTable)
if(!require("dlnm")) install.packages("dlnm");library(dlnm)
if(!require("gam")) install.packages("gam");library(gam)18 시간관련 자료 다루기
시계열 자료를 분석을 위해 필요한 부분은 날짜에 대한 데이터 변환입니다. 보통의 데이터가 오류를 없애기 위해 날짜를 문자나 숫자로 입력합니다. 이를 날짜 변수로 변환하는 과정을 실습하도록 하겠습니다. 엑셀과 SAS, R은 숫자로 변환된 날짜를 다시 날짜로 변환할 때 기준점이 다릅니다. 그렇게 때문에 주의가 필요합니다.
필요한 라이브러리를 불러 오겠습니다.
본 자료는 실제 자료를 기반으로 가상의 랜덤 값을 준 것 입니다. 가상의 자료이니 본 분석 결과를 사실로 받아 들이시지는 마시기 바랍니다.
hm0 <- readRDS("data/hm0.rds")
names(hm0) [1] "date" "temp_mean" "temp_max" "cloud_mean" "sun_day"
[6] "wind_mean" "rhum_mean" "temp" "tempf" "humidity"
[11] "heat_index" "event" "period" "region_id" "inclusion"
[16] "year" "circulation" "cvd" "endo" "gi"
[21] "heat" "infection" "injury" "mental" "nmd"
[26] "respiratory" "skin" "uro" "tmr" "occp1"
[31] "occp2" "occp3" "occp4" "occp5" "occp6"
[36] "occp7" "occp8" "occp9" "occp13" "occp99"
[41] "inside" "outside" "tomr" 18.1 Date create
오늘 날짜와 지금 시간, 초를 생성해 보겠습니다. 어떻게 나오나요?
today()
now()
Sys.time()
ymd, mdy, dmy를 사용하겠습니다.
tibble(ymd("2020-11-10"),
ymd("20201110"),
ymd("2020-Nov-10"),
mdy("November 10th, 2020"),
dmy("10-Nov-2020")) %>% t() [,1]
ymd("2020-11-10") "2020-11-10"
ymd("20201110") "2020-11-10"
ymd("2020-Nov-10") "2020-11-10"
mdy("November 10th, 2020") "2020-11-10"
dmy("10-Nov-2020") "2020-11-10"모든 날짜는 숫자로 저장되고 있습니다. 2020년 11월 10일은 숫자로는 18576입니다. 기준 날짜인 1970년 1월 1일부터 몇일이 지났는지를 적어 놓은 것입니다. 그래서 아래의 함수를 통해 숫자를 날짜로 변환 할 수 있습니다.
ymd("2020-11-10")%>% as.numeric()[1] 18576as.Date(18576, origin = "1970-01-01", tz = "UTC")[1] "2020-11-10"2020년 11월 10일부터 100일이 지난 날의 날짜는 무엇일까요? 아래와 같이 해볼 수 있습니다.
n100days <- ymd("2020-11-10")%>% as.numeric() +100
as.Date(n100days, origin = "1970-01-01")[1] "2021-02-18"여기에 hour, minute, second를 추가할 수 있습니다 .
tibble(year = 2020,
month = 11,
day = 10,
hour = 11,
minute= 15,
sec = 30) %>%
mutate(now = make_datetime(
year, month, day, hour, minute, sec
)) # A tibble: 1 × 7
year month day hour minute sec now
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dttm>
1 2020 11 10 11 15 30 2020-11-10 11:15:30그럼 그냥 문자로 사용하면 되지 왜 date 형식으로 변환을 하는 걸까요? Date-time components 가 필요하기 때문이기도 합니다. 그리고 연산이 가능해 집니다.
datetime <- ymd_hms("2020-11-10 11:30:55")
year(datetime)[1] 2020month(datetime)[1] 11mday(datetime) # 몇 일인지[1] 10yday(datetime) # 1월 1일 부터 몇일이 지난는지[1] 315wday(datetime) # 월요일 1, 일요일 7[1] 3몇월인지 무슨 요일인지 생성해 보겠습니다.
lubridate::month(datetime, label = TRUE)[1] 11월
12 Levels: 1월 < 2월 < 3월 < 4월 < 5월 < 6월 < 7월 < 8월 < ... < 12월lubridate::wday(datetime, label = TRUE)[1] 화
Levels: 일 < 월 < 화 < 수 < 목 < 금 < 토몇가지 연산을 위해 필요한 days, hours, minutes, seconds, weeks, months, years가 있습니다. 상식적 수순에서 이해할 수 있는데요,
- 시간간격
- days: 일수를 나타냅니다. days(100)는 100일의 시간 간격을 나타냅니다.
- hours: 시간을 나타냅니다. hours(24)는 24시간의 시간 간격을 나타냅니다.
- minutes: 분을 나타냅니다. minutes(60)은 60분의 시간 간격을 나타냅니다.
- seconds: 초를 나타냅니다. seconds(60)은 60초의 시간 간격을 나타냅니다.
- weeks: 주를 나타냅니다. weeks(2)는 2주의 시간 간격을 나타냅니다.
- months: 월을 나타냅니다. 주어진 월 수만큼의 평균 일 수를 계산합니다. months(1)은 한 달의 시간 간격을 나타냅니다.
- years: 연도를 나타냅니다. years(1)은 1년의 시간 간격을 나타냅니다.
우리는 날짜 데이터를 생성하고, 다양한 형식으로 날짜를 변환하는 방법을 배웠습니다. 이번에는 특정 날짜로부터 10일 전의 날짜를 계산하는 방법에 대해 알아보겠습니다. R에서 날짜 계산을 용이하게 해주는 lubridate 패키지를 사용하여 이러한 계산을 수행할 수 있습니다. 예를 들어, ’2023년 11월 28일’로부터 10일 전의 날짜를 찾고 싶다면, 다음과 같이 할 수 있습니다:
18.2 Tutor: heat wave and death
hm0 자료를 가지고 몇가지 실습을 해보도록 하겠습니다.
head(hm0) %>%
data.table() date temp_mean temp_max cloud_mean sun_day wind_mean rhum_mean temp
1: 2002-01-01 -4.4222222 4.6 1.6333333 7.577778 2.966667 55.62222 -1.5
2: 2002-01-02 -11.6222222 -4.7 0.1777778 8.033333 3.566667 40.97778 -8.7
3: 2002-01-03 -8.2222222 1.9 1.3111111 8.011111 3.266667 45.83333 -3.6
4: 2002-01-04 0.3777778 7.6 5.9666667 2.811111 3.800000 66.02222 2.7
5: 2002-01-05 -2.6777778 4.0 0.5333333 7.855556 2.833333 45.53333 0.5
6: 2002-01-06 -3.6222222 8.4 1.2333333 8.055556 1.477778 54.55556 1.2
tempf humidity heat_index event period region_id inclusion year circulation
1: 29.3 44.5 -1.51 0 0 42 1 2002 0
2: 16.4 27.1 -8.69 0 0 42 1 2002 10
3: 25.6 24.1 -3.58 0 0 42 1 2002 10
4: 36.8 47.9 2.68 0 0 42 1 2002 16
5: 32.8 29.4 0.45 0 0 42 1 2002 8
6: 34.2 32.4 1.23 0 0 42 1 2002 0
cvd endo gi heat infection injury mental nmd respiratory skin uro tmr occp1
1: 0 0 0 0 0 0 0 0 0 0 0 27 0
2: 0 1 23 0 5 6 6 26 76 10 17 27 0
3: 1 7 10 0 8 8 2 12 63 10 12 25 0
4: 1 7 12 0 5 4 2 17 57 16 12 29 0
5: 0 1 13 0 1 7 1 16 36 8 5 25 0
6: 0 0 0 0 0 0 0 0 0 0 0 23 1
occp2 occp3 occp4 occp5 occp6 occp7 occp8 occp9 occp13 occp99 inside outside
1: 0 0 0 2 3 0 0 0 22 0 27 0
2: 0 1 0 3 3 1 0 1 18 0 27 0
3: 0 0 1 0 6 0 0 0 18 0 25 0
4: 0 0 1 4 3 2 0 1 18 0 29 0
5: 1 0 1 1 4 0 0 0 18 0 25 0
6: 0 0 0 0 5 1 0 0 16 0 23 0
tomr
1: 5
2: 9
3: 7
4: 11
5: 7
6: 718.2.1 date create
날짜를 생성해 보겠습니다.
class(hm0$date)[1] "Date"class(hm0$year)[1] "numeric"2002년부터 2015년까지의 자료 입니다.
s1 <- hm0 %>%
group_by(year) %>%
count()
s1# A tibble: 14 × 2
# Groups: year [14]
year n
<dbl> <int>
1 2002 24820
2 2003 24820
3 2004 24892
4 2005 24824
5 2006 24824
6 2007 24835
7 2008 24894
8 2009 24822
9 2010 24856
10 2011 25200
11 2012 25267
12 2013 25220
13 2014 25213
14 2015 25212월, 일, 숫자형 월을 생성해 보겠습니다.
s2 <- hm0 %>% mutate(months = lubridate::month(date, label = TRUE, abbr = TRUE)) %>%
mutate(weekday = lubridate::wday(date, label = TRUE, abbr = TRUE)) %>%
mutate(month = substr(date, 6, 7)) %>%
filter(year(date) <= 2015) 월별 햇볕이 있었던 날의 수를 표로 나타내 보겠습니다.
s2 %>% group_by(months) %>%
summarize(`number of sunny days` = sum(sun_day, na.rm=TRUE)) %>%
htmlTable()| months | number of sunny days | |
|---|---|---|
| 1 | 1월 | 163634.1 |
| 2 | 2월 | 160058 |
| 3 | 3월 | 196690 |
| 4 | 4월 | 196703.6 |
| 5 | 5월 | 216350.7 |
| 6 | 6월 | 173973.8 |
| 7 | 7월 | 132143 |
| 8 | 8월 | 154961.4 |
| 9 | 9월 | 160004.7 |
| 10 | 10월 | 197543.7 |
| 11 | 11월 | 150447.1 |
| 12 | 12월 | 153073.2 |
월별 분포를 그림으로 그려보겠습니다.
library(ggplot2)
s2 %>% filter(year ==2015) %>%
group_by(months) %>%
ggplot(aes(x= months, y = sun_day, color = months)) +
geom_point()+
ggtitle('sunny day counts by months') +
ylab ('sunny day counts') +
guides(color = FALSE) 
이번에는 최고 기온가 일병 사망 숫자와의 관계를 살폐보겠습니다. 고온과 사망이므로 filter(month %in% c("06", "07", "08", "09"))을 이용하겠습니다. 월을 만들어 놓았으니 filter를 이용해서 필요한 월만 가져오면 되겠습니다. tomr 은 total occupation mortality 입니다. count 입니다. 35도가 넘으면 사망 숫자가 상승하고 있습니다.
s2 %>% filter(month %in% c("06", "07", "08", "09"),
temp_max >15) %>%
ggplot(aes(x = lag(temp_max), y = tomr)) +
geom_point(aes(color = months), alpha = 0.5) +
stat_smooth(method = "lm", formula = y ~ poly(x, 3), size = 1) +
theme_minimal() +
xlab('max temperature') + ylab('total death counts') +
labs(caption ="*row assocation, no lag with linear",
title=" Association between max temperature \n and daily total death counts "#, subtitle="Beta = -0.014, p = 0.158 before 2009\n Beta = 0.002, p = 0.011 from 2009"
)
18.3 home work
18.3.1 home0
2020-11-10로부터 100일전 날짜는?
2020-11-10로 부터 100일전 날짜의 요일은?
1978-01-21 에 태어난 사람이 2020-12-31까지 몇일을 살았을까?
1919-03-01 은 무슨 요일일까?
18.3.2 home work1
6, 7, 8, 9월에 고온과 농업인의 사망과의 관련성을 그려보세요. 농업인의 사망은 occp6 입니다.
18.3.3 home work2
12, 1, 2월에 기온과 농업인의 사망과의 관련성을 그려보세요. 기온은 temp_mean을 농업인의 사망은 occp6 입니다.
18.3.4 home work3
12, 1, 2월에 기온과 관리직의 사망과의 관련성을 그려보세요. 기온은 temp_mean을 관리직의 사망은 occp1 입니다.
18.4 Time series data analysis
기온등의 시계열적 변수가 질병에 어떻게 영향을 미치는지 알고자 할 때 어떻게 해야 할 까요? 우선 전체 직업군의 사망수 ’tomr’을 타겟 변수로 지정하겠습니다. 그리고 월별 사망률을 그려보겠습니다.
s2 <- s2 %>%
mutate(dz = tomr)
library(plotly)
s2 %>% group_by(months) %>%
plot_ly( x= ~ months, y = ~dz,
type = "box",
color = ~ months) %>%
layout(title = 'daily total mortality counts by months',
yaxis = list(title ='daily death counts'))월별 사망수의 차이가 관찰되고 있습니다. 이번에는 주별 차이를 관찰해 보겠습니다.
s2 %>% group_by(weekday) %>%
plot_ly( x= ~ weekday, y = ~dz,
type = "box",
color = ~ weekday) %>%
layout(title = 'daily total mortality counts by weekdays',
yaxis = list(title ='daily mortality counts'))고온에 대해서 알아보기 위해 6월부터 9월까지의 데이터로만 해보겠습니다. 몇다 방정식으로 모델을 구성했는지 관찰해 보세요
s2 %>% filter(month %in% c("06", "07", "08", "09"),
temp_max >15) %>%
ggplot(aes(x = lag(temp_max, 2), y = dz)) +
geom_point(aes(color = months), alpha = 0.5) +
stat_smooth(method = "lm", formula = y ~ poly(x, 3), size = 1) +
theme_minimal() +
xlab('max temperature') + ylab('total death counts') +
labs(caption ="*row assocation, no lab with linear",
title=" Association between max temperature \n and daily total death counts ")
기본적인 분석을 시행해 보겠습니다. 단순 회귀분석을 시행한다고 생각하고 진행해 보겠습니다. glm을 이용해 모델을 만들고, summary()를 통해 통계 값을 구해보겠습니다.
s2 <- s2%>% mutate(index = row_number())
mod1 <- glm(data = s2,
dz ~ temp_max + weekday + months + ns(index) +
region_id, family = quasipoisson())
summary(mod1)$coefficient[1,] Estimate Std. Error t value Pr(>|t|)
2.494459e+00 4.660113e-03 5.352786e+02 0.000000e+00 위에서 무엇이 단순 회귀 분서과 시계열 분석의 차이가 될까요? ns(index) 부분이 될 것입니다. 좀더 자세히 보겠습니다. 아래는 최고 온도의 변호와 3일의 간격을 두고, 주말 효과는 자유도 6으로 각 주가 각기 다른 효과를 준다고 가정하고, 월은 월을 dummy 변수로 취급하고, index에는 10을 주어서 시계열적 변화를 보정한다 이런 뜻 정도로 이해할 수 있습니다. 여기서 각 3, 6, 10 등의 숫자를 무엇을 넣어야 할지가 시계열 분석에서 고민해야 할 지점 중 하나입니다. 이는 다음 시간에 simulation하는 과정을 통해 알아보고 우선 진행해 보겠습니다.
g.mod1 <- glm(data = s2,
dz ~ lag(temp_max, 3) +
ns(week(date), 6) +
as.factor(months) +
ns(index, 10) + region_id, family = quasipoisson())
summary(g.mod1)$coefficient Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.689605306 2.030906e-02 132.4337587 0.000000e+00
lag(temp_max, 3) 0.008380698 2.140835e-04 39.1468619 0.000000e+00
ns(week(date), 6)1 -0.119623009 2.338702e-02 -5.1149320 3.140191e-07
ns(week(date), 6)2 -0.253319594 2.683492e-02 -9.4399234 3.752232e-21
ns(week(date), 6)3 -0.271267053 3.052486e-02 -8.8867588 6.320564e-19
ns(week(date), 6)4 -0.108926422 3.202452e-02 -3.4013445 6.706277e-04
ns(week(date), 6)5 -0.085539954 3.986917e-02 -2.1455163 3.191228e-02
ns(week(date), 6)6 -0.011664855 3.283241e-02 -0.3552848 7.223765e-01
as.factor(months).L 0.035538490 3.812434e-02 0.9321732 3.512477e-01
as.factor(months).Q 0.011562176 1.835534e-02 0.6299081 5.287552e-01
as.factor(months).C -0.005909232 1.247490e-02 -0.4736896 6.357216e-01
as.factor(months)^4 -0.028358118 9.514330e-03 -2.9805690 2.877332e-03
as.factor(months)^5 -0.007483562 8.059884e-03 -0.9284950 3.531515e-01
as.factor(months)^6 -0.001183933 6.422780e-03 -0.1843334 8.537520e-01
as.factor(months)^7 0.003316206 3.614635e-03 0.9174389 3.589134e-01
as.factor(months)^8 0.012798534 4.307070e-03 2.9715178 2.963519e-03
as.factor(months)^9 -0.003655662 2.963199e-03 -1.2336878 2.173201e-01
as.factor(months)^10 -0.001302690 3.038167e-03 -0.4287750 6.680873e-01
as.factor(months)^11 0.003369847 2.982248e-03 1.1299686 2.584902e-01
ns(index, 10)1 -0.058476904 6.546277e-03 -8.9328495 4.171157e-19
ns(index, 10)2 -0.057460708 8.608243e-03 -6.6750798 2.474642e-11
ns(index, 10)3 -0.220889502 7.787630e-03 -28.3641496 8.900638e-177
ns(index, 10)4 -0.093998725 8.275388e-03 -11.3588293 6.785182e-30
ns(index, 10)5 -0.114886433 7.986495e-03 -14.3850886 6.621038e-47
ns(index, 10)6 -0.140188072 8.105597e-03 -17.2952194 5.450567e-67
ns(index, 10)7 -0.013922461 7.906394e-03 -1.7609118 7.825421e-02
ns(index, 10)8 -0.167982193 6.611226e-03 -25.4086293 2.729150e-142
ns(index, 10)9 0.154066429 1.314857e-02 11.7173555 1.052933e-31
ns(index, 10)10 -0.038049443 5.825723e-03 -6.5312827 6.529763e-11
region_id -0.000312727 1.516356e-06 -206.2358381 0.000000e+0018.4.1 outside 사망에 대해서
몇가지를 탐색해보니, 여름철 야외 사망에 대해 연구해 보면 좋겠다는 생각이 드네요.
s3 <- s2 %>%
# filter(region_id %in%c(46)) %>%
filter(month %in% c( '07', '08', '09')) %>%
mutate(index = row_number()) 위에서 얼마의 lag time을 주어야 할지 시물레이션을 돌리는 방법을 생각했었습니다. 한 방법은 가능한한 모든 것을 다 해보고 가장 적당하다고 생각되는 것을 해보는 것입니다.
library(dlnm)
cb1<- crossbasis(s3$temp_max, lag=14, argvar=list("lin"), arglag = list(fun="poly", degree=3))
model1<-glm(data=s3,
outside ~ cb1 + weekday + months + ns(index,10) + as.factor(region_id)#+
#cloud_mean + sun_day #+ wind_mean + rhum_mean
, family=quasipoisson())
summary(model1)
Call:
glm(formula = outside ~ cb1 + weekday + months + ns(index, 10) +
as.factor(region_id), family = quasipoisson(), data = s3)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.820e+01 8.141e-01 -71.488 < 2e-16 ***
cb1v1.l1 1.096e-02 5.851e-04 18.726 < 2e-16 ***
cb1v1.l2 -7.797e-02 5.881e-03 -13.259 < 2e-16 ***
cb1v1.l3 1.470e-01 1.408e-02 10.441 < 2e-16 ***
cb1v1.l4 -8.375e-02 9.276e-03 -9.029 < 2e-16 ***
weekday.L -5.868e-03 6.136e-03 -0.956 0.338954
weekday.Q 2.371e-02 6.155e-03 3.853 0.000117 ***
weekday.C 2.768e-02 6.144e-03 4.506 6.61e-06 ***
weekday^4 -8.264e-03 6.164e-03 -1.341 0.180078
weekday^5 1.279e-03 6.155e-03 0.208 0.835338
weekday^6 1.140e-02 6.188e-03 1.842 0.065418 .
months.L 2.777e-02 4.420e-03 6.283 3.34e-10 ***
months.Q -2.379e-02 5.404e-03 -4.403 1.07e-05 ***
ns(index, 10)1 6.024e+01 8.031e-01 75.009 < 2e-16 ***
ns(index, 10)2 6.128e+01 8.204e-01 74.693 < 2e-16 ***
ns(index, 10)3 6.041e+01 8.124e-01 74.368 < 2e-16 ***
ns(index, 10)4 6.133e+01 8.170e-01 75.063 < 2e-16 ***
ns(index, 10)5 5.921e+01 8.139e-01 72.748 < 2e-16 ***
ns(index, 10)6 5.872e+01 8.165e-01 71.918 < 2e-16 ***
ns(index, 10)7 5.921e+01 8.152e-01 72.634 < 2e-16 ***
ns(index, 10)8 4.090e+01 5.737e-01 71.284 < 2e-16 ***
ns(index, 10)9 1.123e+02 1.543e+00 72.776 < 2e-16 ***
ns(index, 10)10 2.327e+01 3.315e-01 70.180 < 2e-16 ***
as.factor(region_id)26 -1.671e+00 2.638e-02 -63.357 < 2e-16 ***
as.factor(region_id)27 -1.448e+00 2.412e-02 -60.036 < 2e-16 ***
as.factor(region_id)28 -1.501e+00 1.914e-02 -78.429 < 2e-16 ***
as.factor(region_id)29 -2.247e+00 3.379e-02 -66.487 < 2e-16 ***
as.factor(region_id)30 -1.984e+00 3.032e-02 -65.452 < 2e-16 ***
as.factor(region_id)31 -2.744e+00 4.254e-02 -64.510 < 2e-16 ***
as.factor(region_id)41 3.948e-02 1.198e-02 3.297 0.000978 ***
as.factor(region_id)42 -1.624e+00 1.352e-02 -120.127 < 2e-16 ***
as.factor(region_id)43 -1.776e+00 1.674e-02 -106.135 < 2e-16 ***
as.factor(region_id)44 -1.388e+00 1.448e-02 -95.896 < 2e-16 ***
as.factor(region_id)45 -1.726e+00 1.455e-02 -118.582 < 2e-16 ***
as.factor(region_id)46 -1.445e+00 1.366e-02 -105.782 < 2e-16 ***
as.factor(region_id)47 -8.196e-01 1.206e-02 -67.970 < 2e-16 ***
as.factor(region_id)48 -1.236e+00 1.301e-02 -94.956 < 2e-16 ***
as.factor(region_id)4950 -2.781e+00 2.364e-02 -117.630 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 1.109338)
Null deviance: 363575 on 88109 degrees of freedom
Residual deviance: 97171 on 88072 degrees of freedom
(결측으로 인하여 14개의 관측치가 삭제되었습니다.)
AIC: NA
Number of Fisher Scoring iterations: 9pred1.cb1 <-crosspred(cb1, model1, at=1:100, bylag=0.1, cumul=FALSE)
plot(pred1.cb1, "slices", var=1, col=3, ylab="Relative risk of total death count", #ci.arg=list(density=50, lwd=1),#
main="Temporal effect of max temperature on total death count",
ci.arg=list(density=100, lwd=3),
cex.main =1.5,
xlab="Lag (days)", # family="A",#ylim=c(0.980, 1.02),
col='black') ;grid()
title(main="* by 1 unit increments of max temperature",
family="A",
adj=1, line=0, font.main=1, cex=1)
dlnm distribute lag liner and non-linear model 의 약자로 자세한 것은 http://www.ag-myresearch.com/package-dlnm.html 을 찾아 보시기 바랍니다. 여기서 주목할 점은 14일까지의 lag time을 모두 그려서 lag time이 0 인 날 즉 최고 온도가 높은 날 바로 당이에 사망이 높게 올라간다는 것입니다. 이후 4일 후에는 havesting effect가 보이는 군요. 따라서 lag time은 0 설정하고 위험도를 산출하는 방법을 사용하면 좋을 것 같습니다. 또는 그림 그대로를 받아 들여도 좋습니다.
18.5 Time series data analysis, regression
18.6 introduction
시계열적 자료를 분석하고 나타나는 현상과 특정 요인과 관련성을 탐색해보는 시간입니다. 예를 들어 미세먼지가 높은 날 심혈관 질환이 발생하는가?에 대한 질문에 답하기 위해서 생가할 것이 몇가지 있습니다. 미세먼지가 높은 날이란? 심혈관 질환 사망이 높은 날이란? 이 두가지 요소를 검토하게 됩니다.
그런데 심혈관 질환의 사망은 요일마다 다르고, 계절에 따라 변동하며, 장기 적으로는 점차 증가 또는 감소를 합니다. 그런데 미세먼지도 점차 증가하고 있으니, 단순 상관관계를 보면 미세먼지도 증가 심혈관 사망도 증가하면 양의 관련성을 보이게 됩니다. GDP와 자살의 관계를 보면 어떨까요? 우리나라의 자살률은 증가하고 있습니다. 그런데 GDP도 증가하고 있습니다. 그러니 GDP의 증가와 자살의 증가는 양의 상관관계가 있다고 나옵니다. 맞나요? 네 심혈관 사망, 자살의 증가의 계절적 요소, 장기간 추세(trend)가 아니라 변동이 미세먼지나 GDP의 변동가 어떠한 관계가 있는지가 우리의 궁금증일 것 입니다. 이러한 궁금증을 R을 이용해서 풀어보도록 하겠습니다.
18.7 미세먼지와 심혈관사망
우선 몇가지 시계열 자료 분석의 이해를 돕기 위해 시뮬레이션 자료를 이용해 보겠습니다.
x를 일(day) 로 생각하고, 300일 동안 랜덤 변수 y1과 이에 4.5를 곱한 pm(미세먼지)를 가상으로 만들어 보겠습니다.
set.seed(1)
x <- 1:300
y1 <- 5*rnorm(300, sd=.1)+15
pm <- y1*4.5
plot(x, pm, type='l')
여기에
sin()함수를 통해 계절적 요소를 넣고, 0.03을 곱해 long term trend 가 서서히 증가하는 것으로 가정했습니다.
y2 <- y1*5+ sin(x/2)*5+ x * 0.03
y2[y2< 0]<-0
y3<-round(y2)
plot(y3, type='l')
지연 효과와 특정 이벤트가 있는 날을 넣어 보았습니다. 그리고 dataframe을 만들었습니다.
lag <-6
mean(y3)[1] 79.58667death <- c(rep(c(80,79,81), (lag/3)), y3[1:(length(y3)-lag)])
event <- c(rep(1, 30), rep(1, 30), rep(0, 240))
eventd <- c(rep(40,30), rep(30, 30), rep(0, 240))
death2<-death+eventd+10
gg <- data.frame(x, pm, y3, death, event, death2)
head(gg) x pm y3 death event death2
1 1 66.09048 76 80 1 130
2 2 67.91320 80 79 1 129
3 3 65.61984 78 81 1 131
4 4 71.08938 84 80 1 130
5 5 68.24139 79 79 1 129
6 6 65.65395 74 81 1 131이제 그림을 그려 보겠습니다. 첫 50일에 이벤트가 있어 심혈관 사망이 높고 이후 계절적 요소를 보이며 서서히 증가하고 있습니다. 미세먼지는 random + 계절적 요소로 만들었고요.
plot(x, pm, type="l", col=grey(0.5), ylim=c(50, 140), xlim=c(0, 300))
grid()
lines(x, death2, col=grey(0.7), type="p", cex=0.5)
이제 단순 회귀 분석을 해보겠습니다. 어떠한 관계가 관찰되시나요. event 때 많이 사망하고, 미세먼지와는 관련이 없네요. 분명 미세먼지와 관련있게 시뮬레이션 해서 만든 자료인데요. 맞습니다.
lag과seasonality보정이 않되었네요.
mt3 <- glm(death2 ~ x+sin(x/2)+pm+event)
summary(mt3)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.10455149 6.861474711 13.1319513 2.476263e-31
x 0.02379324 0.003500538 6.7970236 5.915161e-11
sin(x/2) -4.41585403 0.308633540 -14.3077581 1.247252e-35
pm -0.06144597 0.101078610 -0.6079028 5.437196e-01
event 35.05109683 0.757861036 46.2500315 3.230388e-137그림으로 확인해 보겠습니다. 무언가 잘못 예측이 되고 있죠?
plot(x, pm, type="l", col=grey(0.5), ylim=c(50, 140), xlim=c(0, 300))
grid()
lines(x, death2, col=grey(0.7), type="p", cex=0.5)
mp3 <- c( predict(mt3))
lines(x, mp3, col=75)
차라리 이렇게 해보는 것은 어떨까요? 시계열적인 요소를 뺀 상태 (residual) 과 미세먼지가 관련이 있나 보는 것입니다 .
mt2 <- glm(death2 ~ x+sin(x/2)+event)
resid_mt2 <-resid(mt2)
risk.m0<-glm(resid_mt2 ~ pm, family=gaussian)
summary(risk.m0)
Call:
glm(formula = resid_mt2 ~ pm, family = gaussian)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.14114 6.79037 0.61 0.542
pm -0.06128 0.10043 -0.61 0.542
(Dispersion parameter for gaussian family taken to be 14.18007)
Null deviance: 4230.9 on 299 degrees of freedom
Residual deviance: 4225.7 on 298 degrees of freedom
AIC: 1650.9
Number of Fisher Scoring iterations: 2risk.mp0 <- c( predict(risk.m0))
plot(pm, resid_mt2, type='p', cex=0.5)
lines(pm, (risk.mp0), col=25)
저는 이것이 더 직관적인데요. 심혈관사망에서 시계열적으로 변동이 있는 부분을 뺀 나머지 (residual) 이 pm 이 변동할 때 같이 변동하면 관련성이 있다고 보는 것이지요.
자 이제 lag 을 줘서 관찰해 보겠습니다. lag을 주면 초반 데이터 숫자가 맞지 않는데요. 이때 pm의 평균 값으로 결측치를 대신해서 해결해 보겠습니다.
mean(pm)[1] 67.57556lag.pm<-6
pm.lag <- c(rep(67.5, lag.pm), pm[1:(length(pm)-lag.pm)])
resid_mt3 <-resid(mt3)
risk.m1<-glm(resid_mt3 ~ pm.lag, family=gaussian)
summary(risk.m1)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) -76.437599 5.21757250 -14.65003 5.620794e-37
pm.lag 1.131554 0.07720006 14.65743 5.276143e-37risk.mp1 <- c( predict(risk.m1))
plot(pm.lag, resid_mt3, type='p', cex=0.5)
lines(pm.lag, risk.mp1, col=25)
네 이제 pm 과 양의 심혈관 사망에 양의 상관관계가 생겼네요. 우리가 원했던 데이터를 그렇게 만들었었습니다.
> 그림으로 관찰해 보면 빨간색이 lag time을 준것입니다. 누가 더 사망과 관련이 있어 보이나요?
plot(x, resid_mt3, type="l", col=grey(0.5), ylim=c(-15, 40), xlim=c(0, 300))
grid()
lines(x, (pm-50), col=grey(0.7), type="l", cex=0.5)
lines(x, (pm.lag-60), col='red', type="l", cex=0.5)
지금 까지 고려한 것은
sin()으로 계절적 요소,lag으로 지연 효과를 고려해서 시계열적 요소를 없앤다음 (residual), pm과 심혈관사망의 관계를 분석하는 방식으로 해보았습니다. 좀더 쉽게 이것을 해보겠습니다.
library(mgcv)
mgam<- gam(death2 ~ s(x, bs="cc", k=100)+event, family=gaussian)
p <- predict(mgam)
plot(x, pm, type="l", col=grey(0.5), ylim=c(40, 150), xlim=c(0, 300), cex=2)
grid()
lines(x, death2, col=grey(0.7), type="p", cex=0.5)
legend(x=250, y=70, 'PM10')
legend(x=150, y=65, 'pm10. lag')
legend(x=210, y=110, 'Obs_death')
legend(x=10, y=50, 'Residual(Obs_Death - Gam(fitting)')
lines(x, p)
lines(x, (resid(mgam)+50), col='blue')
lines(x, pm.lag-10, col='red')
이것을 회귀 분석으로 구해보겠습니다. k 가 높을 수록 모형은 어떠한 가요? 네 즉 위에 lag time 과 k 값을 어떻게 조정하는 지를 고려해야 합니다. 우선 lag time 어떻게 찾을 까요?
mgam<- gam(death2 ~ s(x, bs="cc", k=100)+event, family=gaussian)
p <- predict(mgam)
risk.pp1 <-glm(death2 ~ p+pm.lag,family=gaussian)
summary(risk.pp1)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) -58.3872771 1.937186312 -30.14025 2.402418e-92
p 1.0000436 0.004566766 218.98286 0.000000e+00
pm.lag 0.8642815 0.028266084 30.57663 9.482839e-94AIC(risk.pp1)[1] 885.3135mgam150<- gam(death2 ~ s(x, bs="cc", k=10)+event)
p150 <- predict(mgam150)
risk.pp150 <-glm(death2 ~ p150+ pm.lag, family=gaussian)
summary(risk.pp150)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) -72.5953393 6.74682944 -10.75992 5.109224e-23
p150 0.9979243 0.01630871 61.18966 2.087420e-170
pm.lag 1.0776412 0.09754736 11.04736 5.349868e-24AIC(risk.pp1, risk.pp150) df AIC
risk.pp1 4 885.3135
risk.pp150 4 1629.2747lag 을 찾아 보겠습니다. pm에 대해 lag을 10일까지, 3차 방정식 형태로, 그리고 이를 통해 회귀 분석을 시행해 보세요.
library(dlnm)
cb1.pm <-crossbasis(pm, lag=10, argvar=list(fun="lin"),
arglag=list(fun="poly", degree=3))
model1 <-glm(death2 ~ cb1.pm+x+event ,
family=gaussian )
pred1.pm <-crosspred(cb1.pm, model1, at=0:100, bylag=0.1, cumul=TRUE)
plot(pred1.pm, "slices", var=1, col=3, ylab="RR", ci.arg=list(density=15,lwd=2),
#cumul = TRUE,
main="Association with a 1-unit increase in PM10")
6일을 lag time으로 하면 좋겠네요.
이제 6일을 lag time으로 설정해서 회귀 분석을 수행하면 된다는 것을 알았습니다. 이제 남은 것은 시계열적 요소를 어떻게 찾고, 어떻게 보정할지, 그리고 이 과정을 어떻게 합리적으로 할지 더 논의 하면 됩니다.