DLNM을 이용한 대기오염원 건강영향평가

# DLNM을 이용한 대기오염원 건강영향평가
----
## **질병관리청 건강보험 청구자료 데이터 교육**
### 방태모
### 2022-07-11

---

# .brand-blue[목차]

- .font130[.content-box-blue[모형 소개]]

<br>

- .font130[.content-box-blue[실무 적용 예시]]

<br>

- .font130[.content-box-blue[실습]]

]

---

# 모형 소개

---

## .brand-blue[모형 소개]

`\begin{align}
g(\mu_t) = \alpha + s({\boldsymbol{A}}_{t}, {\rm{df}}_p, {\rm{lag}}, {\rm{df}}_l ; {\boldsymbol{\beta}}) + {\boldsymbol{\gamma}}{\boldsymbol{M}}_{t}
\end{align}`

- 시계열 회귀모형

- .brand-blue[비선형적 효과]와 .brand-blue[지연 효과(delayed effect)] 반영 가능

- 생태학적 연구(ecological study)

- .brand-blue[환경적 요인(stressor)]으로 Outcome( `$Y_t$` ) 모델링
  - 특정 .brand-blue[예측변수] `$X_t$`.brand-blue[(exposure)]의 효과는 관측 시점 당시 뿐만이 아닌 .brand-blue[지연 효과] 빈번히 존재
]

.pull-right[
<div class="figure">
<img src="./fig/3dplots.png" alt="Fig 1. Gasparrini (2011)" width="400" height="350" />
<p class="caption">Fig 1. Gasparrini (2011)</p>
</div>
]

---

# 실무 적용 예시

---
## .brand-blue[실무 적용 예시]

`\begin{align}
{\rm{log}}(E(Y_{ijt})) = \alpha + {\rm{NS}}({\boldsymbol{A}}_{ijt}, {\rm{df}}_p, {\rm{lag}}, {\rm{df}}_l ; {\boldsymbol{\beta}}_{ij}) + {\boldsymbol{\gamma}}_{ij}{\boldsymbol{M}}_{ijt} + {\boldsymbol{\delta}}_{ij}{\boldsymbol{S}}_{t} 
\end{align}`

.pull-left[
 - 우리나라 6개 대도시에 대해 6종의 대기오염원이 정신과 질환에 미치는 건강영향평가 수행
  
  - 대기오염원(exposure): `${\rm{PM}}_{10}$`, `${\rm{PM}}_{2,5}$`, `${\rm{NO}}_{2}$`, `${\rm{SO}}_{2}$`, `${\rm{O}}_{3}$`, `${\rm{NO}}_{2}$`
  - 기상요인(covariates): 강수량, 상대습도, 풍속, 기온, 일교차
  - 🔗 [데이터 제공](https://github.com/G-ABCwork/Data_Useful)
  
- .brand-blue[Single exposure DLNM] 고려

]

.pull-right[
 - 최적화 알고리즘 개발
  - Maximum lag days( `${\rm{lag}}$` ): `$[7, 8, \cdots, 31]$`
  - Degrees of freedom in predictor space( `${\rm{df}}_p$` ): `$[2, 3, 4, 5]$`
  - Degrees of freedom in additional lag dimension ( `${\rm{df}}_l$` ): `$[2, 3, 4, 5]$`
  
- 다변량 메타분석 수행

- 6개 대도시에 관한 분석 결과 통합
  
  - (Gasparrini, Armstrong, and Kenward 2012)
]

---

# 실습

---

## .brand-blue[패키지 설치 및 로딩]

```r
install.packages(c("tidyverse", "dlnm", 
                   "splines", "stringr", 
                   "lubridate"))
```

```r
library(tidyverse)
library(dlnm)
library(splines)
library(stringr)
library(lubridate)
ggplot2::theme_set(theme_bw())
```
]

- {dlnm} 패키지의 1987-2000년 시카고의 일별 사망건수, 기상요인, 대기오염원에 관한 정보가 담겨있는 `chicagoNMMAPS` 자료 이용

```r
chicago <- chicagoNMMAPS |> 
    as_tibble() |> 
    select(date, time, year, pm10, temp, death)
glimpse(chicago)
```

```
## Rows: 5,114
## Columns: 6
## $ date  <date> 1987-01-01, 1987-01-02, 1987-01-03, 1987-01-04, 1987-01-05, 198…
## $ time  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
## $ year  <dbl> 1987, 1987, 1987, 1987, 1987, 1987, 1987, 1987, 1987, 1987, 1987…
## $ pm10  <dbl> 26.95607, NA, 32.83869, 39.95607, NA, 40.95607, 33.95607, 28.956…
## $ temp  <dbl> -0.2777778, 0.5555556, 0.5555556, -1.6666667, 0.0000000, 4.44444…
## $ death <int> 130, 150, 101, 135, 126, 130, 129, 109, 125, 153, 124, 111, 104,…
```

]

---

## .brand-blue[데이터 출력]

```r
head(chicago)
```

```
## # A tibble: 6 × 6
##   date        time  year  pm10   temp death
##   <date>     <int> <dbl> <dbl>  <dbl> <int>
## 1 1987-01-01     1  1987  27.0 -0.278   130
## 2 1987-01-02     2  1987  NA    0.556   150
## 3 1987-01-03     3  1987  32.8  0.556   101
## 4 1987-01-04     4  1987  40.0 -1.67    135
## 5 1987-01-05     5  1987  NA    0       126
## 6 1987-01-06     6  1987  41.0  4.44    130
```

]

## .brand-blue[EDA]

```r
chicago |> 
    ggplot(aes(x = date, y = pm10)) + 
    geom_line()
```

]

---

## .brand-blue[Cross-basis matrix 만들기]

```r
cb_pm <- crossbasis(chicago$pm10,
                    lag = 14,
                    argvar = list(fun = "ns", df = 3),
                    arglag = list(fun = "ns", df = 3))
summary(cb_pm)
```

```
## CROSSBASIS FUNCTIONS
## observations: 5114 
## range: -3.049835 to 356.1768 
## lag period: 0 14 
## total df:  9 
## 
## BASIS FOR VAR:
## fun: ns 
## knots: 24.09097 37.53545 
## intercept: FALSE 
## Boundary.knots: -3.049835 356.1768 
## 
## BASIS FOR LAG:
## fun: ns 
## knots: 7 
## intercept: TRUE 
## Boundary.knots: 0 14
```

]

---

## .brand-blue[모형 적합]

```r
num_year <- chicago |> 
    select(year) |> 
    unique() |> 
    nrow()
mod <- glm(death ~ cb_pm + temp + ns(time, 7*num_year), 
           family = quasipoisson(), data = chicago)
```

## .brand-blue[예측 수행]

```r
pred_pm <- crosspred(cb_pm, mod, 
                     cen = chicago$pm10 |> 
                         median(na.rm = T),
                     at = 10:80, by = 1)
```

---

## .brand-blue[분석 결과 시각화]

### .black[3D plot]

```r
plot(pred_pm, xlab = "PM10", zlab = "RR",
     theta = 210, phi = 30, lphi = 30, border = "gray40")
```

]

### .black[Contour plot]

```r
plot(pred_pm, "contour", xlab = "PM10", ylab = "Lag",
     key.title = title("RR"))
```

]

---

## .brand-blue[분석 결과 시각화]

### .black[Overall cumulative association plot]

```r
plot(pred_pm, "overall", col = "tomato", 
     lwd = 2, xlab = "PM10", ylab = "RR")
rug(chicago$pm10, quiet = TRUE)
```
]

.pull-right[
<img src="kdca_dlnm_files/figure-html/OCA-show-1.png" width="2100" style="display: block; margin: auto;" />
]

---

## .brand-blue[분석 결과 시각화]

### .black[High PM10 effect]

```r
plot(pred_pm, "slices",
     var = chicago$pm10 |> 
         quantile(0.9, na.rm = TRUE) |> 
         round(0),
     col = "tomato", lwd = 2, ylab = "RR", 
     main = "High PM10 effect (vs 90th quantile)")
```
]

.pull-right[
<img src="kdca_dlnm_files/figure-html/high-show-1.png" width="2100" style="display: block; margin: auto;" />
]

---

## .brand-blue[분석 결과 시각화]

### .black[Low PM10 effect]

```r
plot(pred_pm, "slices",
     var = chicago$pm10 |> 
         quantile(0.1, na.rm = TRUE) |> 
         round(0),
     col = "tomato", lwd = 2, ylab = "RR", 
     main = "Low PM10 effect (vs 10th quantile)")
```

]

.pull-right[
<img src="kdca_dlnm_files/figure-html/low-show-1.png" width="2100" style="display: block; margin: auto;" />

]

---
class: inverse, middle, center

# References

---

## .brand-blue[References]

[1] Gasparrini, Antonio, Benedict Armstrong, and M.G. Kenward. “Distributed Lag Non-Linear Models.” Statistics in Medicine 29 (September 20, 2010): 2224–34. https://doi.org/10.1002/sim.3940.

[2] Gasparrini, Antonio. “Distributed Lag Linear and Non-Linear Models in R: The Package Dlnm.” Journal of Statistical Software 43 (July 1, 2011): 1–20. https://doi.org/10.18637/jss.v043.i08.

---

# Thanks!

<a href="mailto:favorite@kakao.com">
.white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M440 6.5L24 246.4c-34.4 19.9-31.1 70.8 5.7 85.9L144 379.6V464c0 46.4 59.2 65.5 86.6 28.6l43.8-59.1 111.9 46.2c5.9 2.4 12.1 3.6 18.3 3.6 8.2 0 16.3-2.1 23.6-6.2 12.8-7.2 21.6-20 23.9-34.5l59.4-387.2c6.1-40.1-36.9-68.8-71.5-48.9zM192 464v-64.6l36.6 15.1L192 464zm212.6-28.7l-153.8-63.5L391 169.5c10.7-15.5-9.5-33.5-23.7-21.2L155.8 332.6 48 288 464 48l-59.4 387.3z"></path></svg> favorite@kakao.com]
</a>

<a href="https://github.com/be-favorite">
.white[<svg viewBox="0 0 496 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"></path></svg> @be-favorite]
</a>

<a href="https://twitter.com/TaemoBang">
.white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"></path></svg> @TaemoBang]
</a>

<a href="https://github.com/be-favorite/Presentation_archive">
.white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M326.612 185.391c59.747 59.809 58.927 155.698.36 214.59-.11.12-.24.25-.36.37l-67.2 67.2c-59.27 59.27-155.699 59.262-214.96 0-59.27-59.26-59.27-155.7 0-214.96l37.106-37.106c9.84-9.84 26.786-3.3 27.294 10.606.648 17.722 3.826 35.527 9.69 52.721 1.986 5.822.567 12.262-3.783 16.612l-13.087 13.087c-28.026 28.026-28.905 73.66-1.155 101.96 28.024 28.579 74.086 28.749 102.325.51l67.2-67.19c28.191-28.191 28.073-73.757 0-101.83-3.701-3.694-7.429-6.564-10.341-8.569a16.037 16.037 0 0 1-6.947-12.606c-.396-10.567 3.348-21.456 11.698-29.806l21.054-21.055c5.521-5.521 14.182-6.199 20.584-1.731a152.482 152.482 0 0 1 20.522 17.197zM467.547 44.449c-59.261-59.262-155.69-59.27-214.96 0l-67.2 67.2c-.12.12-.25.25-.36.37-58.566 58.892-59.387 154.781.36 214.59a152.454 152.454 0 0 0 20.521 17.196c6.402 4.468 15.064 3.789 20.584-1.731l21.054-21.055c8.35-8.35 12.094-19.239 11.698-29.806a16.037 16.037 0 0 0-6.947-12.606c-2.912-2.005-6.64-4.875-10.341-8.569-28.073-28.073-28.191-73.639 0-101.83l67.2-67.19c28.239-28.239 74.3-28.069 102.325.51 27.75 28.3 26.872 73.934-1.155 101.96l-13.087 13.087c-4.35 4.35-5.769 10.79-3.783 16.612 5.864 17.194 9.042 34.999 9.69 52.721.509 13.906 17.454 20.446 27.294 10.606l37.106-37.106c59.271-59.259 59.271-155.699.001-214.959z"></path></svg> Presentation archive]
</a>

]]