Fork me on GitHub

Update on new package: bdlim

If you are interested in the BDLIM method please see the new package with faster implementation, expanded functionality, and improved user interface: https://anderwilson.github.io/bdlim/

Bayesian Distributed Lag Interaction Models

Ander Wilson

An Example

Load the regimes package, simulate data, and fit BDLIM with all four options (BDLIM-n, BDLIM-b, BSLIM-w, BDLIM-bw). When using BDLIM for data analysis increase the number of iterations.

#load package
library(regimes)
#simulate data from scenario 3
dat <- simregimes(scenario="bdlim2", seed=1234, n=1000)

#estimate model
fit <- bdlim(dat$Y,dat$X,dat$Z,dat$G,"all",niter=1000)

Now display some results using the summary method.

#summary
sfit <- summary(fit)
names(sfit)
## [1] "beta"         "cumulative"   "w"            "bw"           "model"       
## [6] "modelfit"     "call"         "windows"      "coefficients"
sfit
## 
## Model fit statistics:
##               DIC     pD modelprob
## BDLIM_bw 6421.401 23.798     0.937
## BDLIM_b  6419.175 18.870     0.063
## BDLIM_w  6633.567 25.610     0.000
## BDLIM_n  7070.935 19.213     0.000
## 
## 
## Posterior results for BDLIM_bw:
## 
## Beta:
##         mean          sd        q2.5      q97.5 Pr>0    n_eff
## 0  0.1070450 0.008117157  0.09233272  0.1248538    1 18.81189
## 1 -0.2042449 0.006049035 -0.21663814 -0.1929648    0 34.74529
## 
## Cumulative:
##        mean        sd      q2.5     q97.5 Pr>0    n_eff
## 0  2.969692 0.2079919  2.623210  3.389669    1 5.198216
## 1 -6.105803 0.1698880 -6.460251 -5.785855    0 6.970461
## 
## Critical windows identified with weighted exposures, beta*w(t):
##              
## 0 5-31, 36-37
## 1        3-30
## 
## n_eff for beta*w(t): min 25.3, max 500, mean 328.6, median 353.3
## n_eff for w(t): min 123.6, max 1135.8, mean 493.8, median 500
## 
## Coefficients for covariates:
## 
##            mean        sd        q2.5      q97.5  Pr>0      n_eff
## G0   0.95944372 2.0959020 -3.27549204  4.3784869 0.696   5.082563
## G1   0.77557420 1.7151645 -2.39977903  4.3182383 0.708   6.318867
## Z1  -0.05247123 0.2005363 -0.46360920  0.3165683 0.402 500.000000
## Z2  -0.15036258 0.1862671 -0.51480902  0.2037550 0.224 500.000000
## Z3  -0.24169066 0.1834509 -0.61443342  0.1100082 0.092 269.906974
## Z4  -0.78938477 0.1811849 -1.12352644 -0.4309955 0.000 500.000000
## Z5   0.33218046 0.1953857 -0.05976555  0.7028209 0.958 464.298975
## Z6   0.92861930 0.1812740  0.59457956  1.2519625 1.000 500.000000
## Z7  -0.60418044 0.1908247 -0.95433441 -0.2025071 0.002 465.839120
## Z8   1.34554671 0.1986730  0.96716318  1.7377969 1.000 327.566346
## Z9   0.02866055 0.1951945 -0.36374748  0.4162282 0.552 589.279249
## Z10  0.85203406 0.1902137  0.47862619  1.1982570 1.000 432.083096

Now display some results via plots.

plot(sfit, print=TRUE)