Treed Distributed Lag Mixture Models
Daniel Mork, Colorado State University
11/17/2020
Introduction to TDLNM
In this vignette we demonstrate how to apply the tdlmm function found in the dlmtree package. TDLMM is a framework for estimating distributed lag mixture models. This method operates using a set of regression tree-pairs that structure exposure measurements taken at high temporal resolution to estimate main effects of individual exposures and pairwise interaction effects of two exposures. The trees and exposures estimated by the trees are both stochastic to evaluate differences in distributed lag models and select exposures that are associated with the response.
Installation
To install package dlmtree first use the devtools package and run function install_github.
library(devtools)
install_github("danielmork/dlmtree")Data Simulation
For this vignette we will consider a data simulation available in the dlmtree package. The exposure data is produced following the covariance structure of real exposure data.
Scenario 1
In this scenario, we produce a binary reponse to a single exposure. The exposure has a random, contiguous 8-week window of effect within the exposure timeframe (37 weeks).
library(dlmtree)
set.seed(1)
D <- sim.tdlmm(sim = 1, mean.p = 0.5, n = 5000)TDLM (single exposure model)
We first run TDLM, which estimates a distributed lag model for a single exposure. We recommend using \(20\) trees and \(\geq15000\) iterations thinned to every fifth iteration following \(\geq5000\) burn-in iterations. Model convergence can be checked by comparing the consitency or results across multiple runs.
m1 <- tdlnm(y ~ .,
data = D$dat,
exposure.data = D$exposures[[1]],
exposure.splits = 0, # make this a DLM (as opposed to DLNM)
family = "logit", # change to 'gaussian' for continuous response
binomial.size = 1,
n.trees = 20, n.burn = 1000, n.iter = 2000, n.thin = 2)## Preparing data...
## Burn-in % complete
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## MCMC iterations (est time: 32 seconds)
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
##
## Compiling results...Summary
Following the model run, summarize the results using the summary function.
m1_sum <- summary(m1, conf.level = 0.95)
m1_sum## ---
## TDLM summary
##
## Model run info
## - trees
## - 1000 burn-in iterations
## - 2000 post-burn iterations
## - 2 thinning factor
## - 0.95 confidence level
##
## Fixed effect coefficients:
## Mean Lower.Bound Upper.Bound
## *(Intercept) 4.271 0.315 8.477
## *c1 -0.109 -0.167 -0.051
## *c2 -0.144 -0.203 -0.083
## c3 0.011 -0.047 0.069
## c4 -0.036 -0.096 0.023
## *c5 0.129 0.072 0.185
## b1 -0.033 -0.140 0.081
## b2 -0.041 -0.157 0.074
## b3 0.102 -0.017 0.221
## *b4 0.116 0.002 0.233
## b5 0.017 -0.098 0.130
## ---
## * = CI does not contain zero
##
## DLM effect:
## range = [-0.011, 0.088]
## critical windows: 20 21 22 23 24 25 26 27
## Mean Lower.Bound Upper.Bound
## Period 1 -0.004 -0.033 0.018
## Period 2 -0.001 -0.019 0.020
## Period 3 -0.001 -0.017 0.017
## Period 4 -0.002 -0.016 0.017
## Period 5 -0.003 -0.020 0.013
## Period 6 -0.003 -0.018 0.013
## Period 7 -0.004 -0.020 0.012
## Period 8 -0.004 -0.019 0.013
## Period 9 -0.005 -0.020 0.014
## Period 10 -0.007 -0.026 0.008
## Period 11 -0.008 -0.025 0.006
## Period 12 -0.009 -0.029 0.006
## Period 13 -0.009 -0.030 0.006
## Period 14 -0.009 -0.030 0.007
## Period 15 -0.008 -0.029 0.009
## Period 16 -0.006 -0.027 0.014
## Period 17 0.001 -0.020 0.036
## Period 18 0.006 -0.017 0.050
## Period 19 0.025 -0.013 0.085
## *Period 20 0.081 0.043 0.107
## *Period 21 0.083 0.053 0.106
## *Period 22 0.086 0.065 0.107
## *Period 23 0.087 0.068 0.107
## *Period 24 0.088 0.069 0.109
## *Period 25 0.088 0.068 0.111
## *Period 26 0.087 0.064 0.111
## *Period 27 0.069 0.002 0.102
## Period 28 0.053 -0.013 0.096
## Period 29 -0.010 -0.038 0.023
## Period 30 -0.011 -0.037 0.011
## Period 31 -0.011 -0.035 0.009
## Period 32 -0.008 -0.027 0.011
## Period 33 -0.007 -0.025 0.013
## Period 34 -0.006 -0.025 0.015
## Period 35 -0.006 -0.023 0.013
## Period 36 -0.006 -0.026 0.016
## Period 37 -0.005 -0.029 0.021
## ---
## * = CI does not contain zeroPlot the DLM using plot with the model summary.
plot(m1_sum, main = "Estimated DLM", xlab = "Time", ylab = "Effect")
Compare to true generated DLM.
# RMSE
sqrt(mean((D$margDLM - m1_sum$matfit)^2))## [1] 0.02068501# Coverage
mean((m1_sum$cilower < D$margDLM) & (m1_sum$ciupper > D$margDLM))## [1] 0.9459459# True positive effect classification
sum((m1_sum$cilower > 0) & (D$margDLM > 0)) / sum(D$margDLM > 0)## [1] 0.875# False positive effect classification
sum(((m1_sum$cilower > 0) | (m1_sum$ciupper < 0)) & (D$margDLM == 0)) /
sum(D$margDLM == 0)## [1] 0.03448276TDLMM (multi exposure model)
Next we demonstrate, which estimates a distributed lag mixture model for many exposure. This demonstrates TDLMMns, which does not estimate within-exposure interactions. We recommend using \(20\) trees and \(\geq15000\) iterations thinned to every fifth iteration following \(\geq5000\) burn-in iterations. Model convergence can be checked by comparing the consitency or results across multiple runs.
m2 <- tdlmm(y ~ .,
data = D$dat,
exposure.data = D$exposures,
mixture.interactions = "noself", # also try 'all' and 'none'
family = "logit", # change to 'gaussian' for continuous response
binomial.size = 1,
n.trees = 20, n.burn = 1000, n.iter = 2000, n.thin = 2)## Preparing data...
## Burn-in % complete
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## MCMC iterations (est time: 1.8 minutes)
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## Compiling results...Summary
Following the model run, summarize the results using the summary function.
m2_sum <- summary(m2, conf.level = 0.95)## Reconstructing main effects...
## Reconstructing interaction effects...
## e1-e2 e1-e3 e1-e4 e1-e5 e2-e2 e2-e3 e2-e4 e2-e5 e3-e2 e3-e3 e3-e4 e3-e5 e4-e2 e4-e3 e4-e4 e4-e5 Calculating marginal effects...print(m2_sum, cw.only = F) # show results for all exposures and interactions##
## TDLMM:
##
## Model run info
## - 20 trees
## - 1000 burn-in iterations
## - 2000 post-burn iterations
## - 2 thinning factor
## - 5 exposures
## - 10 two-way interactions
## - 0.95 confidence level
##
## Fixed effects:
## Mean Lower Upper
## (Intercept) 4.1288 -0.3074 8.6319
## *c1 -0.1099 -0.1678 -0.0523
## *c2 -0.1461 -0.2017 -0.0895
## c3 0.0089 -0.0487 0.0669
## c4 -0.0366 -0.0977 0.0210
## *c5 0.1296 0.0737 0.1849
## b1 -0.0335 -0.1479 0.0870
## b2 -0.0435 -0.1608 0.0717
## b3 0.1037 -0.0166 0.2136
## *b4 0.1192 0.0096 0.2404
## b5 0.0243 -0.0878 0.1365
## ---
## * = CI does not contain zero
##
## --
## Exposure effects: critical windows
##
## *e1 (1): 20-27
## e2 (0.29):
## e3 (0.41):
## e4 (0.35):
## e5 (0.46):
## --
## Interaction effects: critical windows
##
## e1/e2 (0.49):
## e1/e3 (0.61):
## e1/e4 (0.56):
## e1/e5 (0.81):
## e2/e3 (0.44):
## e2/e4 (0.37):
## e2/e5 (0.43):
## e3/e4 (0.39):
## e3/e5 (0.47):
## e4/e5 (0.43):
## ---
## * = Model selected exposure or interaction
## (0.xx) = Relative signal strengthPlot the marginalized DLM (integrating over other exposures) using plot with the model summary.
plot(m2_sum, exposure1 = "e1",
main = "Estimated marginal DLM for active exposure", xlab = "Time", ylab = "Effect")
plot(m2_sum, exposure1 = "e2",
main = "Estimated marginal DLM for a non-active exposure", xlab = "Time", ylab = "Effect")
Compare to true generated DLM.
# RMSE
sqrt(mean((D$margDLM - m2_sum$DLM$e1$marg.matfit)^2))## [1] 0.0213697# Coverage
mean((m2_sum$DLM$e1$marg.cilower < D$margDLM) &
(m2_sum$DLM$e1$marg.ciupper > D$margDLM))## [1] 0.9459459# True positive effect classification
sum((m2_sum$DLM$e1$marg.cilower > 0) & (D$margDLM > 0)) / sum(D$margDLM > 0)## [1] 0.875# False positive effect classification
sum(((m2_sum$DLM$e1$marg.cilower > 0) |
(m2_sum$DLM$e1$marg.ciupper < 0)) & (D$margDLM == 0)) /
sum(D$margDLM == 0)## [1] 0.03448276# Exposure selection
m2_sum$expSel## e1 e2 e3 e4 e5
## TRUE FALSE FALSE FALSE FALSE# Interaction selection
m2_sum$mixSel## e1-e2 e1-e3 e1-e4 e1-e5 e2-e3 e2-e4 e2-e5 e3-e4 e3-e5 e4-e5
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE# Exposure relative effect size
m2_sum$expVar[2,]## e1 e2 e3 e4 e5
## 0.999750 0.285125 0.410000 0.347625 0.457500# Interaction relative effect size
m2_sum$mixVar[2,]## e1-e2 e1-e3 e1-e4 e1-e5 e2-e3 e2-e4 e2-e5 e3-e4
## 0.4908333 0.6091111 0.5586667 0.8106667 0.4408889 0.3699444 0.4261667 0.3940000
## e3-e5 e4-e5
## 0.4705000 0.4292222Scenario 2
In this scenario, we produce a continuous reponse for a main effect from exposure 1 and an interation between exposures 1 and 2.
set.seed(1)
D <- sim.tdlmm(sim = 2, error = 25, n = 5000)TDLMMns (no self interactions)
Next we demonstrate, which estimates a distributed lag mixture model for many exposure. This demonstrates TDLMMns, which does not estimate within-exposure interactions.
m3 <- tdlmm(y ~ .,
data = D$dat,
exposure.data = D$exposures,
mixture.interactions = "noself", # also try 'all' and 'none'
n.trees = 20, n.burn = 1000, n.iter = 2000, n.thin = 2)## Preparing data...
## Burn-in % complete
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## MCMC iterations (est time: 48 seconds)
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## Compiling results...Summary
Following the model run, summarize the results using the summary function.
m3_sum <- summary(m3, conf.level = 0.95)## Reconstructing main effects...
## Reconstructing interaction effects...
## e1-e2 e1-e3 e1-e4 e1-e5 e2-e2 e2-e3 e2-e4 e2-e5 e3-e2 e3-e3 e3-e4 e3-e5 e4-e2 e4-e3 e4-e4 e4-e5 Calculating marginal effects...print(m3_sum, cw.only = F) # show results for all exposures and interactions##
## TDLMM:
##
## Model run info
## - 20 trees
## - 1000 burn-in iterations
## - 2000 post-burn iterations
## - 2 thinning factor
## - 5 exposures
## - 10 two-way interactions
## - 0.95 confidence level
##
## Fixed effects:
## Mean Lower Upper
## (Intercept) -9.1613 -19.6862 2.0401
## *c1 -1.1401 -1.2814 -1.0087
## *c2 -1.5762 -1.7235 -1.4445
## c3 -0.0060 -0.1473 0.1340
## c4 -0.0814 -0.2331 0.0604
## *c5 1.3049 1.1746 1.4368
## b1 0.1779 -0.0971 0.4604
## b2 -0.2052 -0.4910 0.0867
## *b3 0.8449 0.5581 1.1299
## *b4 1.0598 0.7783 1.3328
## *b5 0.6024 0.3071 0.8690
## ---
## * = CI does not contain zero
##
## --
## Exposure effects: critical windows
##
## *e1 (0.87): 19-27
## *e2 (0.79): 22-28
## e3 (0.21):
## e4 (0.36):
## e5 (0.26):
## --
## Interaction effects: critical windows
##
## *e1/e2 (1):
## 19/23
## 20/20-29
## 21/21-29
## 22/22-29
## 23/23-29
## 24/24-29
## 25/25-28
## 26/26-28
## 27/
## e1/e3 (0.42):
## e1/e4 (0.47):
## e1/e5 (0.43):
## e2/e3 (0.45):
## e2/e4 (0.53):
## e2/e5 (0.53):
## e3/e4 (0.39):
## e3/e5 (0.39):
## e4/e5 (0.39):
## ---
## * = Model selected exposure or interaction
## (0.xx) = Relative signal strengthPlot the DLM using plot with the model summary.
plot(m3_sum, exposure1 = "e1",
main = "Estimated marginal DLM for active exposure 1", xlab = "Time", ylab = "Effect")
plot(m3_sum, exposure1 = "e2",
main = "Estimated marginal DLM for a active exposure 2", xlab = "Time", ylab = "Effect")
plot(m3_sum, exposure1 = "e1", exposure2 = "e2")
Compare to true generated DLM.
# RMSE exposure 1
sqrt(mean((D$margDLM1 - m3_sum$DLM$e1$marg.matfit)^2))## [1] 0.05114957# RMSE exposure 2
sqrt(mean((D$margDLM2 - m3_sum$DLM$e2$marg.matfit)^2))## [1] 0.03032486# Coverage exposure 1
mean((m3_sum$DLM$e1$marg.cilower < D$margDLM1) &
(m3_sum$DLM$e1$marg.ciupper > D$margDLM1))## [1] 0.972973# Coverage exposure 2
mean((m3_sum$DLM$e2$marg.cilower < D$margDLM2) &
(m3_sum$DLM$e2$marg.ciupper > D$margDLM2))## [1] 1# True positive effect classification exposure 1
sum((m3_sum$DLM$e1$marg.cilower > 0) & (D$margDLM1 > 0)) /
sum(D$margDLM1 > 0)## [1] 1# True positive effect classification exposure 2
sum((m3_sum$DLM$e2$marg.cilower > 0) & (D$margDLM2 > 0)) /
sum(D$margDLM2 > 0)## [1] 0.875# False positive effect classification exposure 1
sum(((m3_sum$DLM$e1$marg.cilower > 0) |
(m3_sum$DLM$e1$marg.ciupper < 0)) & (D$margDLM1 == 0)) /
sum(D$margDLM1 == 0)## [1] 0.03448276# False positive effect classification exposure 2
sum(((m3_sum$DLM$e2$marg.cilower > 0) |
(m3_sum$DLM$e2$marg.ciupper < 0)) & (D$margDLM2 == 0)) /
sum(D$margDLM2 == 0)## [1] 0# Exposure selection
m3_sum$expSel## e1 e2 e3 e4 e5
## TRUE TRUE FALSE FALSE FALSE# Interaction selection
m3_sum$mixSel## e1-e2 e1-e3 e1-e4 e1-e5 e2-e3 e2-e4 e2-e5 e3-e4 e3-e5 e4-e5
## TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE# Exposure relative effect size
m3_sum$expVar[2,]## e1 e2 e3 e4 e5
## 0.87275 0.79400 0.21325 0.36000 0.26000# Interaction relative effect size
m3_sum$mixVar[2,]## e1-e2 e1-e3 e1-e4 e1-e5 e2-e3 e2-e4 e2-e5 e3-e4
## 0.9970000 0.4175556 0.4703889 0.4331111 0.4514444 0.5308889 0.5300000 0.3894444
## e3-e5 e4-e5
## 0.3861667 0.3940000TDLMM (all interactions)
This demonstrates TDLMM, which also estimates within-exposure interactions.
m4 <- tdlmm(y ~ .,
data = D$dat,
exposure.data = D$exposures,
mixture.interactions = "all",
n.trees = 20, n.burn = 1000, n.iter = 2000, n.thin = 2)## Preparing data...
## Burn-in % complete
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## MCMC iterations (est time: 1.3 minutes)
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## Compiling results...Summary
Following the model run, summarize the results using the summary function.
m4_sum <- summary(m4, conf.level = 0.95)## Reconstructing main effects...
## Reconstructing interaction effects...
## e1-e1 e1-e2 e1-e3 e1-e4 e1-e5 e2-e1 e2-e2 e2-e3 e2-e4 e2-e5 e3-e1 e3-e2 e3-e3 e3-e4 e3-e5 e4-e1 e4-e2 e4-e3 e4-e4 e4-e5 e5-e1 e5-e2 e5-e3 e5-e4 e5-e5 Calculating marginal effects...print(m4_sum, cw.only = F) # show results for all exposures and interactions##
## TDLMM:
##
## Model run info
## - 20 trees
## - 1000 burn-in iterations
## - 2000 post-burn iterations
## - 2 thinning factor
## - 5 exposures
## - 15 two-way interactions
## - 0.95 confidence level
##
## Fixed effects:
## Mean Lower Upper
## (Intercept) -6.9956 -17.7004 5.6167
## *c1 -1.1437 -1.2790 -1.0044
## *c2 -1.5703 -1.7104 -1.4390
## c3 -0.0030 -0.1510 0.1480
## c4 -0.0782 -0.2248 0.0632
## *c5 1.3020 1.1622 1.4356
## b1 0.1717 -0.1234 0.4715
## b2 -0.2020 -0.4826 0.0880
## *b3 0.8460 0.5590 1.1368
## *b4 1.0583 0.7935 1.3448
## *b5 0.6142 0.3358 0.8885
## ---
## * = CI does not contain zero
##
## --
## Exposure effects: critical windows
##
## *e1 (0.89): 20-26
## *e2 (0.56): 23-25,27
## e3 (0.25):
## e4 (0.32):
## e5 (0.48):
## --
## Interaction effects: critical windows
##
## *e1/e1 (0.88):
## *e1/e2 (0.91):
## e1/e3 (0.51):
## e1/e4 (0.49):
## e1/e5 (0.75):
## e2/e2 (0.52):
## e2/e3 (0.4):
## e2/e4 (0.43):
## e2/e5 (0.5):
## e3/e3 (0.33):
## e3/e4 (0.34):
## e3/e5 (0.36):
## e4/e4 (0.35):
## e4/e5 (0.36):
## e5/e5 (0.38):
## ---
## * = Model selected exposure or interaction
## (0.xx) = Relative signal strengthPlot the DLM using plot with the model summary.
plot(m4_sum, exposure1 = "e1",
main = "Estimated marginal DLM for active exposure 1", xlab = "Time", ylab = "Effect")
plot(m4_sum, exposure1 = "e2",
main = "Estimated marginal DLM for a active exposure 2", xlab = "Time", ylab = "Effect")
plot(m4_sum, exposure1 = "e1", exposure2 = "e2")
Compare to true generated DLM.
# RMSE exposure 1
sqrt(mean((D$margDLM1 - m4_sum$DLM$e1$marg.matfit)^2))## [1] 0.04512961# RMSE exposure 2
sqrt(mean((D$margDLM2 - m4_sum$DLM$e2$marg.matfit)^2))## [1] 0.02910858# Coverage exposure 1
mean((m4_sum$DLM$e1$marg.cilower < D$margDLM1) &
(m4_sum$DLM$e1$marg.ciupper > D$margDLM1))## [1] 1# Coverage exposure 2
mean((m4_sum$DLM$e2$marg.cilower < D$margDLM2) &
(m4_sum$DLM$e2$marg.ciupper > D$margDLM2))## [1] 1# True positive effect classification exposure 1
sum((m4_sum$DLM$e1$marg.cilower > 0) & (D$margDLM1 > 0)) /
sum(D$margDLM1 > 0)## [1] 0.875# True positive effect classification exposure 2
sum((m4_sum$DLM$e2$marg.cilower > 0) & (D$margDLM2 > 0)) /
sum(D$margDLM2 > 0)## [1] 0.5# False positive effect classification exposure 1
sum(((m4_sum$DLM$e1$marg.cilower > 0) |
(m4_sum$DLM$e1$marg.ciupper < 0)) & (D$margDLM1 == 0)) /
sum(D$margDLM1 == 0)## [1] 0# False positive effect classification exposure 2
sum(((m4_sum$DLM$e2$marg.cilower > 0) |
(m4_sum$DLM$e2$marg.ciupper < 0)) & (D$margDLM2 == 0)) /
sum(D$margDLM2 == 0)## [1] 0# Exposure selection
m4_sum$expSel## e1 e2 e3 e4 e5
## TRUE TRUE FALSE FALSE FALSE# Interaction selection
m4_sum$mixSel## e1-e1 e1-e2 e1-e3 e1-e4 e1-e5 e2-e2 e2-e3 e2-e4 e2-e5 e3-e3 e3-e4 e3-e5 e4-e4
## TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## e4-e5 e5-e5
## FALSE FALSE# Exposure relative effect size
m4_sum$expVar[2,]## e1 e2 e3 e4 e5
## 0.88975 0.55900 0.24600 0.32125 0.48400# Interaction relative effect size
m4_sum$mixVar[2,]## e1-e1 e1-e2 e1-e3 e1-e4 e1-e5 e2-e2 e2-e3 e2-e4
## 0.8835357 0.9063571 0.5061429 0.4909643 0.7511786 0.5155000 0.4024643 0.4308571
## e2-e5 e3-e3 e3-e4 e3-e5 e4-e4 e4-e5 e5-e5
## 0.4967857 0.3261429 0.3373571 0.3618214 0.3522500 0.3612500 0.3773929TDLMMadd (no interactions)
This demonstrates TDLMMadd, which is an additive distributed lag model.
m5 <- tdlmm(y ~ .,
data = D$dat,
exposure.data = D$exposures,
mixture.interactions = "none",
n.trees = 20, n.burn = 1000, n.iter = 2000, n.thin = 2)## Preparing data...
## Burn-in % complete
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## MCMC iterations (est time: 20 seconds)
## [0--------25--------50--------75--------100]
## ''''''''''''''''''''''''''''''''''''''''''
## Compiling results...Summary
Following the model run, summarize the results using the summary function.
m5_sum <- summary(m5, conf.level = 0.95)## Reconstructing main effects...
## Reconstructing interaction effects...
## Calculating marginal effects...print(m5_sum, cw.only = F) # show results for all exposures and interactions##
## TDLMM:
##
## Model run info
## - 20 trees
## - 1000 burn-in iterations
## - 2000 post-burn iterations
## - 2 thinning factor
## - 5 exposures
## - 0.95 confidence level
##
## Fixed effects:
## Mean Lower Upper
## (Intercept) -14.5044 -24.5646 0.9092
## *c1 -1.1408 -1.2812 -0.9967
## *c2 -1.5724 -1.7127 -1.4314
## c3 -0.0097 -0.1544 0.1180
## c4 -0.0757 -0.2265 0.0635
## *c5 1.3087 1.1755 1.4450
## b1 0.1842 -0.0720 0.4516
## b2 -0.2000 -0.4977 0.1004
## *b3 0.8427 0.5711 1.1216
## *b4 1.0605 0.7949 1.3217
## *b5 0.6183 0.3516 0.9093
## ---
## * = CI does not contain zero
##
## --
## Exposure effects: critical windows
##
## *e1 (0.92): 19-27
## *e2 (0.77): 23-27
## e3 (0.17):
## e4 (0.25):
## e5 (0.39):
## ---
## * = Model selected exposure or interaction
## (0.xx) = Relative signal strengthPlot the DLM using plot with the model summary.
plot(m5_sum, exposure1 = "e1",
main = "Estimated DLM for active exposure 1", xlab = "Time", ylab = "Effect")
plot(m5_sum, exposure1 = "e2",
main = "Estimated DLM for a active exposure 2", xlab = "Time", ylab = "Effect")
Compare to true generated DLM.
# RMSE exposure 1
sqrt(mean((D$margDLM1 - m5_sum$DLM$e1$marg.matfit)^2))## [1] 0.03684105# RMSE exposure 2
sqrt(mean((D$margDLM2 - m5_sum$DLM$e2$marg.matfit)^2))## [1] 0.02360344# Coverage exposure 1
mean((m5_sum$DLM$e1$marg.cilower < D$margDLM1) &
(m5_sum$DLM$e1$marg.ciupper > D$margDLM1))## [1] 0.972973# Coverage exposure 2
mean((m5_sum$DLM$e2$marg.cilower < D$margDLM2) &
(m5_sum$DLM$e2$marg.ciupper > D$margDLM2))## [1] 1# True positive effect classification exposure 1
sum((m5_sum$DLM$e1$marg.cilower > 0) & (D$margDLM1 > 0)) /
sum(D$margDLM1 > 0)## [1] 1# True positive effect classification exposure 2
sum((m5_sum$DLM$e2$marg.cilower > 0) & (D$margDLM2 > 0)) /
sum(D$margDLM2 > 0)## [1] 0.625# False positive effect classification exposure 1
sum(((m5_sum$DLM$e1$marg.cilower > 0) |
(m5_sum$DLM$e1$marg.ciupper < 0)) & (D$margDLM1 == 0)) /
sum(D$margDLM1 == 0)## [1] 0.03448276# False positive effect classification exposure 2
sum(((m5_sum$DLM$e2$marg.cilower > 0) |
(m5_sum$DLM$e2$marg.ciupper < 0)) & (D$margDLM2 == 0)) /
sum(D$margDLM2 == 0)## [1] 0# Exposure selection
m5_sum$expSel## e1 e2 e3 e4 e5
## TRUE TRUE FALSE FALSE FALSE# Exposure relative effect size
m5_sum$expVar[2,]## e1 e2 e3 e4 e5
## 0.91525 0.77225 0.16700 0.25450 0.39100