pppvalue {bayesmeta} | R Documentation |

## 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 |

`parameter` |
the parameter to be tested. May be the effect
( |

`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 |

`rejection.region` |
the test statistic's rejection region. May be
one of |

`n` |
the number of Monte Carlo replications to be generated. The
default value is |

`prior` |
a logical flag to request |

`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 |

`...` |
further parameters passed to ‘ |

### 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 ( |

`p.value` |
the derived |

`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 |

`call` |
an object of class |

`rejection.region` |
the test statistic's rejection region. |

`replicates` |
a |

`computation.time` |
The computation time (in seconds) used. |

### Author(s)

Christian Roever christian.roever@med.uni-goettingen.de

### References

X.-L. Meng.
Posterior predictive p-values.
*The Annals of Statistics*, **22**(3):1142-1160, 1994.
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.
doi:10.18637/jss.v093.i06.

### See Also

### 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*version 3.4 Index]