Open schuemie opened 3 years ago
Convert the posterior samples to 1D density estimates right before here: https://github.com/OHDSI/EvidenceSynthesis/blob/aee99c04c8d6869fb0bf4d46a3662ea94a768b8a/R/BayesianMetaAnalysis.R#L273
stats::density
does this nicely.
But, since the posterior is 2D (mu, tau), you might want to use the dependence between parameters. In which case, just return the samples themselves (or consider a 2D density estimator).
For the particular use case I have in mind I only care about mu: In a network study, each site produces likelihood profiles for negative controls. For each negative control, we combine the per-site likelihood profiles using our Bayesian meta-analysis to produce a cross-network posterior distribution for the effect size of the negative control. I would then like to use these posterior distributions (without summarizing them as point estimates and standard errors) to fit an empirical null distribution. The fitNullNonNormalLl() function in the EmpiricalCalibration
package already takes likelihood profiles as input, so it would be convenient if we summarize the posterior distributions as such.
I could also modify the fitNullNonNormalLl()
to work with the samples, but at some point these samples will need to be converted into a density estimate, so why not in this package?
So perhaps a reasonable approach would be to use stats:density
with many points, and then use the adaptive profiling procedure in Cyclops
to reduce the number of points to the minimum required.
For your use-case, stats::density
will return a good approximation to the marginal distribution p(\mu | Y)
. Do not use Cyclops
to down-sample as that procedure assumes the function is convex (concave) which your posterior is certainly not (as a mixture of normals).
You can use these posteriors as input to EmpiricalCalibration
but keep in mind that your model does not fully follow from the rules of conditional probability; i.e. what you really want to is p(Y | \mu)
. An alternative approach that may not violate this small hitch is fitting a single, joint Bayesian model to all of the negative controls across all of the data sites at once. We could do this as well.
The
computeBayesianMetaAnalysis()
outputs a posterior distribution for the parameter of interest (i.e. effect size estimate). If we convert that to a likelihood profile we can use it for example to estimate a systematic error distribution.. For this, the output should be a data frame with two fields:point
andvalue
. We can use the adaptive profiling also implemented in Cyclops to do this.