Solve the EBNM problem

Solves the empirical Bayes normal means (EBNM) problem using a specified family of priors. For an article-length introduction to the package, see Willwerscheid and Stephens (2021), cited in References below.

ebnm(
  x,
  s = 1,
  prior_family = c("point_normal", "point_laplace", "point_exponential", "normal",
    "horseshoe", "normal_scale_mixture", "unimodal", "unimodal_symmetric",
    "unimodal_nonnegative", "unimodal_nonpositive", "generalized_binary", "npmle",
    "deconvolver", "flat", "point_mass", "ash"),
  mode = 0,
  scale = "estimate",
  g_init = NULL,
  fix_g = FALSE,
  output = ebnm_output_default(),
  optmethod = NULL,
  control = NULL,
  ...
)

ebnm_output_default()

ebnm_output_all()

Arguments

x

A vector of observations. Missing observations (NAs) are not allowed.

s

A vector of standard errors (or a scalar if all are equal). Standard errors may not be exactly zero, and missing standard errors are not allowed. Two prior families have additional restrictions: when horseshoe priors are used, errors must be homoskedastic; and since function deconv in package deconvolveR takes $z$-scores, the "deconvolver" family requires that all standard errors be equal to 1.

prior_family

A character string that specifies the prior family $G$. See Details below.

mode

A scalar specifying the mode of the prior $g$ or "estimate" if the mode is to be estimated from the data. This parameter is ignored by the NPMLE, the deconvolveR family, and the improper uniform (or "flat") prior. For generalized binary priors, which are bimodal, the mode parameter specifies the mode of the truncated normal component (the location of the point mass is fixed at zero).

scale

A scalar or vector specifying the scale parameter(s) of the prior or "estimate" if the scale parameters are to be estimated from the data. This parameter is ignored by the flat prior and the family of point mass priors.

The interpretation of scale depends on the prior family. For normal and point-normal families, it is a scalar specifying the standard deviation of the normal component. For point-Laplace and point-exponential families, it is a scalar specifying the scale parameter of the Laplace or exponential component. For the horseshoe family, it corresponds to $s\tau$ in the usual parametrization of the horseshoe distribution. For the family of generalized binary priors, it specifies the ratio of the (untruncated) standard deviation of the normal component to its mode. This ratio must be fixed in advance (i.e., argument "estimate" is unavailable for generalized binary priors). For the NPMLE and deconvolveR prior family, scale is a scalar specifying the distance between successive means in the grid of point masses or normal distributions used to estimate $g$. For all other prior families, which are implemented using the function ash in package ashr, it is a vector specifying the parameter mixsd to be passed to ash or "estimate" if mixsd is to be chosen by ebnm. (Note that ebnm chooses mixsd differently from ash: see functions ebnm_scale_normalmix, ebnm_scale_unimix, and ebnm_scale_npmle for details. To use the ash grid, set scale = "estimate" and pass in gridmult as an additional parameter. See ash for defaults and details.)

g_init

The prior distribution $g$. Usually this is left unspecified (NULL) and estimated from the data. However, it can be used in conjuction with fix_g = TRUE to fix the prior (useful, for example, to do computations with the "true" $g$ in simulations). If g_init is specified but fix_g = FALSE, g_init specifies the initial value of $g$ used during optimization. For non-parametric priors, this has the side effect of fixing the mode and scale parameters. If g_init is supplied, it should be an object of class normalmix for normal, point-normal, scale mixture of normals, and deconvolveR prior families, as well as for the NPMLE; class laplacemix for point-Laplace families; class gammamix for point-exponential families; class horseshoe for horseshoe families; class unimix for unimodal_ families; or class tnormalmix for generalized binary priors. An object of class ebnm can also be supplied as argument, provided that field fitted_g contains a prior of the correct class (see Examples below).

fix_g

If TRUE, fix the prior $g$ at g_init instead of estimating it.

output

A character vector indicating which values are to be returned. Function ebnm_output_default() provides the default return values, while ebnm_output_all() lists all possible return values. See Value below.

optmethod

A string specifying which optimization function is to be used.

For parametric families other than the horseshoe and generalized binary (normal, point-normal, point-Laplace, and point-exponential), options include "nlm" (which calls nlm), "lbfgsb" (which calls optim with method = "L-BFGS-B"), and "trust" (which calls into the trust package). Other options are "nohess_nlm", "nograd_nlm", and "nograd_lbfgsb", which use numerical approximations rather than exact expressions for the Hessian; both of the "nograd" functions use numerical approximations for the gradient as well. The default option is "nohess_nlm".

Since the horseshoe, generalized binary, and point mass families only require one parameter to be estimated (at most), the only available optimization method is optimize, and thus the optmethod parameter is ignored by ebnm_horseshoe, ebnm_generalized_binary, and ebnm_point_mass.

For most nonparametric families (scale mixtures of normals; unimodal, symmetric unimodal, nonnegative unimodal, and nonpositive unimodal families; and the NPMLE), optmethod options are provided by package ashr. The default method uses the mix-SQP algorithm implemented in the mixsqp package. See the ash function documentation for other options. For the NPMLE only, it is also possible to specify optmethod = "REBayes", which uses function GLmix in the REBayes package to estimate the NPMLE rather than using the ashr package. Note that REBayes requires installation of the commercial interior-point solver MOSEK; for details, see the documentation for REBayes function KWDual.

The nonparametric exception is the the "deconvolveR" family. Since the deconvolveR package only ever uses nlm, ebnm_deconvolver ignores the optmethod parameter.

control

A list of control parameters to be passed to the optimization function specified by parameter optmethod.

...

Additional parameters. When a unimodal_ prior family is used, these parameters are passed to function ash in package ashr. Although it does not call into ashr, the scale mixture of normals family accepts parameter gridmult for purposes of comparison. When gridmult is set, an ashr-style grid will be used instead of the default ebnm grid. When the "deconvolver" family is used, additional parameters are passed to function deconv in package deconvolveR. Families of generalized binary priors take several additional parameters; see ebnm_generalized_binary. In all other cases, additional parameters are ignored.

Value

An ebnm object. Depending on the argument to output, the object is a list containing elements:

data: A data frame containing the observations x and standard errors s.
posterior: A data frame of summary results (posterior means, standard deviations, second moments, and local false sign rates).
fitted_g: The fitted prior $\hat{g}$ (an object of class normalmix, laplacemix, gammamix, unimix, tnormalmix, or horseshoe).
log_likelihood: The optimal log likelihood attained, $L(\hat{g})$.
posterior_sampler: A function that can be used to produce samples from the posterior. For all prior families other than the horseshoe, the sampler takes a single parameter nsamp, the number of posterior samples to return per observation. Since ebnm_horseshoe returns an MCMC sampler, it additionally takes parameter burn, the number of burn-in samples to discard.

S3 methods coef, confint, fitted, logLik,

nobs, plot, predict, print, quantile,

residuals, simulate, summary, and vcov

have been implemented for ebnm objects. For details, see the respective help pages, linked below under See Also.

Details

Given vectors of data x and standard errors s, ebnm solves the "empirical Bayes normal means" (EBNM) problem for various choices of prior family. The model is $$x_j | \theta_j, s_j \sim N(\theta_j, s_j^2)$$ $$\theta_j \sim g \in G,$$ where $g$, which is referred to as the "prior distribution" for $\theta$, is to be estimated from among some specified family of prior distributions $G$. Several options for $G$ are implemented, some parametric and others non-parametric; see below for examples.

Solving the EBNM problem involves two steps. First, $g \in G$ is estimated via maximum marginal likelihood: $$\hat{g} := \arg\max_{g \in G} L(g),$$ where $$L(g) := \prod_j \int p(x_j | \theta_j, s_j) g(d\theta_j).$$ Second, posterior distributions $p(\theta_j | x_j, s_j, \hat{g})$ and/or summaries such as posterior means and posterior second moments are computed.

Implemented prior families include:

point_normal: The family of mixtures where one component is a point mass at $\mu$ and the other is a normal distribution centered at $\mu$.
point_laplace: The family of mixtures where one component is a point mass at $\mu$ and the other is a double-exponential distribution centered at $\mu$.
point_exponential: The family of mixtures where one component is a point mass at $\mu$ and the other is a (nonnegative) exponential distribution with mode $\mu$.
normal: The family of normal distributions.
horseshoe: The family of horseshoe distributions.
normal_scale_mixture: The family of scale mixtures of normals.
unimodal: The family of all unimodal distributions.
unimodal_symmetric: The family of symmetric unimodal distributions.
unimodal_nonnegative: The family of unimodal distributions with support constrained to be greater than the mode.
unimodal_nonpositive: The family of unimodal distributions with support constrained to be less than the mode.
generalized_binary: The family of mixtures where one component is a point mass at zero and the other is a truncated normal distribution with lower bound zero and nonzero mode. See Liu et al. (2023), cited in References below.
npmle: The family of all distributions.
deconvolver: A non-parametric exponential family with a natural spline basis. Like npmle, there is no unimodal assumption, but whereas npmle produces spiky estimates for $g$, deconvolver estimates are much more regular. See deconvolveR-package for details and references.
flat: The "non-informative" improper uniform prior, which yields posteriors $$\theta_j | x_j, s_j \sim N(x_j, s_j^2).$$
point_mass: The family of point masses $\delta_\mu$. Posteriors are point masses at $\mu$.
ash: Calls into function ash in package ashr. Can be used to make direct comparisons of ebnm and ashr implementations of prior families such as scale mixtures of normals and the NPMLE.

Functions

ebnm_output_default(): Lists the default return values.
ebnm_output_all(): Lists all valid return values.

References

Jason Willwerscheid, Peter Carbonetto and Matthew Stephens (2023). ebnm: an R Package for solving the empirical Bayes normal means problem using a variety of prior families. arXiv:2110.00152.

Yusha Liu, Peter Carbonetto, Jason Willwerscheid, Scott A Oakes, Kay F Macleod, and Matthew Stephens (2023). Dissecting tumor transcriptional heterogeneity from single-cell RNA-seq data by generalized binary covariance decomposition. bioRxiv 2023.08.15.553436.

Examples

theta <- c(rep(0, 100), rexp(100))
s <- 1
x <- theta + rnorm(200, 0, s)

# The following are equivalent:
pn.res <- ebnm(x, s, prior_family = "point_normal")
pn.res <- ebnm_point_normal(x, s)

# Inspect results:
logLik(pn.res)
#> 'log Lik.' -341.071 (df=2)
plot(pn.res)


# Fix the scale parameter:
pl.res <- ebnm_point_laplace(x, s, scale = 1)

# Estimate the mode:
normal.res <- ebnm_normal(x, s, mode = "estimate")
plot(normal.res) # posterior means shrink to a value different from zero


# Use an initial g (this fixes mode and scale for ash priors):
normalmix.res <- ebnm_normal_scale_mixture(x, s, g_init = pn.res)

# Fix g and get different output (including a posterior sampler):
pn.res <- ebnm_point_normal(x, s, g_init = pn.res, fix_g = TRUE,
                            output = ebnm_output_all())

# Sample from the posterior:
pn.samp <- simulate(pn.res, nsim = 100)

# Quantiles and HPD confidence intervals can be obtained via sampling:
set.seed(1)
pn.quantiles <- quantile(pn.res, probs = c(0.1, 0.9))
pn.quantiles[1:5, ]
#>             10%       90%
#> [1,]  0.0000000 1.3023520
#> [2,] -0.4284940 0.1340129
#> [3,] -0.5275731 0.0000000
#> [4,] -0.9936367 0.0000000
#> [5,] -1.1610714 0.0000000
confint(pn.res, level = 0.8, parm = 1:5)
#>          CI.lower   CI.upper
#> [1,] -0.004626851 1.18406890
#> [2,] -0.064080962 0.45932241
#> [3,] -0.313072644 0.06839581
#> [4,] -0.276248083 0.27054739
#> [5,] -0.896940425 0.06111930

# Examples of usage of control parameter:
#  point_normal uses nlm:
pn.res <- ebnm_point_normal(x, s, control = list(print.level = 1))
#> iteration = 0
#> Step:
#> [1] 0 0
#> Parameter:
#> [1] 0.0000000 0.5954234
#> Function Value
#> [1] 157.5583
#> Gradient:
#> [1] -1.451161  1.608956
#> 
#> iteration = 7
#> Parameter:
#> [1] 0.5030917 0.7422532
#> Function Value
#> [1] 157.2833
#> Gradient:
#> [1]  2.123385e-05 -3.331885e-05
#> 
#> Relative gradient close to zero.
#> Current iterate is probably solution.
#> 
#  unimodal uses mixsqp:
unimodal.res <- ebnm_unimodal(x, s, control = list(verbose = TRUE))
#> Running mix-SQP algorithm 0.3-54 on 201 x 7 matrix
#> convergence tol. (SQP):     1.0e-08
#> conv. tol. (active-set):    1.0e-10
#> zero threshold (solution):  1.0e-08
#> zero thresh. (search dir.): 1.0e-14
#> l.s. sufficient decrease:   1.0e-02
#> step size reduction factor: 7.5e-01
#> minimum step size:          1.0e-08
#> max. iter (SQP):            1000
#> max. iter (active-set):     8
#> number of EM iterations:    20
#> Computing SVD of 201 x 7 matrix.
#> Matrix is not low-rank; falling back to full matrix.
#> iter        objective max(rdual) nnz stepsize max.diff nqp nls
#>    1 +6.994875996e-01  -- EM --    7 1.00e+00 1.04e-01  --  --
#>    2 +6.141425072e-01  -- EM --    7 1.00e+00 1.03e-01  --  --
#>    3 +5.635847714e-01  -- EM --    7 1.00e+00 8.83e-02  --  --
#>    4 +5.341728768e-01  -- EM --    7 1.00e+00 6.96e-02  --  --
#>    5 +5.169976722e-01  -- EM --    7 1.00e+00 5.24e-02  --  --
#>    6 +5.067452345e-01  -- EM --    7 1.00e+00 3.88e-02  --  --
#>    7 +5.004267543e-01  -- EM --    7 1.00e+00 2.85e-02  --  --
#>    8 +4.963935909e-01  -- EM --    7 1.00e+00 2.10e-02  --  --
#>    9 +4.937303221e-01  -- EM --    7 1.00e+00 1.55e-02  --  --
#>   10 +4.919174991e-01  -- EM --    7 1.00e+00 1.15e-02  --  --
#>   11 +4.906514399e-01  -- EM --    7 1.00e+00 8.54e-03  --  --
#>   12 +4.897485316e-01  -- EM --    7 1.00e+00 6.37e-03  --  --
#>   13 +4.890938508e-01  -- EM --    7 1.00e+00 4.77e-03  --  --
#>   14 +4.886130180e-01  -- EM --    7 1.00e+00 3.58e-03  --  --
#>   15 +4.882563857e-01  -- EM --    7 1.00e+00 2.70e-03  --  --
#>   16 +4.879898997e-01  -- EM --    7 1.00e+00 2.04e-03  --  --
#>   17 +4.877896543e-01  -- EM --    7 1.00e+00 1.54e-03  --  --
#>   18 +4.876385451e-01  -- EM --    6 1.00e+00 1.17e-03  --  --
#>   19 +4.875241460e-01  -- EM --    6 1.00e+00 8.84e-04  --  --
#>   20 +4.874373211e-01  -- EM --    6 1.00e+00 6.70e-04  --  --
#>    1 +4.873712915e-01 +1.732e-03   6  ------   ------   --  --
#>    2 +4.871340589e-01 +5.513e-04   4 1.00e+00 1.56e-02   4   1
#>    3 +4.871340386e-01 -1.884e-06   4 1.00e+00 1.69e-04   2   1
#> Optimization took 0.00 seconds.
#> Convergence criteria met---optimal solution found.