bhpm was developed for the Precision Drug Theraputics: Risk Prediction in Pharmacoepidemiology project as part of a Rutherford Fund Fellowship at Health Data Research (UK), Medical Research Council (UK) award reference MR/S003967/1 (https://gtr.ukri.org/).
The goal of bhpm is to investigate associations between multiple outcomes and corresponding patient treatments. bhpm implements Bayesian hierarchical models, which allow a stratification of the population into clusters with similar characteristics, and which take advantage of known relationships between clinical outcomes, to determine which outcomes are associated with treatments.
You can install the released version of bhpm from CRAN with:
install.packages("bhpm")And the development version from GitHub with:
# install.packages("devtools")
devtools::install_github("rcarragh/bhpm")This is a basic example which shows how to fit the model:
library(bhpm)
#> Loading required package: coda
#> Warning: package 'coda' was built under R version 4.4.3
data(demo.cluster.data)
mod.fit <- bhpm.pm(demo.cluster.data)#, burnin = 100, iter = 200)
#> Memory Model: HIGH
#> 1000 iterations...
#> 2000 iterations...
#> 3000 iterations...
#> 4000 iterations...
#> 5000 iterations...
#> 6000 iterations...
#> 7000 iterations...
#> 8000 iterations...
#> 9000 iterations...
#> 10000 iterations...
#> 11000 iterations...
#> 12000 iterations...
#> 13000 iterations...
#> 14000 iterations...
#> 15000 iterations...
#> 16000 iterations...
#> 17000 iterations...
#> 18000 iterations...
#> 19000 iterations...
#> 20000 iterations...
#> 21000 iterations...
#> 22000 iterations...
#> 23000 iterations...
#> 24000 iterations...
#> 25000 iterations...
#> 26000 iterations...
#> 27000 iterations...
#> 28000 iterations...
#> 29000 iterations...
#> 30000 iterations...
#> 31000 iterations...
#> 32000 iterations...
#> 33000 iterations...
#> 34000 iterations...
#> 35000 iterations...
#> 36000 iterations...
#> 37000 iterations...
#> 38000 iterations...
#> 39000 iterations...
#> 40000 iterations...
#> 41000 iterations...
#> 42000 iterations...
#> 43000 iterations...
#> 44000 iterations...
#> 45000 iterations...
#> 46000 iterations...
#> 47000 iterations...
#> 48000 iterations...
#> 49000 iterations...
#> 50000 iterations...
#> 51000 iterations...
#> 52000 iterations...
#> 53000 iterations...
#> 54000 iterations...
#> 55000 iterations...
#> 56000 iterations...
#> 57000 iterations...
#> 58000 iterations...
#> 59000 iterations...
#> 60000 iterations...
#> MCMC fitting complete.mod.conv <- bhpm.convergence.diag(mod.fit)
bhpm.print.convergence.summary(mod.conv)
#> Summary Convergence Diagnostics:
#> ================================
#> theta:
#> ------
#> Max Gelman-Rubin diagnostic (4 M/65-84 Bleed Bleed_Outcome4): 1.008329
#> Min Gelman-Rubin diagnostic (2 M/85+, G00-G99 G00-99_Outcome3): 1.000007
#> gamma:
#> ------
#> Max Gelman-Rubin diagnostic (M/85+ G00-G99 G00-99_Outcome2): 1.000626
#> Min Gelman-Rubin diagnostic (M/65-84 Bleed Bleed_Outcome4): 0.999987
#> mu.gamma:
#> ---------
#> Max Gelman-Rubin diagnostic (Bleed): 1.000026
#> Min Gelman-Rubin diagnostic (I00-I99): 1.000005
#> mu.theta:
#> ---------
#> Max Gelman-Rubin diagnostic (3 I00-I99): 1.011972
#> Min Gelman-Rubin diagnostic (2 G00-G99): 1.000005
#> sigma2.gamma:
#> -------------
#> Max Gelman-Rubin diagnostic (I00-I99): 1.000121
#> Min Gelman-Rubin diagnostic (Bleed): 1.000097
#> sigma2.theta:
#> -------------
#> Max Gelman-Rubin diagnostic (3 I00-I99): 1.004777
#> Min Gelman-Rubin diagnostic (4 I00-I99): 1.000110
#> pi:
#> ---
#> Max Gelman-Rubin diagnostic (3 G00-G99): 1.001848
#> Min Gelman-Rubin diagnostic (4 G00-G99): 1.000040
#> mu.gamma.0:
#> -----------
#> Max Gelman-Rubin diagnostic: 1.000064
#> Min Gelman-Rubin diagnostic: 1.000064
#> mu.theta.0:
#> -----------
#> Max Gelman-Rubin diagnostic: 3 1.011682
#> Min Gelman-Rubin diagnostic: 2 1.000093
#> tau2.gamma.0:
#> -------------
#> Max Gelman-Rubin diagnostic: 1.000264
#> Min Gelman-Rubin diagnostic: 1.000264
#> tau2.theta.0:
#> -------------
#> Max Gelman-Rubin diagnostic: 4 1.138070
#> Min Gelman-Rubin diagnostic: 2 1.000545
#> alpha.pi:
#> ----------
#> Max Gelman-Rubin diagnostic: 4 1.000380
#> Min Gelman-Rubin diagnostic: 2 1.000069
#> beta.pi:
#> --------
#> Max Gelman-Rubin diagnostic: 3 1.000475
#> Min Gelman-Rubin diagnostic: 2 1.000188
#>
#> Sampling Acceptance Rates:
#> ==========================
#> theta:
#> ------
#> [1] "Min: 0.207517, Max: 0.652350"mod.ptheta <- bhpm.ptheta(mod.fit)
print(mod.ptheta[mod.ptheta$ptheta.pos > .8 | mod.ptheta$ptheta.neg > 0.8, c("Trt.Grp","Cluster","Outcome","ptheta.pos","ptheta.neg")])
#> Trt.Grp Cluster Outcome ptheta.pos ptheta.neg
#> 5 2 M/0-64 G00-99_Outcome1 0.999515 1.0e-05
#> 6 2 M/0-64 G00-99_Outcome2 0.998815 1.5e-05
#> 14 2 M/65-84 G00-99_Outcome1 0.999845 0.0e+00
#> 23 2 M/85+ G00-99_Outcome1 1.000000 0.0e+00