NOTE: If you are interested in the effects of longitudinally assessed exposure to a mixture please see the Treed Distributed Lag Mixture Model (TDLMM) method in the package dlmtree.
First, install and load the regimes package.
We simulate data to illustrate the required data structures and the method.
The exposure data is a list. Each element of the list is a matrix of values for a single exposure. The columns of each matrix are are the ordered measurements at each time point starting with the first time point in the leftmost column. The rows are subjects. Each element of the list must have the same dimensions (same number of individuals and time points).
## [1] "list"
## [1] 2
## $x1
## [1] 200 37
##
## $x2
## [1] 200 37
The response variable is a vector of length n. This should be continuous. There is also a matrix of covariates that should not include an intercept.
## [1] 200 5
## [1] 200
The main function for estimation is bkmrdlm
. We fit
BKMR-DLM with a quadratic kernel. Use gaussian=TRUE
for a
Gaussian kernel instead of a polynomial kernel. When using BKMR-DLM for
real data analysis increase the number of iterations. This may take a
few minutes or longer depending the the data size and number of
iterations.
Now display numerical results using the summary
method.
This reports any critical windows identified and reports estimates of
the regression coefficients for the covariates.
## [1] "call" "windows" "weights" "sdres" "coefficients"
##
## Call:
##
## bkmrdlm(y = dat$y, x = dat$x, z = dat$z, niter = 2000, gaussian = FALSE,
## polydegree = 2)
##
##
## Significant windows identified:
##
## x1 11-28, 34-37
## x2 11-37
##
##
## Coefficients for covariates:
##
## Post. Mean Post. SD q2.5 q97.5 Pr>0
## z1 0.3082875 0.07793993 0.1565000 0.4601251 1
## z2 0.8848398 0.06504687 0.7520208 1.0050850 1
## z3 1.2097626 0.08163969 1.0600426 1.3738288 1
## z4 -2.7809500 0.07964868 -2.9421545 -2.6229652 0
## z5 -0.3588762 0.08295008 -0.5348712 -0.2036607 0
##
##
## Posterior mean residuals standard deviation: 1.022655
The plot
method displays the estimated weight function
for each exposure. Times where the displayed 95 percent intervals do not
contain zero are considered to be critical windows. The sign of the
weight function does not indicate a direction of the association with
the outcome. The direction of association is determined by the
exposure-response function.
We first illustrate the exposure-response function at the median
value of all other exposures. This is one form of the main effect. The
first step is to predict the exposure response function at a set of
values. This may take some time. The number of points
is
the number of points that the function is predicted at for each
exposure.
After prediction we can visualise the exposure-response function.
library(ggplot2)
p <- ggplot(univarpred_quad, aes(x=E, y=mean, ymin=lower, ymax=upper))
p <- p + geom_ribbon(alpha=.5)
p <- p + geom_line(size=2) + facet_wrap(~reorder(name,m), scales="free_x")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
To visualize interactions we can plot the exposure-response function for each exposure at different quantiles of the other exposures. First we estimate the exposure-response function at five quantiles of the other exposures.
# estimate exposure response
combodat <- NULL
for(q in c(0.10,0.25,0.75,0.90)){
for(m in 1:length(dat$x)){
combodat <- rbind(combodat,predictUnivariate(fit, points=30,crossM=m, qtl=q))
}
}
# some formatting
combodat <- combodat[-which(combodat$m==combodat$cross_m),]
univarpred_quad2 <- univarpred_quad
univarpred_quad2$cross <- univarpred_quad2$name
univarpred_quad2$cross_m <- univarpred_quad2$m
Next we can plot the results.