pppvalue: Posterior predictive p-values

In bayesmeta: Bayesian Random-Effects Meta-Analysis and Meta-Regression

Posterior predictive p-values

Description

Compute posterior or prior predictive p-values from a bayesmeta object.

Usage

  pppvalue(x, parameter = "mu", value = 0.0,
           alternative = c("two.sided", "less", "greater"),
           statistic = "median",
           rejection.region,
           n = 10,
           prior = FALSE,
           quietly = FALSE,
           parallel, seed, ...)

Arguments

x	a bayesmeta object.
parameter	the parameter to be tested. May be the effect ("mu"), the heterogeneity ("tau") or one of the study-specific (\theta_i) parameters denoted by their label or their index.
value	the (null-) hypothesized value.
alternative	the type of alternative hypothesis.
statistic	the figure to be used a the ‘test statistic’, or ‘discrepancy variable’. May be chosen as "t", "Q" or "cdf", or among the row names of the bayesmeta object's ‘...$summary’ element. Or it may be specified as a function. For details, see below.
rejection.region	the test statistic's rejection region. May be one of "upper.tail", "lower.tail" or "two.tailed". If unspecified, it is set automatically based on the ‘alternative’ and ‘statistic’ parameters.
n	the number of Monte Carlo replications to be generated. The default value is n=10, but in practice a substantially larger value should be appropriate.
prior	a logical flag to request prior predictive (instead of posterior predictive) p-values. Prior predictive values are only available for hypotheses concerning the effect (\mu) and heterogeneity (\tau) parameters.
quietly	a logical flag to show (or suppress) output during computation; this may also speed up computations slightly.
parallel	the number of parallel processes to utilize. By default, if multiple (k) cores are detected, then k-1 parallel processes are used.
seed	(optional) an integer random seed value to generate reproducible results.
...	further parameters passed to ‘statistic’, if the ‘statistic’ argument was specified as a function.

Details

Posterior predictive p-values are Bayesian analogues to ‘classical’ p-values (Meng, 1994; Gelman, Meng and Stern, 1996; Gelman et al., 2014). The pppvalue() function allows to compute these values for one- and two-sided hypotheses concerning the effect (\mu) or heterogeneity (\tau) parameter, or one of the study-specific effect parameters (\theta_i) in a random-effects meta-analysis.

Prior predictive p-values have a similar interpretation, but they have a stronger dependence on the prior specification and are only available when the prior is proper; for a more detailed discussion, see Gelman, Meng and Stern (1996, Sec. 4).

The function may also be used to only generate samples (\tau, \mu, \theta_i, y_i) without having to also derive a statistic or a p-value. In order to achieve that, the ‘statistic’ argument may be specified as ‘NA’, and generated samples may be recovered from the ‘...$replicates’ output element.

p-values from Monte Carlo sampling

The computation is based on Monte Carlo sampling and repeated analysis of re-generated data sets drawn from the parameters' (conditional) posterior predictive (or prior) distribution; so the p-value derivation is somewhat computationally expensive. The p-value eventually is computed based on how far in the tail area the actual data are (in terms of the realized ‘test statistic’ or ‘discrepancy’) relative to the Monte-Carlo-sampled distribution. Accuracy of the computed p-value hence strongly depends on the number of samples (as specified through the ‘n’ argument) that are generated. Also, (near-) zero p-values need to be interpreted with caution, and in relation to the used Monte Carlo sample size (n).

‘Test’-statistics or ‘discrepancy variables’

The ‘statistic’ argument determines the statistic to be computed from the data as a measure of deviation from the null hypothesis. If specified as "Q", then Cochran's Q statistic is computed; this is useful for testing for homogeneity (\tau=0). If specified as one of the row names of the ‘x$summary’ element, then, depending on the type of null hypothesis specified through the ‘parameter’ argument, the corresponding parameter's posterior quantity is used for the statistic. If specified as "t", then a t-type statistic is computed (the difference between the corresponding parameter's posterior mean and its hypothesized value, divided by the posterior standard deviation). If specified as "cdf", the parameter's marginal posterior cumulative distribution function evaluated a the hypothesized value (‘value’) is used.

The ‘statistic’ argument may also be specified as an arbitrary function of the data (y). The function's first argument then needs to be the data (y), additional arguments may be passed as arguments (‘...’) to the ‘pppvalue()’ function. See also the examples below.

One- and two-sided hypotheses

Specification of one- or two-sided hypotheses not only has implications for the determination of the p-value from the samples, but also for the sampling process itself. Parameter values are drawn from a subspace according to the null hypothesis, which for a two-sided test is a line, and for a one-sided test is a half-plane. This also implies that one- and two-sided p-values cannot simply be converted into one another.

For example, when specifying pppvalue(..., param="mu", val=0, alt="two.sided"), then first paramater values (\tau, \mu) are drawn from the conditional posterior distribution p(\tau, \mu | y, \sigma, \mu=0), and subsequently new data sets are generated based on the parameters. If a one-sided hypothesis is specified, e.g. via pppvalue(..., param="mu", val=0, alt="less"), then parameters are drawn from p(\tau, \mu | y, \sigma, \mu>0).

For a hypothesis concerning the individual effect parameters \theta_i, conditions are imposed on the corresponding \theta_i. For example, for a specification of pppvalue(..., param=2, val=0, alt="less"), the hypothesis concerns the i=2nd study's effect paramater \theta_2. First a sample is generated from p(\theta_2|y, \sigma, \theta_2 > 0). Then samples of \mu and \tau are generated by conditioning on the generated \theta_2 value, and data y are generated by conditioning on all three.

Unless explicitly specified through the ‘rejection.region’ argument, the test statistic's “rejection region” (the direction in which extreme statistic values indicate a departure from the null hypothesis) is set based on the ‘alternative’ and ‘statistic’ parameters. The eventually used setting can be checked in the output's ‘...$rejection.region’ component.

Computation

When aiming to compute a p-value, it is probably a good idea to first start with a smaller ‘n’ argument to get a rough idea of the p-value's order of magnitude as well as the computational speed, before going over to a larger, more realistic n value. The implementation is able to utilize multiple processors or cores via the parallel package; details may be specified via the ‘parallel’ argument.

Value

A list of class ‘htest’ containing the following components:

statistic	the ‘test statistic’ (or ‘discrepancy’) value based on the actual data.
parameter	the number (n) of Monte Carlo replications used.
p.value	the derived p-value.
null.value	the (null-) hypothesized parameter value.
alternative	the type of alternative hypothesis.
method	a character string indicating what type of test was performed.
data.name	the name of the underlying bayesmeta object.
call	an object of class call giving the function call that generated the htest object.
rejection.region	the test statistic's rejection region.
replicates	a list containing the replicated parameters (\tau, \mu, \theta_i), data (y_i) and statistic, along with an indicator for those samples constituting the distribution's ‘tail area’.
computation.time	The computation time (in seconds) used.

Author(s)

Christian Roever [e-mail]

References

X.-L. Meng. Posterior predictive p-values. The Annals of Statistics, 22(3):1142-1160, 1994. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1214/aos/1176325622")}.

A. Gelman, X.-L. Meng, H. Stern. Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4):733-760, 1996.

A. Gelman, J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, D.B. Rubin. Bayesian data analysis. Chapman & Hall / CRC, Boca Raton, 2014.

C. Roever. Bayesian random-effects meta-analysis using the bayesmeta R package. Journal of Statistical Software, 93(6):1-51, 2020. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v093.i06")}.

See Also

bayesmeta, prop.test.

Examples

## Not run: 
# perform a meta analysis;
# load data:
data("CrinsEtAl2014")
# compute effect sizes (log odds ratios) from count data
# (using "metafor" package's "escalc()" function):
require("metafor")
crins.srr <- escalc(measure="OR",
                    ai=exp.SRR.events,  n1i=exp.total,
                    ci=cont.SRR.events, n2i=cont.total,
                    slab=publication, data=CrinsEtAl2014, subset=c(1,4,6))
# analyze:
bma <- bayesmeta(crins.srr, mu.prior.mean=0, mu.prior.sd=4,
                 tau.prior=function(t){dhalfnormal(t, scale=0.5)})

# compute a 2-sided p-value for the effect (mu) parameter
# (note: this may take a while!):
p <- pppvalue(bma, parameter="mu", value=0, n=100)

# show result:
print(p)

# show the test statistic's distribution
# along with its actualized value:
plot(ecdf(p$replicates$statistic[,1]),
     xlim=range(c(p$statistic, p$replicates$statistic[,1])))
abline(v=p$statistic, col="red")

# show the parameter values
# drawn from the (conditional) posterior distribution:
plot(bma, which=2)
abline(h=p$null.value)                                # (the null-hypothesized mu value)
points(p$replicates$tau, p$replicates$mu, col="cyan") # (the samples)

######################################################################
#  Among the 3 studies, only the first (Heffron, 2003) was randomized.
#  One might wonder about this particular study's effect (theta[1])
#  in the light of the additional evidence and compute a one-sided
#  p-value:

p <- pppvalue(bma, parameter="Heffron", value=0, n=100, alternative="less")
print(p)

######################################################################
#  One may also define one's own 'test' statistic to be used.
#  For example, one could utilize the Bayes factor to generate
#  a p-value for the homogeneity (tau=0) hypothesis:

BF <- function(y, sigma)
{
  bm <- bayesmeta(y=y, sigma=sigma,
                  mu.prior.mean=0, mu.prior.sd=4,
                  tau.prior=function(t){dhalfnormal(t, scale=0.5)},
                  interval.type="central")
  # (central intervals are faster to compute;
  #  interval type otherwise is not relevant here)
  return(bm$bayesfactor[1,"tau=0"])
}
# NOTE: the 'bayesmeta()' arguments above should probably match
#       the specifications from the original analysis

p <- pppvalue(bma, parameter="tau", statistic=BF, value=0, n=100,
              alternative="greater", rejection.region="lower.tail",
              sigma=bma$sigma)
print(p)

######################################################################
#  If one is only interested in generating samples (and not in test
#  statistics or p-values), one may specify the 'statistic' argument
#  as 'NA'.
#  Note that different 'parameter', 'value' and 'alternative' settings
#  imply different sampling schemes.

p <- pppvalue(bma, parameter="mu", statistic=NA, value=0,
              alternative="less", n=100)

plot(bma, which=2)
abline(h=p$null.value)
points(p$replicates$tau, p$replicates$mu, col="cyan")

## End(Not run)

bayesmeta documentation built on July 9, 2023, 5:12 p.m.