glm.logml.npp: Log marginal likelihood of a GLM under normalized power prior...

In hdbayes: Bayesian Analysis of Generalized Linear Models with Historical Data

Log marginal likelihood of a GLM under normalized power prior (NPP)

Description

Uses bridge sampling to estimate the logarithm of the marginal likelihood of a GLM under the normalized power prior (NPP).

Usage

glm.logml.npp(post.samples, bridge.args = NULL)

Arguments

post.samples	output from glm.npp() giving posterior samples of a GLM under the normalized power prior (NPP), with an attribute called 'data' which includes the list of variables specified in the data block of the Stan program.
bridge.args	a list giving arguments (other than samples, log_posterior, data, lb, ub) to pass onto bridgesampling::bridge_sampler().

Value

The function returns a list with the following objects

model

"NPP"

logml

the estimated logarithm of the marginal likelihood

bs

an object of class bridge or bridge_list containing the output from using bridgesampling::bridge_sampler() to compute the logarithm of the marginal likelihood of the normalized power prior (NPP)

References

Duan, Y., Ye, K., and Smith, E. P. (2005). Evaluating water quality using power priors to incorporate historical information. Environmetrics, 17(1), 95–106.

Gronau, Q. F., Singmann, H., and Wagenmakers, E.-J. (2020). bridgesampling: An r package for estimating normalizing constants. Journal of Statistical Software, 92(10).

Examples

  if(requireNamespace("parallel")){
    data(actg019)
    data(actg036)
    ## take subset for speed purposes
    actg019 = actg019[1:100, ]
    actg036 = actg036[1:50, ]

    library(parallel)
    ncores    = 2
    data.list = list(data = actg019, histdata = actg036)
    formula   = cd4 ~ treatment + age + race
    family    = poisson()
    a0        = seq(0, 1, length.out = 11)
    if (instantiate::stan_cmdstan_exists()) {
      ## call created function
      ## wrapper to obtain log normalizing constant in parallel package
      logncfun = function(a0, ...){
        hdbayes::glm.npp.lognc(
          formula = formula, family = family, a0 = a0, histdata = data.list[[2]],
          ...
        )
      }

      cl = makeCluster(ncores)
      clusterSetRNGStream(cl, 123)
      clusterExport(cl, varlist = c('formula', 'family', 'data.list'))
      a0.lognc = parLapply(
        cl = cl, X = a0, fun = logncfun, iter_warmup = 500,
        iter_sampling = 1000, chains = 1, refresh = 0
      )
      stopCluster(cl)
      a0.lognc = data.frame( do.call(rbind, a0.lognc) )

      ## sample from normalized power prior
      d.npp = glm.npp(
        formula = formula,
        family = family,
        data.list = data.list,
        a0.lognc = a0.lognc$a0,
        lognc = matrix(a0.lognc$lognc, ncol = 1),
        chains = 1, iter_warmup = 500, iter_sampling = 1000,
        refresh = 0
      )
      glm.logml.npp(
        post.samples = d.npp,
        bridge.args = list(silent = TRUE)
      )
    }
  }

hdbayes documentation built on Sept. 11, 2024, 5:34 p.m.