normalmixture: Compute normal mixtures In bayesmeta: Bayesian Random-Effects Meta-Analysis

Description

This function allows to derive density, distribution function and quantile function of a normal mixture with fixed mean (μ) and random standard deviation (σ).

Usage

 ```1 2 3 4 5 6``` ``` normalmixture(density, cdf = Vectorize(function(x){integrate(density,0,x)\$value}), mu = 0, delta = 0.01, epsilon = 0.0001, rel.tol.integrate = 2^16*.Machine\$double.eps, abs.tol.integrate = rel.tol.integrate, tol.uniroot = rel.tol.integrate) ```

Arguments

 `density` the σ mixing distribution's probability density `function`. `cdf` the σ mixing distribution's cumulative distribution `function`. `mu` the normal mean (μ). `delta, epsilon` the parameters specifying the desired accuracy for approximation of the mixing distribution, and with that determining the number of σ support points being used internally. Smaller values imply greater accuracy and greater computational burden (Roever and Friede, 2017). `rel.tol.integrate, abs.tol.integrate, tol.uniroot` the `rel.tol`, `abs.tol` and `tol` ‘accuracy’ arguments that are passed to the `integrate()` or `uniroot()` functions for internal numerical integration or root finding (see also the help there).

Details

When a normal random variable

X|mu,sigma ~ Normal(mu, sigma^2)

has a fixed mean μ, but a random standard deviation

sigma|phi ~ G(phi)

following some probability distribution G(phi), then the marginal distribution of X,

X|mu,phi

is a mixture distribution (Lindsay, 1995; Seidel, 2010).

The mixture distribution's probability density function etc. result from integration and often are not available in analytical form. The `normalmixture()` function approximates density, distribution function and quantile function to some pre-set accuracy by a discrete mixture of normal distributions based on a finite number of σ values using the ‘DIRECT’ algorithm (Roever and Friede, 2017).

Either the “`density`” or “`cdf`” argument needs to be supplied. If only “`density`” is given, then the CDF is derived via integration, if only “`cdf`” is supplied, then the density function is not necessary.

In meta-analysis applications, mixture distributions arise e.g. in the context of prior predictive distributions. Assuming that a study-specific effect theta[i] a priori is distributed as

theta[i]|mu,tau ~ Normal(mu, tau^2)

with a prior distribution for the heterogeneity τ,

tau|phi ~ G(phi)

yields a setup completely analogous to the above one.

Since it is sometimes hard to judge what constitutes a sensible heterogeneity prior, it is often useful to inspect the implications of certain settings in terms of the corresponding prior predictive distribution of

theta[i]|mu,phi

indicating the a priori implied variation between studies due to heterogeneity alone based on a certain prior setup (Spiegelhalter et al., 2004, Sec. 5.7.3). Some examples using different heterogeneity priors are illustrated below.

Value

A `list` containing the following elements:

 `delta, epsilon` The supplied design parameters. `mu` the normal mean. `bins` the number of bins used. `support` the matrix containing lower, upper and reference points for each bin and its associated probability. `density` the mixture's density `function(x)`. `cdf` the mixture's cumulative distribution `function(x)` (CDF). `quantile` the mixture's quantile `function(p)` (inverse CDF). `mixing.density` the mixing distribution's density `function()` (if supplied). `mixing.cdf` the mixing distribution's cumulative distribution `function()`.

Author(s)

Christian Roever [email protected]

References

B.G. Lindsay. Mixture models: theory, geometry and applications. Vol. 5 of CBMS Regional Conference Series in Probability and Statistics, Institute of Mathematical Statistics, Hayward, CA, USA, 1995.

C. Roever, T. Friede. Discrete approximation of a mixture distribution via restricted divergence. Journal of Computational and Graphical Statistics, 26(1):217-222, 2017.

C. Roever. Bayesian random-effects meta-analysis using the bayesmeta R package. arXiv preprint 1711.08683 (submitted for publication), 2017.

W.E. Seidel. Mixture models. In M. Lovric (ed.) International Encyclopedia of Statistical Science, Springer, Heidelberg, pp. 827-829, 2010.

D.J. Spiegelhalter, K.R. Abrams, J.P.Myles. Bayesian approaches to clinical trials and health-care evaluation. Wiley & Sons, 2004.

`bayesmeta`.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95``` ```################################################################## # compare half-normal mixing distributions with different scales: nm05 <- normalmixture(cdf=function(x){phalfnormal(x, scale=0.5)}) nm10 <- normalmixture(cdf=function(x){phalfnormal(x, scale=1.0)}) # (this corresponds to the case of assuming a half-normal prior # for the heterogeneity tau) # check the structure of the returned object: str(nm05) # show density functions: # (these would be the marginal (prior predictive) distributions # of study-specific effects theta[i]) x <- seq(-1, 3, by=0.01) plot(x, nm05\$density(x), type="l", col="blue", ylab="density") lines(x, nm10\$density(x), col="red") abline(h=0, v=0, col="grey") # show cumulative distributions: plot(x, nm05\$cdf(x), type="l", col="blue", ylab="CDF") lines(x, nm10\$cdf(x), col="red") abline(h=0:1, v=0, col="grey") # determine 5 percent and 95 percent quantiles: rbind("HN(0.5)"=nm05\$quantile(c(0.05,0.95)), "HN(1.0)"=nm10\$quantile(c(0.05,0.95))) ################################################################## # compare different mixing distributions # (half-normal, half-Cauchy, exponential and Lomax): nmHN <- normalmixture(cdf=function(x){phalfnormal(x, scale=0.5)}) nmHC <- normalmixture(cdf=function(x){phalfcauchy(x, scale=0.5)}) nmE <- normalmixture(cdf=function(x){pexp(x, rate=2)}) nmL <- normalmixture(cdf=function(x){plomax(x, shape=4, scale=2)}) # show densities (logarithmic y-axis): x <- seq(-1, 3, by=0.01) plot(x, nmHN\$density(x), col="green", type="l", ylab="density", ylim=c(0.005, 6.5), log="y") lines(x, nmHC\$density(x), col="red") lines(x, nmE\$density(x), col="blue") lines(x, nmL\$density(x), col="cyan") abline(v=0, col="grey") # show CDFs: plot(x, nmHN\$cdf(x), col="green", type="l", ylab="CDF", ylim=c(0,1)) lines(x, nmHC\$cdf(x), col="red") lines(x, nmE\$cdf(x), col="blue") lines(x, nmL\$cdf(x), col="cyan") abline(h=0:1, v=0, col="grey") # add "exponential" x-axis at top: axis(3, at=log(c(0.5,1,2,5,10,20)), lab=c(0.5,1,2,5,10,20)) # show 95 percent quantiles: abline(h=0.95, col="grey", lty="dashed") abline(v=nmHN\$quantile(0.95), col="green", lty="dashed") abline(v=nmHC\$quantile(0.95), col="red", lty="dashed") abline(v=nmE\$quantile(0.95), col="blue", lty="dashed") abline(v=nmL\$quantile(0.95), col="cyan", lty="dashed") rbind("half-normal(0.5)"=nmHN\$quantile(0.95), "half-Cauchy(0.5)"=nmHC\$quantile(0.95), "exponential(2.0)"=nmE\$quantile(0.95), "Lomax(4,2)" =nmL\$quantile(0.95)) ##################################################################### # a normal mixture distribution example where the solution # is actually known analytically: the Student-t distribution. # If Y|sigma ~ N(0,sigma^2), where sigma = sqrt(k/X) # and X|k ~ Chi^2(df=k), # then the marginal Y|k is Student-t with k degrees of freedom. # define CDF of sigma: CDF <- function(sigma, df){pchisq(df/sigma^2, df=df, lower.tail=FALSE)} # numerically approximate normal mixture (with k=5 d.f.): k <- 5 nmT1 <- normalmixture(cdf=function(x){CDF(x, df=k)}) # in addition also try a more accurate approximation: nmT2 <- normalmixture(cdf=function(x){CDF(x, df=k)}, delta=0.001, epsilon=0.00001) # check: how many grid points were required? nmT1\$bins nmT2\$bins # show true and approximate densities: x <- seq(-2,10,le=400) plot(x, dt(x, df=k), type="l") abline(h=0, v=0, col="grey") lines(x, nmT1\$density(x), col="red", lty="dashed") lines(x, nmT2\$density(x), col="blue", lty="dotted") # show ratios of true and approximate densities: plot(x, nmT1\$density(x)/dt(x, df=k), col="red", type="l", log="y", ylab="density ratio") abline(h=1, v=0, col="grey") lines(x, nmT2\$density(x)/dt(x, df=k), col="blue") ```