Hyperparameter Estimation with openEBGM

2026-09-03

Background

DuMouchel used an “Empirical Bayes” (EB) method. That is, the data are a driving force behind the choice of the prior distribution. One outcome of this is a known posterior distribution that relies on a set of hyperparameters derived from the prior distribution. These hyperparameters are estimated by their most likely value in the context of Empirical Bayesian methods. One process by which this can occur involves maximizing the likelihood function of the marginal distributions of the counts (or in our case, minimizing the negative log-likelihood function).

Optimization

Global optimization is a broad field. There are many existing R packages with minimization routines which may be used in the estimation of these hyperparameters. The hyperparameter estimation functions offered by openEBGM utilize the following:

openEBGM’s hyperparameter estimation functions use local optimization algorithms. The Newton-like approaches use algorithms implemented in R’s stats package and allow the user to choose multiple starting points to improve the chances of finding a global optimum. The user is encouraged to explore a variety of optimization approaches because the accuracy of a global optimization result is extremely difficult to verify and other approaches might work better in some cases.

References

  1. Meng X-L, Rubin D (1993). “Maximum likelihood estimation via the ECM algorithm: A general framework.” Biometrika, 80(2), 267-278.

  2. DuMouchel W, Pregibon D (2001). “Empirical Bayes Screening for Multi-item Associations.” In Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’01, pp. 67-76. ACM, New York, NY, USA. ISBN 1-58113-391-X.

  3. Millar, Russell B (2011). “Maximum Likelihood Estimation and Inference.” John Wiley & Sons, Ltd, 111-112.

Data Squashing

DuMouchel & Pregibon (2) discuss a methodology they call “data squashing”, which “compacts” the dataset to reduce the amount of computation needed to estimate the hyperparameters. openEBGM provides an implementation for data squashing.

The actual counts (\(N\)) and expected counts (\(E\)) are used to estimate the hyperparameters of the prior distribution. A large contingency table will have many cells, resulting in computational difficulties for the optimization routines needed for estimation. Data squashing (DuMouchel et al., 2001) transforms a set of 2-dimensional points to a smaller number of 3-dimensional points (triples). The idea is to reduce a large number of points \((N, E)\) to a smaller number of points \((N_k, E_k, W_k)\), where \(k = 1,...,M\) and \(W_k\) is the weight of the \(k^{th}\) “superpoint”. To minimize information loss, only points close to each other should be squashed.

For a given \(N\), squashData() combines points with similar \(E\)s into bins using a specified bin size and uses the average \(E\) within each bin as the \(E\) for that bin’s “superpoint”. The new superpoints are weighted by bin size. For example, the points (1, 1.1) and (1, 1.3) could be squashed to the superpoint (1, 1.2, 2).

An example is given below using unstratified expected counts:

library(openEBGM)
data(caers)

processed <- processRaw(caers)
processed[1:4, 1:4]
#>                                             var1                    var2 N
#> 1  7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS           ACUTE PSYCHOSIS 1
#> 2  7 ZEN 7-HYDROXYMITRAGYNINE CHEWABLE TABLETS           HOSPITALISATION 1
#> 3               7-HYDROXYMITRAGYNINE 7-OH KRATOM                   DEATH 1
#> 4                                     ARTRI KING BLOOD PRESSURE ABNORMAL 1
#>              E
#> 1 0.0003602305
#> 2 0.1253602305
#> 3 0.0407060519
#> 4 0.0021613833

squashed <- squashData(processed) #Using defaults
head(squashed)
#>   N            E weight
#> 1 1 0.0003602305     50
#> 2 1 0.0003602305     50
#> 3 1 0.0003602305     50
#> 4 1 0.0003602305     50
#> 5 1 0.0003602305     50
#> 6 1 0.0003602305     50

nrow(processed)
#> [1] 13441
nrow(squashed)
#> [1] 1202

As shown above, the squashed data set has 8.94% of the observations as the full dataset. Using this squashed data, we can then estimate the hyperparameters in a far more efficient manner.

squashData() can be used iteratively for each count (\(N\)):

squash1 <- squashData(processed)
squash2 <- squashData(squash1, count = 2, bin_size = 8)

Or autoSquash() can be used to squash all counts at once:

squash3 <- autoSquash(processed)
ftable(squash3[, c("N", "weight")])
#>    weight   1   2   8   9  24
#> N                            
#> 1         100   0   0   1 520
#> 2          76   0  55   0   0
#> 3          51  58   0   0   0
#> 4          70   0   0   0   0
#> 5          31   0   0   0   0
#> 6          22   0   0   0   0
#> 7          14   0   0   0   0
#> 8           6   0   0   0   0
#> 9           5   0   0   0   0
#> 10          6   0   0   0   0
#> 11          5   0   0   0   0
#> 12          2   0   0   0   0
#> 13          2   0   0   0   0
#> 16          1   0   0   0   0
#> 17          1   0   0   0   0
#> 19          1   0   0   0   0
#> 20          1   0   0   0   0
#> 23          1   0   0   0   0
#> 24          1   0   0   0   0

Likelihood Functions

As previously mentioned, the hyperparameters are estimated by minimizing the negative log-likelihood function. There are actually 4 different functions, depending on the use of data squashing and zero counts. All 4 functions, however, are based on the marginal distribution of the counts, which are mixtures of two negative binomial distributions. The hyperparameters are denoted by the vector \(\theta=(\alpha_1,\beta_1,\alpha_2,\beta_2,P)\), where \(P\) is the mixture fraction.

The most commonly used likelihood function is negLLsquash(), which is used when squashing data and not using zero counts. negLLsquash() is not called directly, but rather by some optimization function:

theta_init <- c(alpha1 = 1, beta1 = 1, alpha2 = 2, beta2 = 2, p = .2)
theta_hats <- stats::nlm(negLLsquash, p = theta_init,
  ni = squashed$N, ei = squashed$E, wi = squashed$weight, N_star = 1
)$estimate
#> Warning in log(alpha_beta[2]): NaNs produced
#> Warning in log(alpha_beta[2] + ei): NaNs produced
#> Warning in log(alpha_beta[2]): NaNs produced
#> Warning in log(alpha_beta[2] + ei): NaNs produced
#> Warning in log(1 - theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
#> Warning in log(alpha_beta[2]): NaNs produced
#> Warning in log(alpha_beta[2] + ei): NaNs produced
#> Warning in log(alpha_beta[2]): NaNs produced
#> Warning in log(alpha_beta[2] + ei): NaNs produced
#> Warning in log(1 - theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
#> Warning in log(alpha_beta[2]): NaNs produced
#> Warning in log(alpha_beta[2] + ei): NaNs produced
#> Warning in log(theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
#> Warning in log(theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
#> Warning in log(theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
#> Warning in log(theta[5]): NaNs produced
#> Warning in stats::nlm(negLLsquash, p = theta_init, ni = squashed$N, ei =
#> squashed$E, : NA/NaN replaced by maximum positive value
theta_hats
#> [1] 3.77387733 0.51353686 3.73731349 3.64827500 0.04818549

For unconstrained optimization, it is often better to work on the log and logit scales and then backtransform:

theta_init[1:4] <- log(theta_init[1:4])
theta_init[5] <- log(theta_init[5] / (1 - theta_init[5]))
theta_hats_transformed <- stats::nlm(negLLsquash, p = theta_init,
  ni = squashed$N, ei = squashed$E, wi = squashed$weight, N_star = 1,
  transformed = TRUE
)$estimate
theta_hats <- theta_hats_transformed
theta_hats[1:4] <- exp(theta_hats[1:4])
theta_hats[5] <- stats::plogis(theta_hats[5])
theta_hats
#> [1] 3.77418874 0.51358277 3.73737163 3.64834218 0.04818672

The N_star argument allows the user to choose the smallest value of \(N\) used for hyperparameter estimation. The user must be careful to match N_star with the actual minimum count in ni. If the user wishes to use a larger value for N_star, the vectors supplied to arguments ni, ei, and wi must be filtered. Here, we are using all counts except zeroes, so no filtering is needed. In general, N_star = 1 should be used whenever practical.

The other likelihood functions are negLL(), negLLzero(), and negLLzeroSquash(). Make sure to use the appropriate function for your choice of data squashing and use of zero counts.

Luckily, openEBGM offers some wrapper functions to help with this optimization task.

Special Wrapper Functions

Hyperparameters can be calculated by exploring the parameter space of the likelihood function using either the full data set of \(N\)s and \(E\)s or the squashed set. The methodology implemented by this package essentially maximizes the likelihood function (or more specifically and equivalently, minimizes the negative log-likelihood function). Starting points must be chosen to begin the exploration. DuMouchel (1999, 2001) provides a “lightly justified” set of initial hyperparameters. However, openEBGM’s functions support a large set of starting choices to help reach convergence and reduce the chance of false convergence. We begin by defining some starting points:

theta_init <- data.frame(
  alpha1 = c(1, 2, 3),
  beta1 = c(1, 2, 3),
  alpha2 = c(2, 4, 5),
  beta2 = c(2, 4, 5),
  p = c(.1, 0.2, 0.3)
)

The following functions (autoHyper(), exploreHypers(), and hyperEM()) automatically perform the transformations and backtransformations mentioned above.

autoHyper()

Now that the initial guesses for the hyperparameters have been defined, the function autoHyper() can be used to determine the actual hyperparameter estimates. autoHyper() performs a “verification check” by requiring the optimization routine to converge at least twice within the bounds of the parameter space. The estimate corresponding to the smallest negative log-likelihood is chosen as a tentative result. By default, this estimate must be similar to at least one other convergent solution. If the algorithm fails to consistently converge, an error message is returned.

system.time(
  hyper_estimates_full <- autoHyper(data = processed, theta_init = theta_init, 
    squashed = FALSE
  )
)
#> Trying method 'nlminb'...
#>   Working on initial guess number 1 of 3...
#>   Working on initial guess number 2 of 3...
#>   Working on initial guess number 3 of 3...
#>    user  system elapsed 
#>    8.02    0.02    8.70

squashed <- squashData(processed, count = 1, bin_size = 25, keep_pts = 10)
squashed <- squashData(squashed, count = 2, bin_size = 10, keep_pts = 10)
system.time(
  hyper_estimates_squashed <- autoHyper(data = squashed, theta_init = theta_init)
)
#> Trying method 'nlminb'...
#>   Working on initial guess number 1 of 3...
#>   Working on initial guess number 2 of 3...
#>   Working on initial guess number 3 of 3...
#>    user  system elapsed 
#>    0.40    0.01    0.52

hyper_estimates_full
#> $method
#> [1] "nlminb"
#> 
#> $estimates
#>     alpha1      beta1     alpha2      beta2          P 
#> 3.77518082 0.51366516 3.73386169 3.64522042 0.04817939 
#> 
#> $conf_int
#> NULL
#> 
#> $num_close
#> [1] 2
#> 
#> $theta_hats
#>   guess_num   a1_hat    b1_hat   a2_hat   b2_hat      p_hat code converge
#> 1         1 3.775181 0.5136652 3.733862 3.645220 0.04817939    0     TRUE
#> 2         2 3.775167 0.5136651 3.733887 3.645248 0.04818000    0     TRUE
#> 3         3 3.775175 0.5136644 3.733860 3.645219 0.04817937    0     TRUE
#>   in_bounds minimum
#> 1      TRUE 2995.19
#> 2      TRUE 2995.19
#> 3      TRUE 2995.19
hyper_estimates_squashed
#> $method
#> [1] "nlminb"
#> 
#> $estimates
#>     alpha1      beta1     alpha2      beta2          P 
#> 3.77687131 0.51432212 3.75499924 3.66462742 0.04828858 
#> 
#> $conf_int
#> NULL
#> 
#> $num_close
#> [1] 2
#> 
#> $theta_hats
#>   guess_num   a1_hat    b1_hat   a2_hat   b2_hat      p_hat code converge
#> 1         1 3.776871 0.5143221 3.754999 3.664627 0.04828858    0     TRUE
#> 2         2 3.776847 0.5143192 3.755000 3.664633 0.04828893    0     TRUE
#> 3         3 3.777163 0.5143511 3.755022 3.664659 0.04828866    0     TRUE
#>   in_bounds  minimum
#> 1      TRUE 2994.068
#> 2      TRUE 2994.068
#> 3      TRUE 2994.068

As seen above, the process is much faster when utilizing the squashed data, with estimates that are nearly identical. Of course, the amount of efficiency increase depends on the parameter values used in the call to squashData() and the size of the original data set. Another factor affecting run time is the number of starting points used. In general, you should use five or fewer starting points and limit the number of data points to a maximum of about 20,000 if possible. Notice that we squashed the same data again, which is allowed if we use a different value for count each time.

autoHyper() utilizes multiple minimization techniques to verify convergence. It first attempts the stats::nlminb() function, which implements a quasi-Newton unconstrained optimization technique using PORT routines. If nlminb() fails to consistently converge, autoHyper() next attempts stats::nlm(), which is a non-linear maximization algorithm based on the Newton method. Finally, if the first two approaches fail, a quasi-Newton method (also known as a variable metric algorithm) is used, which “…uses function values and gradients to build up a picture of the surface to be optimized.” (source: R documentation for stats::optim) This routine is implemented by the stats::optim() function with the method = BFGS argument.

autoHyper() can also return standard errors and confidence intervals for the hyperparameter estimates if needed:

autoHyper(squashed, theta_init = theta_init, conf_ints = TRUE)$conf_int
#> Trying method 'nlminb'...
#>   Working on initial guess number 1 of 3...
#>   Working on initial guess number 2 of 3...
#>   Working on initial guess number 3 of 3...
#>        pt_est     SE  LL_95  UL_95
#> a1_hat 3.7769 2.2950 0.0000 8.2751
#> b1_hat 0.5143 0.2453 0.0336 0.9950
#> a2_hat 3.7550 1.0022 1.7907 5.7193
#> b2_hat 3.6646 0.9139 1.8734 5.4559
#> p_hat  0.0483 0.0143 0.0202 0.0764

exploreHypers()

autoHyper() is a semi-automated approach to attempt estimation of the hyperparameters. The user is encouraged to explore various optimization approaches as no single approach will always work (also, the wrapper functions in openEBGM are not the only optimization functions available in R). One way to explore is with openEBGM’s exploreHypers() function, which is actually called by autoHyper() behind-the-scenes. exploreHypers() can also return standard errors for the hyperparameter estimates if needed:

exploreHypers(data = squashed, theta_init = theta_init, std_errors = TRUE)
#>   Working on initial guess number 1 of 3...
#>   Working on initial guess number 2 of 3...
#>   Working on initial guess number 3 of 3...
#> $estimates
#>   guess_num   a1_hat    b1_hat   a2_hat   b2_hat      p_hat code converge
#> 1         1 3.776871 0.5143221 3.754999 3.664627 0.04828858    0     TRUE
#> 2         2 3.776847 0.5143192 3.755000 3.664633 0.04828893    0     TRUE
#> 3         3 3.777163 0.5143511 3.755022 3.664659 0.04828866    0     TRUE
#>   in_bounds  minimum
#> 1      TRUE 2994.068
#> 2      TRUE 2994.068
#> 3      TRUE 2994.068
#> 
#> $std_errs
#>   guess_num    a1_se     b1_se    a2_se     b2_se       p_se
#> 1         1 2.295040 0.2452684 1.002227 0.9139328 0.01433702
#> 2         2 2.295058 0.2452673 1.002229 0.9139369 0.01433744
#> 3         3 2.295436 0.2453012 1.002246 0.9139586 0.01433784

exploreHypers() offers three gradient-based optimization methods and is basically just a wrapper around commonly used functions from the stats package mentioned earlier: nlminb() (the default), nlm(), & optim().

Expectation-Maximization Approach with hyperEM()

hyperEM() implements a version of the EM algorithm referred to by Meng & Rubin (1) as the Expectation/Conditional Maximization (ECM) algorithm. hyperEM() uses a single starting point to find a local maximum likelihood estimate for \(\theta\) either by finding roots of the score functions (partial derivatives of the log-likelihood function) or by using stats::nlminb() to directly minimize the negative log-likelihood function. The method described by Millar (3) is used to accelerate the estimate of \(\theta\) every 100 iterations of the algorithm.

data(caers)
processed <- processRaw(caers)
squashed <- squashData(processed, count = 1, bin_size = 25, keep_pts = 10)
squashed <- squashData(squashed, count = 2, bin_size = 10, keep_pts = 10)
hyperEM_ests <- hyperEM(squashed, theta_init_vec = c(1, 1, 2, 2, .1),
  conf_ints = TRUE, LL_tol = 1e-6, track = TRUE
)
#> change in LL: 1.47e+00 | theta: 3.377 0.597 2.074 2.347 0.065 
#> change in LL: 9.53e-02 | theta: 4.346 0.557 2.602 2.657 0.042 
#> change in LL:  2.9e-02 | theta: 4.378 0.548 2.834 2.839 0.041 
#> change in LL: 1.71e-02 | theta: 4.287 0.543 2.994 2.979 0.042 
#> change in LL: 1.08e-02 | theta: 4.202 0.539 3.123 3.094 0.043 
#>     Current iterations: 50 
#> change in LL: 7.01e-03 | theta: 4.132 0.535 3.229 3.189 0.044 
#> change in LL: 4.62e-03 | theta: 4.075 0.532 3.316 3.267 0.045 
#> change in LL: 3.09e-03 | theta: 4.027 0.53  3.388 3.332 0.045 
#> change in LL: 2.08e-03 | theta: 3.988 0.528 3.449 3.387 0.046 
#> change in LL: 3.88e-02 | theta: 3.786 0.514 3.75  3.661 0.048 
#>     Current iterations: 100 
#> change in LL: 4.22e-07 | theta: 3.78  0.514 3.752 3.662 0.048 
#> change in LL: 1.74e-07 | theta: 3.779 0.514 3.752 3.662 0.048 
#> change in LL: 1.14e-07 | theta: 3.779 0.514 3.753 3.663 0.048 
#> change in LL: 7.03e-08 | theta: 3.778 0.514 3.754 3.663 0.048 
#> change in LL:  4.1e-08 | theta: 3.778 0.514 3.754 3.664 0.048 
#>     Current iterations: 150 
#> change in LL: 2.17e-08 | theta: 3.778 0.514 3.754 3.664 0.048 
#> change in LL: 9.31e-09 | theta: 3.778 0.514 3.755 3.664 0.048 
#> change in LL: 1.59e-09 | theta: 3.778 0.514 3.755 3.665 0.048 
#> change in LL: -3.09e-09 | theta: 3.778 0.514 3.755 3.665 0.048 
#> change in LL: -5.64e-09 | theta: 3.778 0.514 3.755 3.665 0.048 
#>     Current iterations: 200 
#> 
#>    Iterations used: 208 
#> 
#>    Timing:  
#>    user  system elapsed 
#>    3.39    0.09    4.11
str(hyperEM_ests)
#> List of 9
#>  $ estimates : num [1:5] 3.7775 0.5144 3.7555 3.6651 0.0483
#>  $ conf_int  :'data.frame':  5 obs. of  4 variables:
#>   ..$ pt_est: num [1:5] 3.7775 0.5144 3.7555 3.6651 0.0483
#>   ..$ SE    : num [1:5] 2.2957 0.2453 1.0025 0.9142 0.0143
#>   ..$ LL_95 : num [1:5] 0 0.0336 1.7907 1.8733 0.0202
#>   ..$ UL_95 : num [1:5] 8.2769 0.9953 5.7203 5.4569 0.0764
#>  $ maximum   : num -2994
#>  $ method    : chr "score"
#>  $ elapsed   : num 4.11
#>  $ iters     : num 208
#>  $ score     : num [1:5] -4.91e-05 -9.99e-04 9.45e-04 -1.62e-03 1.08e-03
#>  $ score_norm: num 0.00238
#>  $ tracking  :'data.frame':  209 obs. of  7 variables:
#>   ..$ iter  : int [1:209] 0 1 2 3 4 5 6 7 8 9 ...
#>   ..$ logL  : num [1:209] -3067 -3042 -3025 -3017 -3013 ...
#>   ..$ alpha1: num [1:209] 1 1.69 1.86 2.08 2.29 ...
#>   ..$ beta1 : num [1:209] 1 0.864 0.787 0.737 0.702 ...
#>   ..$ alpha2: num [1:209] 2 2.22 1.87 1.7 1.65 ...
#>   ..$ beta2 : num [1:209] 2 1.95 1.95 1.98 2.03 ...
#>   ..$ P     : num [1:209] 0.1 0.156 0.155 0.142 0.127 ...

Setting track = TRUE tracks the log-likelihood and hyperparameter estimates at each iteration, which we can plot to study the behavior of the algorithm:

library(ggplot2)
library(tidyr)
pdat <- pivot_longer(hyperEM_ests$tracking,
  cols = -iter, names_to = "metric", values_to = "value"
)
pdat$metric <- factor(pdat$metric, levels = unique(pdat$metric), ordered = TRUE)
ggplot(pdat, aes(x = iter, y = value)) +
  geom_line(linewidth = 1.1, col = "blue") +
  facet_grid(metric ~ ., scales = "free") +
  ggtitle("Convergence Assessment",
    subtitle = "Dashed red line indicates accelerated estimate"
  ) +
  labs(x = "Iteration Count", y = "Estimate") +
  geom_vline(
    xintercept = c(100, 200), linewidth = 1, linetype = 2, col = "red"
  )

Specialized Optimization Packages

Other packages that specialize in optimization, such as DEoptim, can alternatively be used along with the likelihood functions to estimate the hyperparameters:

library(DEoptim)
set.seed(123456)
theta_hat <- DEoptim(negLLsquash,
  lower = rep(1e-03, 5),
  upper = c(rep(5, 4), .999),
  control = DEoptim.control(
    itermax = 2000,
    reltol = 1e-04,
    steptol = 200,
    NP = 100,
    CR = 0.85,
    F = 0.75,
    trace = 25
  ),
  ni = squashed$N, ei = squashed$E, wi = squashed$weight
)
#> Iteration: 25 bestvalit: 2994.539209 bestmemit:    3.957522    3.908413    3.355530    0.512015    0.948343
#> Iteration: 50 bestvalit: 2994.074653 bestmemit:    3.676489    3.603541    3.842106    0.522877    0.951507
#> Iteration: 75 bestvalit: 2994.068495 bestmemit:    3.764289    3.672363    3.757569    0.511906    0.951742
#> Iteration: 100 bestvalit: 2994.068360 bestmemit:    3.752504    3.662191    3.774799    0.513865    0.951743
#> Iteration: 125 bestvalit: 2994.068352 bestmemit:    3.754901    3.664609    3.777355    0.514387    0.951710
#> Iteration: 150 bestvalit: 2994.068352 bestmemit:    3.754906    3.664537    3.776929    0.514325    0.951712
#> Iteration: 175 bestvalit: 2994.068352 bestmemit:    3.754976    3.664606    3.776856    0.514319    0.951712
#> Iteration: 200 bestvalit: 2994.068352 bestmemit:    3.754972    3.664603    3.776872    0.514320    0.951712
#> Iteration: 225 bestvalit: 2994.068352 bestmemit:    3.754973    3.664604    3.776869    0.514320    0.951712
(theta_hat <- as.numeric(theta_hat$optim$bestmem))
#> [1] 3.7549733 3.6646042 3.7768695 0.5143198 0.9517117

Note that the two components in the mixture distribution switched here, which is not a problem.

Although the user is encouraged to explore various approaches for maximum likelihood hyperparameter estimation, autoHyper() will often give reasonable results with minimal effort. Once the hyperparameters have been estimated, they can be used in the calculation of the \(EBGM\) and quantile scores by applying them to the posterior distribution. This process can be found in the Empirical Bayes Metrics with openEBGM vignette.