지연 효과와 비선형 관계를 모두 고려한 모델인 DLNM의 확장 버전, the Extended DLNM에 대해 소개하고, R의 dlnm 패키지를 이용하여 the Extended DLNM을 적합시키는 방법에 대해 소개합니다.

1. DLNM이란?

어떤 노출 요인(exposure, \(X\))은 관측 당시 뿐만 아니라, 일정 시간이 지난 후 그 효과가 발생할 수 있습니다. 즉, ‘지연 효과’를 가집니다. 모형에 지연 효과를 포함한다면 더 효과적인 예측을 진행할 수 있으며, 이를 위해 시차에 대한 차원이 추가적으로 필요합니다.

지연 효과를 고려한 모델로 DLM(Distributed Lag Models)가 존재합니다. 하지만 DLM은 선형적 효과를 다루는데 더 적합하므로, 비선형적 효과를 표현하는데 한계가 있습니다. DLNM은 DLM과 달리 비선형적 효과를 반영한, 더 flexible한 모형입니다.

DLNM은 cross-basis(교차 기저)를 정의하여 \(X\)(exposure)와 그것의 시차(lag)에 따른 \(Y\)(response)의 관계를 고려함으로써 지연 효과까지 반영합니다. 즉, cross-basis란 exposure-lag-response의 관계를 3차원 상의 공간에서 표현한 것입니다. cross-basis를 정의하기 위해 exposure-response 차원과 lag-response 차원에서 사용할 basis function(기저 함수)을 선택해야 하는데, 어떤 basis function을 선택하느냐에 따라 선형적/비선형적 관계를 고려할지 결정할 수 있습니다.

2. the Extended DLNM이란?

DLNM과 the Extended DLNM의 주요 차이점은 matrix of exposure histories의 정의입니다. matrix of exposure histories는 \(n\)개의 관측치 각각에 대한 lag \(\ell\)에서의 exposure 값들로 이루어진 행렬입니다. 즉, exposure 값들로 이루어진 \(n \times (L - \ell_0 + 1)\) 크기의 행렬을 말합니다. 이 matrix of exposure histories는 study design과 수집된 정보에 따라 다르게 구성됩니다.

DLNM을 사용하는 일반적인 시계열 데이터에서 matrix of exposure histories의 \([t, \ell]\) 값은 \([t+1, \ell+1]\) 값과 동일합니다. \(n\)개의 관측치가 하나의 변수에 대한 여러 시점의 관측값이기 때문입니다.

예를 들어, 다음과 같은 데이터를 살펴봅시다. 아래 시계열 데이터는 1987-01-01부터 1987-01-15까지의 미세먼지(pm10)를 기록한 것입니다.

This is dlnm 2.4.7. For details: help(dlnm) and vignette('dlnmOverview').

         date     pm10
1  1987-01-01 26.95607
2  1987-01-02       NA
3  1987-01-03 32.83869
4  1987-01-04 39.95607
5  1987-01-05       NA
6  1987-01-06 40.95607
7  1987-01-07 33.95607
8  1987-01-08 28.95607
9  1987-01-09 32.34877
10 1987-01-10       NA
11 1987-01-11 14.95607
12 1987-01-12 18.95607
13 1987-01-13 45.95607
14 1987-01-14 35.95607
15 1987-01-15 26.93720

이 데이터의 matrix of exposure histories는 아래와 같습니다. 이 경우, \([t, \ell]\) 값은 \([t+1, \ell+1]\) 값과 동일합니다.

               lag0     lag1     lag2     lag3     lag4     lag5
1987-01-10       NA 32.34877 28.95607 33.95607 40.95607       NA
1987-01-11 14.95607       NA 32.34877 28.95607 33.95607 40.95607
1987-01-12 18.95607 14.95607       NA 32.34877 28.95607 33.95607
1987-01-13 45.95607 18.95607 14.95607       NA 32.34877 28.95607
1987-01-14 35.95607 45.95607 18.95607 14.95607       NA 32.34877
1987-01-15 26.93720 35.95607 45.95607 18.95607 14.95607       NA

하지만 \(n\)개의 관측치가 서로 관련이 없는 경우, 이는 성립하지 않습니다. 이런 데이터는 the Extended DLNM을 사용해야 합니다.

예를 들어, 아래의 데이터는 각 환자들별로 매 시차의 약물 복용량을 기록한 것입니다. 이 경우, \([t, \ell]\) 값은 \([t+1, \ell+1]\) 값과 동일하지 않습니다.

  lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1   37   37   37   37   37   37   37   40   40   40    40    40    40    40
2    0    0    0    0    0    0    0   55   55   55    55    55    55    55
3    0    0    0    0    0    0    0    0    0    0     0     0     0     0

3. R `dlnm` 패키지

R의 dlnm 패키지를 이용하여 the Extended DLNM을 적합시키는 방법에 대해 소개하겠습니다. 사실, 함수의 input으로 시계열 데이터 대신 matrix of exposure histories를 입력한다는 부분을 제외하면 기본 DLNM을 적합시키는 방법과 크게 다르지 않습니다.

먼저 R에서 dlnm 패키지를 불러옵니다.

# install.packages("dlnm")
library(dlnm)

3.1 Data

the Extended DLNM의 R 실습에서 사용될 데이터는 dlnm 패키지에 포함되어 있는 drug와 nested 데이터셋입니다.

drug 데이터는 환자별 시간에 따른 약물 복용량의 효과에 대한 데이터입니다. 해당 실험은 200명의 무작위 추출된 환자들에 대해 진행되었으며, 각 환자들은 4주 중 무작위로 선정된 2주동안 투약을 받게 되고, 약물 복용량은 매주 달라집니다. 아래 데이터를 보면, 각 환자별로 주별 복용량(day1.7, day8.14, day15.21, day22.28)이 기록되어 있으며, out에 28일째에 측정된 outcome 값이 기록되어 있습니다.

head(drug, 5)

  id out sex day1.7 day8.14 day15.21 day22.28
1  1  46   M      0       0       40       37
2  2  50   F      0      47       55        0
3  3   7   F     56      22        0        0
4  4  70   M     91       0        0       87
5  5  -3   F      0      42       28        0

nested 데이터는 시간에 따른 노출 요인(exposure)과 암 사이의 연관성에 대한 nested case-control study 데이터입니다. 300명의 case 집단과 300명의 control 집단에 대한 데이터가 포함되어 있습니다. 아래 데이터의 case에 case(1)/control(0) 여부가 기록되어 있으며, exp15-exp60은 15세부터 65세까지 5년 간격으로 평균 노출 요인(exposure)을 구한 값입니다. 만약, 환자의 나이(age)가 컬럼에 해당하는 나이보다 적다면 해당 컬럼에는 NA가 입력됩니다. 예를 들어, 아래 데이터의 4번 환자는 52세이기 때문에 exp55 이후의 컬럼에는 NA가 입력되어 있습니다.

head(nested, 5)

  id case age riskset exp15 exp20 exp25 exp30 exp35 exp40 exp45 exp50 exp55
1  1    1  81     240     5    84    34    45   128    81    14    52    11
2  2    1  69     129    11     8    25     6     8    12    19    60    16
3  3    1  73     180    14    15     7    69    10   143    18    19    44
4  4    0  52      19    10    16     5    30    24    33    14   122    NA
5  4    0  66      96    10    16     5    30    24    33    14   122     2
  exp60
1    16
2    10
3    23
4    NA
5    11

The matrix of exposure histories

DLNM을 적합시키기 전, 데이터의 형태를 matrix of exposure histories 형태로 바꿔주어야 합니다.

drug 데이터의 matrix of exposure histories에서 lag0에 해당하는 값은 28일째의 약물 복용량입니다. 즉, lag0-lag6은 마지막 주의 약물 복용량, lag7-lag13은 셋째 주의 약물 복용량입니다. drug 데이터의 주별 약물 복용량 값을 7번씩 반복하여 다음과 같은 matrix of exposure histories를 만들어줍니다.

Qdrug <- as.matrix(drug[,rep(7:4, each=7)])
colnames(Qdrug) <- paste("lag", 0:27, sep="")
Qdrug[1:3,1:14]

  lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1   37   37   37   37   37   37   37   40   40   40    40    40    40    40
2    0    0    0    0    0    0    0   55   55   55    55    55    55    55
3    0    0    0    0    0    0    0    0    0    0     0     0     0     0

nested 데이터의 matrix of exposure histories에서 lag0에 해당하는 시점은 환자의 나이마다 다릅니다. 예를 들어, 52세 환자의 lag0 시점은 52년의 값이고, 65세 환자의 lag0 시점은 65년의 값입니다. 이런 경우, matrix of exposure histories를 만들기 까다로운데, dlnm 패키지에서는 이런 경우의 matrix of exposure histories를 만들어주는 exphist()라는 함수를 제공합니다. exphist() 함수는 다음과 같이 이용합니다.

exphist(exp, times, lag, fill=0)

exp에는 exposure profile(관측 첫 시점부터 매 시점의 exposure 값)을 입력합니다. times에는 lag0에 해당하는 시점을 입력합니다. times에 입력된 시점부터 시점을 거슬러가며 각 시차별 exposure 값이 구해집니다. lag에는 maximum lag 또는 lag range를 입력합니다. fill에는 어떤 시차에 해당하는 exposure 값이 존재하지 않을 때 채울 값을 입력합니다.

exphist() 함수를 이용하여 다음과 같은 matrix of exposure histories를 만들어줍니다. nested 데이터의 평균 노출 요인 값에 0, 0, 0을 추가하고 각 데이터를 5번씩 반복하여 1년부터 65년까지에 해당하는 exposure profile을 만들고 이 값을 exp에 넣습니다. times에는 각 환자별 나이(age)를 넣습니다. nested 데이터를 이용한 분석에서는 lag3부터 lag40까지의 데이터만 사용할 예정이기 때문에 lag에는 다음과 같이 lag range를 벡터로 넣어줍니다. nested 데이터의 matrix of exposure histories는 다음과 같습니다.

Qnest <- t(apply(nested, 1, function(sub) exphist(rep(c(0,0,0,sub[5:14]),
                                                      each=5), sub["age"], lag=c(3,40))))
colnames(Qnest) <- paste("lag", 3:40, sep="")
Qnest[1:3,1:11]

  lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1    0    0    0    0    0    0    0     0     0     0     0
2    0   10   10   10   10   10   16    16    16    16    16
3    0    0    0    0    0   23   23    23    23    23    44

3.2 A simple DLM

drug 데이터와 dlnm 패키지의 함수를 이용하여 약물 복용량과 outcome 사이의 관계에 대해 분석해보고자 합니다.

먼저, 약물 복용량에 대한 cross-basis matrix를 만들어야 합니다. crossbasis() 함수를 이용합니다. crossbasis()의 첫 번째 인수에는 matrix of exposure histories를 넣어줍니다. lag에는 lag-response 차원에서 시차를 얼마나 고려할지 lag period를 입력합니다. maximum lag 또는 lag range를 입력해야하며, minimum lag는 기본값이 0으로 지정되어 있습니다. 이때, lag period는 matrix of exposure histories의 시차 범위와 일치해야 합니다. argvar에는 exposure-response 차원에서 사용할 basis function, arglag에는 lag-response 차원에서 사용할 basis function을 입력합니다.

drug 데이터의 cross-basis는, exposure-response 차원에서 simple linear function을 이용하고, lag-response 차원에서 natural cubic spline을 이용합니다. 이 때, natural cubic spline에서 lag9, 18을 knots로 이용하는데, 이는 lag9, 18을 기준으로 구간을 나눠 각 구간에서 cubic regression model을 적합시킨다는 의미입니다.

summary() 함수를 통해 cross-basis matrix의 세부 사항을 확인할 수 있습니다.

cbdrug <- crossbasis(Qdrug, lag=27, argvar=list("lin"),
                     arglag=list(fun="ns",knots=c(9,18)))
summary(cbdrug)

CROSSBASIS FUNCTIONS
observations: 200 
range: 0 to 100 
lag period: 0 27 
total df:  4 

BASIS FOR VAR:
fun: lin 
intercept: FALSE 

BASIS FOR LAG:
fun: ns 
knots: 9 18 
intercept: TRUE 
Boundary.knots: 0 27

위에서 생성된 cross-basis 객체인 cbdrug를 회귀식에 포함시키고, 성별에 대한 효과를 보정하여 단순 선형 회귀 분석을 진행합니다. 모형을 통해 추정된 약물 복용량 및 그 시차에 대한 효과를 crosspred() 함수를 통해 확인할 수 있습니다. crosspred()의 첫 번째 인수에는 cross-basis 객체를, 두 번째 인수에는 cross-basis 객체를 사용한 모델을 입력합니다. 다음 코드에서, crosspred()의 at은 약물 복용량이 0:20*5(즉, 0, 5, 10, 15, …, 100)일 때의 outcome을 예측하라는 뜻입니다.

mdrug <- lm(out~cbdrug+sex, drug)
pdrug <- crosspred(cbdrug, mdrug, at=0:20*5)

crosspred 객체인 pdrug에 저장된 effect summaries는 다음과 같이 추출될 수 있습니다.

with(pdrug,cbind(allfit,alllow,allhigh)["50",])

  allfit   alllow  allhigh 
30.29584 20.12871 40.46298

위 코드는 약물 복용량이 50일때의 overall cumulative effects 추정치(allfit)와 95% 신뢰구간(alllow, allhigh) 추출한 것입니다. 이때, overall cumulative effects는 lag period인 28일동안 약물 복용량이 50으로 유지되었을 때 outcome의 전체적인 증가량 또는 약물 복용량 50이 28일 뒤 미치는 총 영향으로 해석될 수 있습니다.

위에서처럼 all-로 시작하는 객체들을 통해서는 전반적인 추정값을 확인할 수 있고, mat-으로 시작하는 객체들을 통해서는 특정 약물 복용량과 시차의 조합에 따른 추정 결과를 확인할 수 있습니다. 다음 코드는 lag3에서의 약물 복용량이 20일 때, outcome의 증가량을 추출한 것입니다. ?crosspred를 통해 crosspred()의 더 많은 기능을 확인할 수 있습니다.

pdrug$matfit["20","lag3"]

[1] 1.118139

plot() 함수를 이용해 추정 결과를 시각화할 수 있습니다.

par(mfrow=c(1,3))
plot(pdrug, zlab="Effect", xlab="Dose", ylab="Lag (days)")
plot(pdrug, var=60, ylab="Effect at dose 60", xlab="Lag (days)", ylim=c(-1,5))
plot(pdrug, lag=10, ylab="Effect at lag 10", xlab="Dose", ylim=c(-1,5))

첫 번째 그래프는 회귀 모델을 통해 추정한 exposure-lag-response 관계를 3차원 공간에 그려놓은 것입니다. 약물 복용량과 시차의 변화에 따라 effect가 어떻게 달라지는지 확인할 수 있습니다. 그래프에 따르면 약물 복용량의 효과는 복용 후 첫 번째 날에 나타나고 15-20일 후에 사라지는 경향이 있음을 알 수 있습니다.

두 번째와 세 번째 그래프는 약물 복용량이 60일 때의 lag-response curve와 lag가 10일 때 exposure-response curve를 그린 것입니다. 각각 var=60, lag=10을 지정하여 3차원상의 첫 번째 그래프의 단면을 자른 것과 같습니다. cross-basis matrix를 만들 때 basis function을 지정했던대로, lag-response 차원에서는 natural cubic spline, exposure-response 차원에서는 simple linear function 형태로 나타나는 것을 확인할 수 있습니다.

3.3 A more complex DLNM

neseted 데이터와 dlnm 패키지의 함수를 이용하여 노출 요인(exposure)에 대한 장기간의 노출이 암에 어떤 영향을 미치는지 분석하고자 합니다.

최근 3년(lag0-2)간의 노출은 암에 영향을 주지 않는다고 가정해봅시다. 따라서 matrix of exposure histories에는 lag3부터의 데이터만 있으면 됩니다.

먼저, cross-basis matrix를 만듭니다. crossbasis() 함수를 이용합니다. lag3부터 lag40까지의 데이터만 이용할 것이므로, lag에 c(3,40)을 넣어줍니다. exposure-response 차원에서는 basis 함수로 quadratic spline을 사용하고, exposure-lag 차원에서는 natural cubic spline을 사용합니다. quadratic spline의 자유도(df)와 차수(degree)를 입력해야하며, single knot는 별도로 지정하지 않으면 중앙값으로 지정됩니다. natural cubic spline은 intercept=F를 입력하여 intercept를 제외합니다. 위에서 lag0-2는 고려하지 않는다고 하였으므로, intercept를 제외함으로써 위 가정과 일치하게 시차 차원에서의 null effect를 예측할 수 있습니다.

cbnest <- crossbasis(Qnest, lag=c(3,40), argvar=list("bs",degree=2,df=3),
                     arglag=list(fun="ns",knots=c(10,30),intercept=F))
summary(cbnest)

CROSSBASIS FUNCTIONS
observations: 600 
range: 0 to 1064 
lag period: 3 40 
total df:  9 

BASIS FOR VAR:
fun: bs 
knots: 15 
degree: 2 
intercept: FALSE 
Boundary.knots: 0 1064 

BASIS FOR LAG:
fun: ns 
knots: 10 30 
intercept: FALSE 
Boundary.knots: 3 40

앞에서와 마찬가지로, cross-basis 객체를 회귀 모형에 포함시킵니다. nested 데이터는 nested case-control study 데이터이기 때문에, conditional logistic regression을 이용합니다. survival 패키지의 clogit() 함수를 이용합니다.

crosspred 객체를 만들 때, cen=0을 입력하여 reference value를 0으로 지정합니다. 즉, 노출 요인에 따른 암의 OR은 exposure=0일 때를 기준으로 계산됩니다.

library(survival)
mnest <- clogit(case~cbnest+strata(riskset), nested)
pnest <- crosspred(cbnest, mnest, cen=0, at=0:20*5)

effect summaries는 pnest를 통해 확인할 수 있습니다. 이번에는 allRR-이나 matRR-로 시작하는 객체를 추출하여 OR 추정값을 확인할 수 있습니다.

with(pnest,cbind(allRRfit,allRRlow,allRRhigh)["50",])

  allRRfit   allRRlow  allRRhigh 
 32.061676   2.561702 401.276652

pnest$matRRfit["50","lag5"]

[1] 1.058661

앞에서와 마찬가지로, plot() 함수를 통해 추정 결과를 시각화할 수 있습니다.

par(mfrow=c(1,3))
plot(pnest, zlab="OR", xlab="Exposure", ylab="Lag (years)")
plot(pnest, var=50, ylab="OR for exposure 50", xlab="Lag (years)", xlim=c(0,40))
plot(pnest, lag=5, ylab="OR at lag 5", xlab="Exposure", ylim=c(0.95,1.15))

첫 번째 그래프는 노출 요인과 암 사이의 exposure-lag-response 관계를 3차원 상에서 시각화한 것입니다. 그래프에 따르면 노출 초기엔 암의 위험성을 증가시키다가 점차 감소한다는 사실을 확인할 수 있습니다.

두 번째와 세 번째 그래프는 exposure이 50일 때의 lag-response curve와 lag가 5일 때 exposure-response curve를 그린 것입니다. 두 번째 그래프를 통해 노출 후 10-15년이 지났을 때 암의 위험성이 가장 크게 증가하며, 30년이 지났을 때쯤부턴 거의 원 상태로 되돌아온다는 사실을 확인할 수 있습니다. 세 번째 그림을 통해 노출 요인의 양이 20으로 증가할 때까지는 위험성이 크게 증가하지만 그 이후로는 위험성이 천천히 증가함을 확인할 수 있습니다.

3.4 Extended prediction summaries

crosspred()의 at과 lag 인자를 통해 특정 exposure 값 또는 lag에서의 effect summaries를 얻을 수 있습니다. 이때, nested 데이터의 matrix of exposure histories를 만들 때 사용한 함수인 exphist()를 사용한다면, 특정 exposure history에서의 effect summaries를 얻을 수 있습니다.

예를 들어, ’3.3 A more complex DLNM’에서 이용한 nested 데이터를 다시 분석해봅시다. 5년 동안의 exposure 값이 10이고, 그 이후 5년 간은 0, 그 이후 10년간은 13일 때의 effect summary를 구하고자 할 때, 아래와 같이 exposure history를 만들 수 있습니다.

expnested <- rep(c(10,0,13), c(5,5,10))
hist <- exphist(expnested, time=length(expnested), lag=c(3,40))
hist

   lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13 lag14 lag15 lag16
20   13   13   13   13   13   13   13     0     0     0     0     0    10    10
   lag17 lag18 lag19 lag20 lag21 lag22 lag23 lag24 lag25 lag26 lag27 lag28
20    10    10    10     0     0     0     0     0     0     0     0     0
   lag29 lag30 lag31 lag32 lag33 lag34 lag35 lag36 lag37 lag38 lag39 lag40
20     0     0     0     0     0     0     0     0     0     0     0     0

time 인자에는 expnested의 길이를 입력하고, lag는 lag3부터 40까지만 고려한다고 입력합니다. 이때, expnested의 길이는 20이기 때문에 lag20-40의 exposure 값은 0으로 입력됩니다.

이렇게 만들어진 exposure history를 crosspred()의 at 인자에 넣어줍니다.

pnesthist <- crosspred(cbnest, mnest, cen=0, at=hist)
with(pnesthist, c(allRRfit,allRRlow,allRRhigh))

      20       20       20 
3.503928 1.240109 9.900351

위의 코드를 통해 해당 exposure history에서의 OR 추정치는 3.5임을 알 수 있습니다.

이 방법을 통해 여러 개의 time-varying exposure histories를 만들어 effect summaries를 구할 수 있습니다. ’3.2 A simple DLM’에서 이용한 drug 데이터를 다시 분석해봅시다. 환자의 2주 동안의 약물 복용량이 10이고, 그 후 1주간은 50이고, 그 후 1주간은 복용하지 않았을 때, 각 시점에서의 effect summary를 구해봅시다. 아래와 같이 각 시점에서의 exposure history을 만들어줍니다.

expdrug <- rep(c(10,50,0,20),c(2,1,1,2)*7)
dynhist <- exphist(expdrug, lag=27)
dynhist[1:10,]

   lag0 lag1 lag2 lag3 lag4 lag5 lag6 lag7 lag8 lag9 lag10 lag11 lag12 lag13
1    10    0    0    0    0    0    0    0    0    0     0     0     0     0
2    10   10    0    0    0    0    0    0    0    0     0     0     0     0
3    10   10   10    0    0    0    0    0    0    0     0     0     0     0
4    10   10   10   10    0    0    0    0    0    0     0     0     0     0
5    10   10   10   10   10    0    0    0    0    0     0     0     0     0
6    10   10   10   10   10   10    0    0    0    0     0     0     0     0
7    10   10   10   10   10   10   10    0    0    0     0     0     0     0
8    10   10   10   10   10   10   10   10    0    0     0     0     0     0
9    10   10   10   10   10   10   10   10   10    0     0     0     0     0
10   10   10   10   10   10   10   10   10   10   10     0     0     0     0
   lag14 lag15 lag16 lag17 lag18 lag19 lag20 lag21 lag22 lag23 lag24 lag25
1      0     0     0     0     0     0     0     0     0     0     0     0
2      0     0     0     0     0     0     0     0     0     0     0     0
3      0     0     0     0     0     0     0     0     0     0     0     0
4      0     0     0     0     0     0     0     0     0     0     0     0
5      0     0     0     0     0     0     0     0     0     0     0     0
6      0     0     0     0     0     0     0     0     0     0     0     0
7      0     0     0     0     0     0     0     0     0     0     0     0
8      0     0     0     0     0     0     0     0     0     0     0     0
9      0     0     0     0     0     0     0     0     0     0     0     0
10     0     0     0     0     0     0     0     0     0     0     0     0
   lag26 lag27
1      0     0
2      0     0
3      0     0
4      0     0
5      0     0
6      0     0
7      0     0
8      0     0
9      0     0
10     0     0

exphist() 함수의 time 인자가 지정되지 않았을 경우에는 모든 time point에 대하여 exposure history가 만들어집니다. 이렇게 만들어진 exposure histories를 crosspred()의 at 인자에 넣어줍니다.

pdyndrug <- crosspred(cbdrug, mdrug, at=dynhist)

아래 코드로 그린 plot을 통해 각 시점에서의 exposure history의 effect summary를 확인할 수 있습니다.

plot(pdyndrug,"overall", ylab="Effect", xlab="Time (days)", ylim=c(-10,27),
     xlim=c(1,50), yaxt="n")
axis(2, at=-1:5*5)
par(new=TRUE)
plot(expdrug, type="h", xlim=c(1,50), ylim=c(0,300), axes=F, ann=F)
axis(4, at=0:6*10, cex.axis=0.8)
mtext("Dose", 4, line=-1.5, at=30, cex=0.8)

3.5 Applying user-defined functions

cross-basis 객체를 생성할 때, exposure와 lag 차원에서 사용할 basis function을 사용자 정의 함수로 지정할 수 있습니다. 이때, 사용자 정의 함수의 첫 번째 인자는 반드시 \(X\)값이어야하고, return 값은 변환된 벡터 혹은 행렬이어야 합니다.

’3.3 A more complex DLNM’의 세 번째 plot을 보면, exposure-response 차원에서의 quadratic spline은 log 함수와 비슷한 형태임을 알 수 있습니다. cross-basis 객체를 생성할 때, exposure-response 차원에서의 basis function을 quadratic spline이 아닌 log 함수를 이용해봅시다.

다음과 같이 log 함수를 정의합니다.

mylog <- function(x) log(x+1)

위에서 정의한 mylog 함수를 crossbasis() 함수에 넣어줍니다.

cbnest2 <- crossbasis(Qnest, lag=c(3,40), argvar=list("mylog"),
                      arglag=list(fun="ns",knots=c(10,30),intercept=F))
summary(cbnest2)

CROSSBASIS FUNCTIONS
observations: 600 
range: 0 to 1064 
lag period: 3 40 
total df:  3 

BASIS FOR VAR:
fun: mylog 

BASIS FOR LAG:
fun: ns 
knots: 10 30 
intercept: FALSE 
Boundary.knots: 3 40

3.3에서 정의한 cross-basis cbnest와 비교하여 자유도가 9에서 3으로 감소한 것을 확인할 수 있습니다.

아래 plot을 통해 추정 결과를 비교해봅시다.

mnest2 <- clogit(case~cbnest2+strata(riskset), nested)
pnest2 <- crosspred(cbnest2, mnest2, cen=0, at=0:20*5)

par(mfrow=c(1,3))
plot(pnest2, zlab="OR", xlab="Exposure", ylab="Lag (years)")
plot(pnest2, var=50, ylab="OR for exposure 50", xlab="Lag (years)", xlim=c(0,40))
lines(pnest, var=50, lty=2)
plot(pnest2, lag=5, ylab="OR at lag 5", xlab="Exposure", ylim=c(0.95,1.15))
lines(pnest, lag=5, lty=2)

점선은 3.3에서 추정한 결과를 나타낸 것입니다. exposure-response 차원에서 basis function으로 log 함수를 사용했을 때와 quadratic spline을 사용했을 때의 결과가 유사함을 확인할 수 있습니다.

’3.2 A simple DLM’의 두 번째 plot을 보면, lag-response 차원에서의 natural cubic spline은 지수적으로 감소하는 형태를 보입니다. 다음과 같이 exponential decay 함수를 정의합니다.

fdecay <- function(x,scale=5) {
  basis <- exp(-x/scale)
  attributes(basis)$scale <- scale
  return(basis)
}

다음과 같이 crossbasis() 함수에 인자로 넣어주고, plot을 그려 결과를 확인합니다.

cbdrug2 <- crossbasis(Qdrug, lag=27, argvar=list("lin"),
                      arglag=list(fun="fdecay",scale=6))
summary(cbdrug2)

CROSSBASIS FUNCTIONS
observations: 200 
range: 0 to 100 
lag period: 0 27 
total df:  1 

BASIS FOR VAR:
fun: lin 
intercept: FALSE 

BASIS FOR LAG:
fun: fdecay 
scale: 6

mdrug2 <- lm(out~cbdrug2+sex, drug)
pdrug2 <- crosspred(cbdrug2, mdrug2, at=0:20*5)

par(mfrow=c(1,3))
plot(pdrug2, zlab="Effect", xlab="Dose", ylab="Lag (days)")
plot(pdrug2, var=60, ylab="Effect at dose 60", xlab="Lag (days)", ylim=c(-1,5))
lines(pdrug, var=60, lty=2)
plot(pdrug2, lag=10, ylab="Effect at lag 10", xlab="Dose", ylim=c(-1,5))
lines(pdrug, lag=10, lty=2)

점선은 3.2에서 추정한 결과를 나타낸 것입니다. lag-response 차원에서 basis function으로 exponential decay 함수를 사용했을 때와 natural cubic spline을 사용했을 때의 결과가 유사함을 확인할 수 있습니다.

References

Gasparrini A, Armstrong B, Kenward MG. Distributed lag non-linear models. Stat Med. 2010 Sep 20;29(21):2224-34. doi: 10.1002/sim.3940. PMID: 20812303; PMCID: PMC2998707.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2998707/pdf/sim0029-2224.pdf
https://cran.r-project.org/web/packages/dlnm/vignettes/dlnmExtended.pdf

Citation

BibTeX citation:

@online{kwon2024,
  author = {Kwon, Suhyeon},
  title = {The {Extended} {DLNM} {소개}},
  date = {2024-06-05},
  url = {https://blog.zarathu.com/posts/2024-06-05-theExtendedDLNM/},
  langid = {en}
}

For attribution, please cite this work as:

Kwon, Suhyeon. 2024. “The Extended DLNM 소개.” June 5, 2024. https://blog.zarathu.com/posts/2024-06-05-theExtendedDLNM/.