Created
August 13, 2020 16:15
-
-
Save alexpkeil1/e527e6aa6f2a26171d52049d65505644 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.