Calculate plugin pD from a JAGS model with univariate likelihood for studies with repeated measurements

Uses results from MBNMA JAGS models to calculate pD via the plugin method (Spiegelhalter et al. 2002) . Can only be used for models with known standard errors or covariance matrices. Currently only functions with univariate likelihoods. Function is identical in MBNMAdose and MBNMAtime packages.

pDcalc(
  obs1,
  obs2,
  fups = NULL,
  narm,
  NS,
  theta.result,
  resdev.result,
  likelihood = "normal",
  type = "time"
)

Arguments

obs1

A matrix (study x arm) or array (study x arm x time point) containing observed data for y (normal likelihood) or r (binomial or poisson likelihood) in each arm of each study. This will be the same array used as data for the JAGS model.

obs2

A matrix (study x arm) or array (study x arm x time point) containing observed data for se (normal likelihood), n (binomial likelihood) or E (poisson likelihood) in each arm of each study. This will be the same array used as data for the JAGS model.

fups

A numeric vector of length equal to the number of studies, containing the number of follow-up mean responses reported in each study. Required for time-course MBNMA models (if type="time")

narm

A numeric vector of length equal to the number of studies, containing the number of arms in each study.

NS

A single number equal to the number of studies in the dataset.

theta.result

A matrix (study x arm) or array (study x arm x time point) containing the posterior mean predicted means/probabilities/rate in each arm of each study. This will be estimated by the JAGS model.

resdev.result

A matrix (study x arm) or array (study x arm x time point) containing the posterior mean residual deviance contributions in each arm of each study. This will be estimated by the JAGS model.

likelihood

A character object of any of the following likelihoods:

  • normal

  • binomial (does not work with time-course MBNMA models)

  • poisson (does not work with time-course MBNMA models)

type

The type of MBNMA model fitted. Can be either "time" or "dose"

Value

A single numeric value for pD calculated via the plugin method.

Details

Method for calculating pD via the plugin method proposed by Spiegelhalter (Spiegelhalter et al. 2002) . Standard errors / covariance matrices must be assumed to be known. To obtain values for theta.result and resdev.result these parameters must be monitored when running the MBNMA model (using parameters.to.save).

For non-linear time-course MBNMA models residual deviance contributions may be skewed, which can lead to non-sensical results when calculating pD via the plugin method. Alternative approaches are to use pV as an approximation or pD calculated by Kullback-Leibler divergence (Plummer 2008) .

References

Plummer M (2008). “Penalized loss functions for Bayesian model comparison.” Biostatistics, 9(3), 523-39. ISSN 1468-4357 (Electronic) 1465-4644 (Linking), https://pubmed.ncbi.nlm.nih.gov/18209015/.

Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A (2002). “Bayesian measures of model complexity and fit.” J R Statistic Soc B, 64(4), 583-639.

Examples

# \donttest{
# Using the triptans data
network <- mbnma.network(triptans)
#> Values for `agent` with dose = 0 have been recoded to `Placebo`
#> agent is being recoded to enforce sequential numbering

# Fit a dose-response MBNMA, monitoring "psi" and "resdev"
result <- mbnma.run(network, fun=dloglin(), method="random",
              parameters.to.save=c("psi", "resdev"))
#> `likelihood` not given by user - set to `binomial` based on data provided
#> `link` not given by user - set to `logit` based on assigned value for `likelihood`
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 182
#>    Unobserved stochastic nodes: 190
#>    Total graph size: 4074
#> 
#> Initializing model
#> 
#> Warning: error/missing in parameter psi in parameters.to.save, 
#>    Be aware of the output results.
#> Warning: error/missing in parameter resdev in parameters.to.save, 
#>    Be aware of the output results.


#### Calculate pD for binomial data ####

# Prepare data for pD calculation
r <- result$model$data()$r
n <- result$model$data()$n
narm <- result$model$data()$narm
NS <- result$model$data()$NS

psi <- result$BUGSoutput$median$psi
resdevs <- result$BUGSoutput$median$resdev

# Calculate pD via plugin method
pD <- pDcalc(obs1=r, obs2=n, narm=narm, NS=NS,
          theta.result=psi, resdev.result=resdevs,
          likelihood="binomial", type="dose")
# }