Overview

This vignette details the posthoc TDLMM analysis to adjust for expected changes in co-occuring exposures, as described in Section 6.3 of the corresponding manuscript. In particular, we calculate

\begin{matrix} E [Y | {\tilde{X}}_{t} = E {X_{t} | x_{m t} = x_{m (0.75)}}, {\tilde{X}}_{[t]} = \bar{x}, z = z_{0}] \\ - E [Y | {\tilde{X}}_{t} = E {X_{t} | x_{m t} = x_{m (0.25)}}, {\tilde{X}}_{[t]} = \bar{x}, z = z_{0}] . \end{matrix}

$\begin{multline}\label{eq:exp_dlm} \mbox{E}\left[Y \Big| \widetilde{\mathbf{X}}_t = \mbox{E}\left\{\mathbf{X}_t \big| x_{mt}=x_{m(0.75)}\right\}, \widetilde{\mathbf{X}}_{[t]}=\overline{\mathbf{x}}, \mathbf{z}=\mathbf{z}_0\right]\\ -\mbox{E}\left[Y \Big| \widetilde{\mathbf{X}}_t = \mbox{E}\left\{\mathbf{X}_t \big| x_{mt}=x_{m(0.25)}\right\}, \widetilde{\mathbf{X}}_{[t]}=\overline{\mathbf{x}}, \mathbf{z}=\mathbf{z}_0\right]. \end{multline}$ where

{\tilde{X}}_{t} = {{\tilde{x}}_{1 t}, \dots, {\tilde{x}}_{M t}}

$\widetilde{\mathbf{X}}_t=\{\widetilde{x}_{1t},\ldots,\widetilde{x}_{Mt}\}$ defines the values of all exposures at time

t

$t$ while

{\tilde{X}}_{[t]}

$\widetilde{\mathbf{X}}_{[t]}$ is the collection of exposure measurements at all time points except

t

$t$ . We set each element of

{\tilde{X}}_{[t]}

$\widetilde{\mathbf{X}}_{[t]}$ equal to the empirical mean for that exposure (

\bar{x}

$\overline{\mathbf{x}}$ ) to isolate the effect only at the time of interest. The value of

z_{0}

$\mathbf{z}_0$ , which we set to a population baseline level, does not influence the expected change because it does not interact with exposure measurements.

Data

The exposure data is taken directly from our analysis, where each row is a single observation day (37 days per individual) and columns are the different exposure types. Each exposure has been scaled by its IQR. We gather the 25 $^{th}$ and 75 $^{th}$ percentiles of each exposure for predictions.

# load("path/of/exposureData")
head(exposureData)                # Data has been scaled by IQR

ABCDEFGHIJ0123456789

	log_pm <dbl>	no2 <dbl>	so2 <dbl>	co <dbl>	temp <dbl>
1	6.299177	3.597932	5.091201	3.978367	1.077205
2	6.449174	3.579307	2.852399	4.936389	1.094033
3	6.064634	3.587742	3.866215	4.712110	1.111355
4	6.895106	3.580315	2.973472	4.815050	1.094453
5	5.670315	3.590588	4.208343	4.612842	1.170800
6	5.664042	3.591482	4.315801	4.480957	1.096927

expMean <- colMeans(exposureData) # Empirical means
expIQRLevels <-                   # 25/75 percentiles of each exposure
  lapply(exposureData, function(i) quantile(i, probs = c(0.25, 0.75)))
do.call(cbind.data.frame, expIQRLevels)

ABCDEFGHIJ0123456789

	log_pm <dbl>	no2 <dbl>	so2 <dbl>	co <dbl>	temp <dbl>
25%	5.184863	1.872663	0.9961919	2.094969	1.377652
75%	6.184863	2.872663	1.9961919	3.094969	2.377652

Predicted levels of co-occuring exposures

To estimate the values of co-exposures at time $t$ , given by $\mbox{E}[\mathbf{X}_t \big|x_{mt}=x_{m(q)}]$ , we fit penalized spline models (cubic splines with 5 degrees of freedom) between pairs of exposures, estimated using restricted maximum likelihood using R package mgcv. Specifically, we considered measurements for exposure $m$ as the only predictor and fit separate models for each co-exposure (e.g. log PM $_{2.5}$ was used as the predictor for NO $_2$ , SO $_2$ , CO, and temperature in four separate models). This model matches our assumption of a simple, smooth relationship between pairs of exposures. Other models could also be substituted for the penalized spline model. Using the model fits, we predicted all co-exposures at the 25 $^{th}$ and 75 $^{th}$ percentiles of exposure $m$ .

library(mgcv)

## Loading required package: nlme

## This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.

# predLevels[[predictor]][[response]] =
#  use level of 'predictor' to estimate 'response'.
predLevels <- lapply(exposureData, function(i) list())
# Loop over pairs of exposures
for (predExp in names(predLevels)) {
  for (respExp in names(predLevels)) {
    cat("\nPredicting", respExp, "at 25/75 percentiles of", predExp)
    if (predExp == respExp) { 
      # If predictor/response exposure same, set to 25/75 percentiles.
      predLevels[[predExp]][[respExp]] <- expIQRLevels[[predExp]]
      
    } else { 
      # Predict `response` (respExp) at 25/75 percentiles of `predictor` (predExp).
      # Other methods could be substituted here to achieve prediction
      formula <- as.formula(paste(respExp, "~s(", predExp, ", k = 5, bs = 'cr')"))
      model <- bam(formula, dat = exposureData)
      newdat <- data.frame(expIQRLevels[[predExp]])
      names(newdat) <- predExp
      predLevels[[predExp]][[respExp]] <- predict(model, newdata = newdat)
    }
    cat(":", predLevels[[predExp]][[respExp]])
  }
}

## 
## Predicting log_pm at 25/75 percentiles of log_pm: 5.184863 6.184863
## Predicting no2 at 25/75 percentiles of log_pm: 2.219529 2.442719
## Predicting so2 at 25/75 percentiles of log_pm: 1.409915 1.701812
## Predicting co at 25/75 percentiles of log_pm: 2.489573 2.830155
## Predicting temp at 25/75 percentiles of log_pm: 1.874715 1.904258
## Predicting log_pm at 25/75 percentiles of no2: 5.527204 5.752684
## Predicting no2 at 25/75 percentiles of no2: 1.872663 2.872663
## Predicting so2 at 25/75 percentiles of no2: 1.296038 1.698897
## Predicting co at 25/75 percentiles of no2: 2.243837 3.001701
## Predicting temp at 25/75 percentiles of no2: 2.105268 1.592558
## Predicting log_pm at 25/75 percentiles of so2: 5.584255 5.877585
## Predicting no2 at 25/75 percentiles of so2: 2.183866 2.540504
## Predicting so2 at 25/75 percentiles of so2: 0.9961919 1.996192
## Predicting co at 25/75 percentiles of so2: 2.439064 2.908714
## Predicting temp at 25/75 percentiles of so2: 1.91681 1.831812
## Predicting log_pm at 25/75 percentiles of co: 5.490305 5.865365
## Predicting no2 at 25/75 percentiles of co: 1.973449 2.686335
## Predicting so2 at 25/75 percentiles of co: 1.315538 1.646576
## Predicting co at 25/75 percentiles of co: 2.094969 3.094969
## Predicting temp at 25/75 percentiles of co: 2.173791 1.607157
## Predicting log_pm at 25/75 percentiles of temp: 5.692187 5.728444
## Predicting no2 at 25/75 percentiles of temp: 2.739964 1.947469
## Predicting so2 at 25/75 percentiles of temp: 1.700598 1.357958
## Predicting co at 25/75 percentiles of temp: 3.122098 2.342968
## Predicting temp at 25/75 percentiles of temp: 1.377652 2.377652

Adjusting for changes in co-exposures

Load the output from tdlmm() and run summary() with option keep.mcmc = TRUE to extract the MCMC samples taken from the model posterior distribution.
For each exposure at each time point, we add the following:
1. Difference of MCMC main effects from the upper and lower quartiles of concentration for the exposure of interest
2. Difference of MCMC main effects for all co-exposures at the predicted concentration corresponding to the upper and lower quartiles of concentration for exposure of interest
3. Difference of same-time MCMC interaction effects taken at upper and lower quartiles of concentration in exposure of interest and corresponding predicted concentrations in co-exposures
4. Difference of same-time MCMC interaction effects taken at corresponding predicted concentration in pairs of co-exposures
5. Difference of MCMC interaction effects at time of interest in exposure of interest using upper and lower quartiles for exposure of interest combined with other time interactions for co-exposures set at their empirical means
6. Difference of MCMC interaction effects across other time interactions for co-exposures set at their empirical means (which equals zero)
Summarize MCMC results into mean effect at each time point, upper and lower credible values, and critical windows where credible interval does not contain zero.

library(dlmtree)
# load("path/to/tdlmmResult")
# Use summary function to get all MCMC samples from TDLMM
modelSummary <- summary(tdlmmResult, keep.mcmc = TRUE, verbose = FALSE)
expDLM <- list() # For output
for (exposure in names(modelSummary$DLM)) {
  # MCMC results for main effect of interest
  mcmc <- modelSummary$DLM[[exposure]]$mcmc * diff(expIQRLevels[[exposure]])
  # Add MCMC results for expected changes in co-exposures
  for (coexposure in names(modelSummary$DLM)) {
    if (exposure != coexposure)
      mcmc <- mcmc + 
        modelSummary$DLM[[coexposure]]$mcmc * diff(predLevels[[exposure]][[coexposure]])
  }
  # Add MCMC results due to interactions: 
  # 1) same time interactions, use expected changes in main and co-exposure
  # 2) diff time interactions, use expected change in main, mean in co-exposure
  # 3) same time interactions (b/t co-exp), use expected change in co-exposures
  # 4) diff time interactions (b/t co-exp), cancels (set at mean)
  for (mix in 1:length(modelSummary$MIX)) {
    # Interaction between main and co-exposure
    if (modelSummary$MIX[[mix]]$rows == exposure) {
      mcmc <- mcmc +
        # add same time interactions using expected changes
        t(sapply(1:37, function(k) {
          modelSummary$MIX[[mix]]$mcmc[k,k,]
        })) * 
        (predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][2] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][2] -
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][1] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][1]) +
        # add different time interactions using IQR change in primary exposure
        # (other exposures set to mean)
        t(sapply(1:37, function(k) {
          colSums(modelSummary$MIX[[mix]]$mcmc[k,-k,]) # interaction sum removing time k
        })) * 
        diff(predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]]) *
        expMean[modelSummary$MIX[[mix]]$cols]
    # Interaction between main and co-exposure
    } else if (modelSummary$MIX[[mix]]$cols == exposure) {
      mcmc <- mcmc +
        # add same time interactions using expected changes
        t(sapply(1:37, function(k) {
          modelSummary$MIX[[mix]]$mcmc[k,k,]
        })) * 
        (predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][2] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][2] -
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][1] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][1]) +
        # add different time interactions using IQR change in primary exposure
        # (other exposure set to mean)
        t(sapply(1:37, function(k) {
          colSums(modelSummary$MIX[[mix]]$mcmc[-k,k,]) # interaction sum removing time k
        })) * 
        diff(predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]]) *
        expMean[modelSummary$MIX[[mix]]$rows]
    # Interaction between two co-exposures
    } else {
      mcmc <- mcmc +
        # add same time interactions using expected changes
        t(sapply(1:37, function(k) {
          modelSummary$MIX[[mix]]$mcmc[k,k,]
        })) * 
        (predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][2] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][2] -
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$rows]][1] * 
         predLevels[[exposure]][[modelSummary$MIX[[mix]]$cols]][1])
      # no different time interactions (they are set at mean, so equals zero)
    }
  }
  # Summarize main effects, adjusted for changes in co-exposures
  expDLM[[exposure]] <- 
    data.frame("Name" = exposure,
               "Week" = 1:modelSummary$nLags,
               "Effect" = apply(mcmc, 1, mean),
               # Credible interval at levels 0.025 and 0.975
               "Lower" =  apply(mcmc, 1, quantile, probs = .025),
               "Upper" =  apply(mcmc, 1, quantile, probs = .975))
  expDLM[[exposure]]$CW <- # Critical window where CI does not contain zero
    (expDLM[[exposure]]$Lower > 0 | expDLM[[exposure]]$Upper < 0)
}
rm(modelSummary, tdlmmResult)

Plot of results

library(ggplot2)
expDLM <- do.call(rbind.data.frame, expDLM)
ggplot(expDLM, aes(x = Week, y = Effect, ymin = Lower, ymax = Upper)) +
  geom_hline(yintercept = 0, color = "red", size = 1) +
  geom_ribbon(fill = "grey", col = "grey", size = 1) +
  geom_line(size = 1) +
  facet_wrap(~Name, nrow = 2) +
  theme_bw(base_size = 12) +
  theme(strip.background = element_rect(fill = "#DDDDDD")) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  ylab("Change in outcome for IQR change in exposure \nadjusting for expected changes in co-exposures") +
  xlab("Week of gestation") +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())

TDLMM Vignette: Adjust for Changes in Co-Exposures

Daniel Mork

6/30/2021

Overview

Data

Predicted levels of co-occuring exposures

Adjusting for changes in co-exposures

Plot of results