This vignette demonstrates the implementation of treed distributed lag mixture model (TDLMM). More details can be found in Mork and Wilson (2023) <doi: 10.1111/biom.13568>.
Load data
Simulated data is available on GitHub. It can be loaded with the following code.
sbd_dlmtree <- get_sbd_dlmtree()Data preparation
# Response and covariates
sbd_cov <- sbd_dlmtree %>%
select(bwgaz, ChildSex, MomAge, GestAge, MomPriorBMI, Race,
Hispanic, MomEdu, SmkAny, Marital, Income,
EstDateConcept, EstMonthConcept, EstYearConcept)
# Exposure data
sbd_exp <- list(PM25 = sbd_dlmtree %>% select(starts_with("pm25_")),
TEMP = sbd_dlmtree %>% select(starts_with("temp_")),
SO2 = sbd_dlmtree %>% select(starts_with("so2_")),
CO = sbd_dlmtree %>% select(starts_with("co_")),
NO2 = sbd_dlmtree %>% select(starts_with("no2_")))
sbd_exp <- sbd_exp %>% lapply(as.matrix)Fitting the model
tdlmm.fit <- dlmtree(formula = bwgaz ~ ChildSex + MomAge + MomPriorBMI +
Race + Hispanic + SmkAny + EstMonthConcept,
data = sbd_cov,
exposure.data = sbd_exp,
family = "gaussian", dlm.type = "linear", mixture = TRUE,
mixture.interactions = "noself",
n.burn = 2500, n.iter = 10000, n.thin = 5)
#> Preparing data...
#>
#> Running TDLMM:
#> Burn-in % complete
#> [0--------25--------50--------75--------100]
#> ''''''''''''''''''''''''''''''''''''''''''
#> MCMC iterations (est time: 2.9 minutes)
#> [0--------25--------50--------75--------100]
#> ''''''''''''''''''''''''''''''''''''''''''
#> Compiling results...Model fit summary
# Marginalization with co-exposure fixed at exact levels for each exposure
tdlmm.sum.exact <- summary(tdlmm.fit, marginalize = c(3, 2, 1, 2, 1))
#> Specified co-exposure values:
#> - PM25 : 3
#> - TEMP : 2
#> - SO2 : 1
#> - CO : 2
#> - NO2 : 1
#>
#> Reconstructing main effects...
#> Reconstructing interaction effects...
#> 0%...25%...50%...75%...100%
#> Calculating marginal effects...
#> Calculating fixed effects...
# Marginalization with co-exposure fixed at 25th percentile
tdlmm.sum.percentile <- summary(tdlmm.fit, marginalize = 25)
#> Specified co-exposure values:
#> - PM25 : 5.183927
#> - TEMP : 1.374847
#> - SO2 : 0.9961942
#> - CO : 2.09893
#> - NO2 : 1.872199
#>
#> Reconstructing main effects...
#> Reconstructing interaction effects...
#> 0%...25%...50%...75%...100%
#> Calculating marginal effects...
#> Calculating fixed effects...
# Marginalization with co-exposure fixed at the empirical means (default)
tdlmm.sum <- summary(tdlmm.fit, marginalize = "mean", log10BF.crit = 0.5)
#> Specified co-exposure values:
#> - PM25 : 6.842084
#> - TEMP : 3.187414
#> - SO2 : 1.992197
#> - CO : 3.219101
#> - NO2 : 3.077839
#>
#> Reconstructing main effects...
#> Reconstructing interaction effects...
#> 0%...25%...50%...75%...100%
#> Calculating marginal effects...
#> Calculating fixed effects...
tdlmm.sum
#> ---
#> TDLMM summary
#>
#> Model run info:
#> - bwgaz ~ ChildSex + MomAge + MomPriorBMI + Race + Hispanic + SmkAny + EstMonthConcept
#> - sample size: 10,000
#> - family: gaussian
#> - 20 trees (alpha = 0.95, beta = 2)
#> - 2500 burn-in iterations
#> - 10000 post-burn iterations
#> - 5 thinning factor
#> - 5 exposures measured at 37 time points
#> - 10 two-way interactions (no-self interactions)
#> - 1 kappa sparsity prior
#> - 0.95 confidence level
#>
#> Fixed effects:
#> Mean Lower Upper
#> *(Intercept) 0.172 0.044 0.299
#> *ChildSexM -2.063 -2.085 -2.041
#> MomAge 0.000 -0.001 0.002
#> *MomPriorBMI -0.020 -0.022 -0.019
#> RaceAsianPI 0.028 -0.066 0.115
#> RaceBlack 0.033 -0.063 0.126
#> Racewhite 0.017 -0.070 0.095
#> *HispanicNonHispanic 0.248 0.224 0.272
#> *SmkAnyY -0.395 -0.442 -0.348
#> *EstMonthConcept2 0.082 0.011 0.154
#> *EstMonthConcept3 0.123 0.033 0.219
#> *EstMonthConcept4 0.156 0.039 0.275
#> *EstMonthConcept5 0.233 0.104 0.363
#> *EstMonthConcept6 0.175 0.039 0.313
#> *EstMonthConcept7 0.211 0.080 0.343
#> *EstMonthConcept8 0.210 0.087 0.327
#> *EstMonthConcept9 0.303 0.187 0.412
#> *EstMonthConcept10 0.185 0.073 0.288
#> *EstMonthConcept11 0.120 0.027 0.211
#> EstMonthConcept12 -0.010 -0.075 0.055
#> ---
#> * = CI does not contain zero
#>
#> --
#> Exposure effects: critical windows
#> * = Exposure selected by Bayes Factor
#> (x.xx) = Relative effect size
#>
#> *PM25 (0.81): 4-5,11-20
#> *TEMP (0.6): 4-20
#> *NO2 (0.84): 8-19,23-25
#> --
#> Interaction effects: critical windows
#>
#> PM25/TEMP (0.93):
#> 12/6-16
#> 13/6-16
#> 14/6-16
#> 15/6-16
#> 16/6-12,14-16
#> 17/5-18
#> 18/5-19
#> 19/5-18
#> PM25/NO2 (0.95):
#> 11/11-21,23
#> 12/8-25
#> 13/7-25
#> 14/7-24
#> 15/6-25
#> 16/6-24
#> 17/7-24
#> 18/7-24
#> 19/7-24
#> 20/7-12,18
#> ---
#> residual standard errors: 0.007Main exposure effect
p1 <- plot(tdlmm.sum, exposure1 = "PM25", main = "PM2.5")
p2 <- plot(tdlmm.sum, exposure1 = "TEMP", main = "Temperature")
p3 <- plot(tdlmm.sum, exposure1 = "SO2", main = "SO2")
p1
p2
p3
Adjusting for expected changes in co-occurring exposures
Here we consider going from the 25th to the 75th percentile in each exposure while adjusting for the expected changes in other exposures due to their correlations with the exposure of interest.
library(ggplot2)
dlm_coexp <- adj_coexposure(sbd_exp, tdlmm.fit, contrast_perc = c(0.25, 0.75))
#> PM25 TEMP SO2 CO NO2
#> 25% 5.183927 1.374847 0.9961942 2.098930 1.872199
#> 75% 6.187326 2.378621 1.9960009 3.096027 2.875369
#>
#> Predicting PM25 at lower/upper values of PM25: 5.183927 6.187326
#> Predicting TEMP at lower/upper values of PM25: 1.886921 1.937889
#> Predicting SO2 at lower/upper values of PM25: 1.403909 1.672338
#> Predicting CO at lower/upper values of PM25: 2.477085 2.784547
#> Predicting NO2 at lower/upper values of PM25: 2.208032 2.409519
#> Predicting PM25 at lower/upper values of TEMP: 5.654899 5.701115
#> Predicting TEMP at lower/upper values of TEMP: 1.374847 2.378621
#> Predicting SO2 at lower/upper values of TEMP: 1.686905 1.353916
#> Predicting CO at lower/upper values of TEMP: 3.07044 2.314228
#> Predicting NO2 at lower/upper values of TEMP: 2.72529 1.952642
#> Predicting PM25 at lower/upper values of SO2: 5.589447 5.87508
#> Predicting TEMP at lower/upper values of SO2: 1.888874 1.841457
#> Predicting SO2 at lower/upper values of SO2: 0.9961942 1.996001
#> Predicting CO at lower/upper values of SO2: 2.446352 2.898562
#> Predicting NO2 at lower/upper values of SO2: 2.179354 2.554288
#> Predicting PM25 at lower/upper values of CO: 5.482744 5.85876
#> Predicting TEMP at lower/upper values of CO: 2.0984 1.659518
#> Predicting SO2 at lower/upper values of CO: 1.310125 1.617886
#> Predicting CO at lower/upper values of CO: 2.09893 3.096027
#> Predicting NO2 at lower/upper values of CO: 1.971065 2.673015
#> Predicting PM25 at lower/upper values of NO2: 5.529577 5.732295
#> Predicting TEMP at lower/upper values of NO2: 2.121927 1.660641
#> Predicting SO2 at lower/upper values of NO2: 1.297695 1.702796
#> Predicting CO at lower/upper values of NO2: 2.25247 2.968066
#> Predicting NO2 at lower/upper values of NO2: 1.872199 2.875369
ggplot(dlm_coexp, aes(x = Time, y = Effect, ymin = Lower, ymax = Upper)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line() +
facet_wrap(~Name) +
theme_bw() +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0))