bwqs_sim1.R

# simulation 1 from appendix of 
#  ”Per- and poly-fluoroalkyl substances and bone mineral density: 
#    Results from the Bayesian weighted quantile sum regression” by Colicino et al.
#
#  simulate data and run the BWQS model

library(rstan)
library(qgcomp)
parallel:::setDefaultClusterOptions(setup_strategy = "sequential") # kludge for bug in parallel package with r4.0
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

dgm <- function(N=100, sigma=1.0, psi=0.8, w = c(0.5, 0.2, 0.2, 0.05, 0.05)){
  SigmaX = matrix(
    c(1.0000,0.6500,0.4225,1.2675,0.4225,
      0.6500,1.0000,0.4225,1.2675,0.4225,
      0.4225,0.4225,1.0000,1.2675,0.4225,
      1.2675,1.2675,1.2675,9.0000,1.9500,
      0.4225,0.4225,0.4225,1.9500,1.0000)
    , nrow=5, byrow=TRUE)
  muX = c(0,0,0,1,3)
  X = qgcomp:::.rmvnorm(N, muX, SigmaX)
  mu = X %*% cbind(psi*w)
  y = mu + rnorm(N, 0, sigma)
  list(y=as.numeric(y),X=X, N=N, p=5, dprior=rep(1, 5))
}


stanmod <- "
data{
  int<lower=0> N;
  int<lower=0> p;
  vector[p] dprior;
  matrix[N,p] X;
  vector[N] y;
}
parameters{
 real psi;
 simplex[p] w;
 real b0;
 real<lower=0> sigma;
}
model{
  vector[N] wqs = X*w;
  w ~ dirichlet(dprior);
  y ~ normal(b0 + psi * wqs, sigma);
}
"

sim1 = dgm(1000, sig=1, psi=0.8)

# check that code works
testfit = stan(model_code = stanmod, data=sim1, chains=1, iter=1)

# sample from model
ft = stan(fit=testfit, data=sim1, chains=4, iter=1000)
ft

Note that this simulation is not of the full BWQS approach, which would first quantize the exposure variables of interest.

Reference: Colicino E, Pedretti NF, Busgang SA, et al. Per- and poly-fluoroalkyl substances and bone mineral density: Results from the Bayesian weighted quantile sum regression. Environmental Epidemiology. 2020;4(3):e092.