A Constrained Conditional Likelihood Approach for Estimating the Means of Selected Populations
Claudio Fuentes Vik Gopal


Polarized Heavy Quarkonium Production in the Color Evaporation Model
    Ramona Vogt
    December 30, 2023
====================================================================










Given p independent normal populations, we consider the problem of estimating
the mean of those populations, that based on the observed data, give the
strongest signals. We explicitly condition on the ranking of the sample means,
and consider a constrained conditional maximum likelihood (CCMLE) approach,
avoiding the use of any priors and of any sparsity requirement between the
population means.
Our results show that if the observed means are too close together, we should
in fact use the grand mean to estimate the mean of the population with the
larger sample mean. If they are separated by more than a certain threshold, we
should shrink the observed means towards each other. As intuition suggests,
it is only if the observed means are far apart that we should conclude that the
magnitude of separation and consequent ranking are not due to chance. 




Unlike other methods, our approach does not need to pre-specify the number of selected populations and the proposed
CCMLE is able to perform simultaneous inference.








Our method, which
is conceptually straightforward, can be easily adapted to incorporate other
selection criteria.










§ INTRODUCTION

Consider a scenario where p ≥ 2 independent normal populations are
available, and from each one of them, we obtain a sample of size n. In this
context, practitioners are sometimes interested in estimating the true means of
those populations that yielded the k largest sample means in the experiment. 










A naive solution to the problem is to estimate the selected population means
with the corresponding sample means. Such an approach, however, is known to be
problematic.  <cit.> showed that the resulting estimator is
biased for the case k=1. This bias is particularly evident when all p
populations are identically distributed. 

In terms of optimality, <cit.> both
showed that the estimator is minimax only when p=2.





In order to improve on the naive estimator, several alternatives have been
proposed in the literature, including <cit.>,
<cit.>, and <cit.>. These papers
propose estimators that perform better in terms of the Mean Squared Error (MSE).
<cit.> considered a bias correction approach for the
problem, obtaining estimators that perform well in terms of frequentist risk.
Following up on his own idea, <cit.> introduced 
ω-estimators which are essentially a weighted average of the order
statistics. <cit.> considered a two-stage procedure, assuming we
can obtain a second sample from the selected population, to produce unbiased
estimators of the selected means. Despite these results, performance theorems
are scarce, with the exception of  <cit.> and 
<cit.>. The latter proposes an
empirical Bayes estimator and shows that it performs better in terms of the
Bayes risk with respect to any normal prior. The former paper 
focuses on admissibility and obtains a generalized Bayes estimator under a
harmonic prior.


More recently, <cit.> make an implicit assumption of sparsity,
that many effect sizes (absolute population means) θ_i=0. By adapting 
the theory developed in <cit.>, the authors in <cit.> 
perform post-selection inference with the Lasso. <cit.>
approached the problem by estimating the first and second order bias of the
naive estimator. Their results, which are similar in performance to the
empirical Bayes approach in <cit.>, extend to the non-Gaussian
setting.




In this paper we propose a new estimator, that is based entirely on the
likelihood function after incorporating the selection process.  We motivate and
define this new estimator in section <ref>, where we also provide
a neat result for p=2 that yields some insight into how the estimator works.
In section <ref>, we discuss the main computational hurdle in
computing this estimator, and provide some direction on overcoming it. Following
that section, we summarize the results of a simulation study that highlights the
benefits and flaws of this new estimator. We conclude the paper with a brief
discussion of our main results.



§ CONDITIONAL LIKELIHOOD ESTIMATION





 §.§ Defining the Conditional Likelihood


For simplicity, suppose that we obtain a single observation X_i ∼
N(μ_i, σ^2) from each population, and that the common variance
σ^2=1. Then, for μ = (μ_1, …, μ_p)^T, the unconditional
likelihood is given by L_0(μ) = ∏_i=i^p ϕ(x_i - μ_i),
where ϕ denotes the density function of the standard normal distribution.
However, once we observe X_1=x_1, X_2=x_2, …, X_p=x_p, we rank the
observations in order to identify the populations corresponding to the largest
x_i's and estimate the respective means μ_i.











Hence, the conditional likelihood of interest is 

    L(μ) = 
    ∏_i=i^p ϕ(x_i - μ_i)/P_μ(X_1 > X_2 … > X_p),

where the notation P_μ(·) explicitly states that the probability
under consideration depends on μ. Note that, for equation
(<ref>), we do not have to worry about the labels attached to the
groups. In other words, there is no loss of generality in assuming that the
ordering of the means is X_1 > X_2 > … > X_p. It follows from equation
(<ref>) that the log of the conditional likelihood is




















    l(μ) =  C - 1/2∑_i=1^p (x_i - μ_i)^2 - 
    logP_μ(X_1 > X_2 … > X_p)


 

 §.§ Constrained Conditional MLE



For the case p=2, note that μ = (μ_1, μ_2) and 









    P_μ(X_1 > X_2) = P_μ(X_1 - X_2 > 0) = 1 - Φ( μ_2 - μ_1/√(2))

For a fixed μ_2, P_μ(X_1 > X_2) approaches 0 as μ_1↓-∞.
This means that the -log P_μ(X_1 > X_2) term in equation (<ref>)
increases to ∞, causing the conditional likelihood to be unbounded.  In
other words, there is no global maximum for the expression in 
(<ref>) for -∞ < μ_1, μ_2 < ∞. This phenomenon is
of course, not specific to p=2.

The occurrence of an unbounded likelihood is not without precedent in the
statistics literature. Possibly the most widely studied models in which this
occurs is a normal mixture model with unequal variances
<cit.>. One of the solutions in that model was to constrain
the parameter space to where there are local modes, and that is what we shall
attempt to do here as well.

Since we expect that the populations with the largest sample means would be
the ones with the largest population means, we now aim to maximise the
conditional likelihood in equation (<ref>) subject to 
μ∈Θ, where Θ = {μ : μ_1 ≥μ_2 ≥μ_3 …≥μ_p }.
We thus define the Constrained Conditional Maximum Likelihood Estimator (CCMLE) to be 

    μ̂ = μ∈Θ  
    l(μ)


This is clearly a case of constrained statistical inference. Much of the theory
regarding this approach can be found in  <cit.>.














 §.§ A Closed Form Result













































Consider the CCMLE when p=2 and the variance is known. Thus, we have
realisations x_1 and x_2 from N(μ_i, σ^2) for i=1,2, and have
observed that x_1 > x_2. If x_1 - x_2 > 2 σ/√(π), then
the unique CCMLE is an interior point of Θ = {(μ_1, μ_2): μ_1 ≥μ_2 }. If x_1 - x_2 ≤2 σ/√(π), the CCMLE is
μ̂ = (x̅, x̅), where x̅ = (x_1 + x_2)/2.








Let us first define g(·) to denote the inverse Mills ratio:

    g(x) = ϕ(x)/1 - Φ(x)

We shall use certain properties of g in this proof – in particular, note that 
g is convex and monotone increasing <cit.>.

Recall that the log-likelihood function is given by

    l(μ) = -log 2π - logσ^2 - 1/2σ^2∑_i=1^2 (x_i - μ_i)^2 - 
    log( 1 - Φ(μ_2 - μ_1/√(2 σ^2)) )





















Setting the partial derivatives to zero, we find that, if a solution exists, it
must satisfy







    g ( √(2/σ^2) (x̅ - μ_1) ) = 
    √(2/σ^2) ( x_1 - μ_1 )

Observe that equation (<ref>) depends only on μ_1. Let us define the
function on the left to be h_1(·) and the function on the right to be
h_2(·). Thus, a stationary point exists if h_1 and h_2 intersect. Note
that we are finding the root of a transcendental equation; no closed form
solution exists. A plot of these two functions can be seen in Figure
<ref>. Several quantities on the plot can be derived. For instance,
from equation (<ref>), we can compute that h_1(x̅) = g(0) =
√(2/π) and h'_1(x̅) = -g'(0) √(2 /σ^2)≈
-0.900/ √(σ^2).  Moreover, the height H is (x_1 - x_2)/√(2
σ^2).
 


The figure shows that an intersection exists in {μ_1: μ_1
≥x̅} if and only if H > √(2/π). In other words, 
if and only if x_1 - x_2 > 2 σ/√(π).

Suppose that there is more than 1 root in {μ_1: μ_1 ≥x̅}.
Then, at some point, it must be that h'_1(μ_1)< - √(2 / σ^2).
However h'_1(0) ≈ -0.900 / √(σ^2), and h_1 is convex.
Hence h'_1(μ_1) > -0.900 / √(σ^2) for all μ_1 ≥x̅. It follows from this contradiction that the stationary point (if it
exists) is unique.

Our final task is to show that, when x_1 - x_2 ≤ 2 σ/√(π),  the
CCMLE is on the boundary of Θ at (x̅, x̅). Within the
constrained parameter space Θ, we have μ_2 ≤μ_1 and therefore 







    1 - Φ( μ_2 - μ_1/√(2 σ^2)) ≥ 0.5


Hence, for all μ∈Θ, 






    l(μ) ≤ C - 1/2 σ^2∑_i=1^2 (x_i - μ_i)^2 + log 2

It is then easy to see that as μ_1 and/or μ_2 go to +∞ and/or
-∞, it holds that l(μ) → -∞.

When x_1 - x_2 < 2 σ/√(π), we already know that there are no stationary
points in Θ. This information, coupled with the understanding that
l(μ) decreases without bound as μ moves away from the
boundary, leads to the conclusion that the CCMLE must be on the boundary {μ: μ_1 = μ_2 }.

Now consider any point μ such that μ_1 = μ_2. If we can show that
the directional derivative in the direction of (x̅, x̅) is always
positive, we are done, because it implies that we can always increase the
log-likelihood by taking a suitable step in the direction of (x̅,
x̅). To simplify notation, let  

    γ = 1/√(2 σ^2) g ( μ_2 - μ_1/√(2
    σ^2))

Thus we can write the gradient of l as ∇ l = 
( 1/σ^2(x_1 - μ_1) - γ, 1/σ^2
(x_2 - μ_2) + γ)^T. The direction we need to consider is 
u = K(x̅ - μ_1, x̅ - μ_2)^T,
where K is a positive normalising constant, that makes |u| = 1. 

For any point in Θ, the directional derivative is

    ∇ l^T u∝ ((1/σ^2)(x_1 - μ_1) - γ)(x̅ - μ_1) + 
    ((1/σ^2)(x_2 - μ_2) + γ)(x̅ - μ_2)

Now for any point on the edge of Θ, we let μ = μ_1 = μ_2 and
simplify expression (<ref>) to show that 

    ∇ l^T u∝2/σ^2(x̅ - μ)^2 ≥ 0





 §.§ Sample Computations



It is not straightforward to generalize the proof in Proposition <ref> to
higher dimensions. However, the same intuitive results and numerical approaches apply. 

In Table <ref> below, we present some sample CCMLE
computations for various combinations of observed values when p=4. In
performing the computations, we assume that σ^2 is known, and is equal to 1.
The purpose is to underline how, whenever a subset of observed sample means are
close together, the CCMLE procedure will shrink them towards each other. This is
apparent in row 1 of the table. Unlike the p=2 case, however, the critical
distance at which they are collapsed onto one another is no longer
2/√(π). This latter phenomenon can be observed in row 2, where the
separation is less than 2/√(π), but the estimated values are not exactly
equal.






















§ COMPUTATION OF THE CCMLE





 §.§ Calculation of Probabilities

As part of numerically maximising the conditional likelihood in equation
(<ref>), we need to repeatedly compute P_μ(X_1 > X_2 …
> X_p) for different μ vectors.  In this section, we outline how it is
possible to compute this integral for small values of p using a trick of
conditioning, followed by a nested application of the 
function in R.



For the case p=2, it is clear from equation (<ref>) that the
probability of interest can be computed using the usual approximations to the
standard normal distribution function.

For the case p=3, we can condition on the random variable in the middle to
yield a univariate integral.

    P_μ(X_1 > X_2 > X_3)    =   ∫ P_μ(X_1 > X_2 > X_3 | X_2 = x_2) f(x_2)  x_2 
       =   ∫ P_μ(X_1 > x_2) P(x_2 > X_3) f(x_2)  x_2


























































The same approach enabled us to compute the probabilities up to p=7 without
any further optimisation in R. For higher dimensions, we advocate computing the
integral using a lower level language such as C, and switching to a sparse grid
method <cit.> instead of persisting with cubature
techniques.









 §.§ Obtaining A Good Starting Point




























We now focus on obtaining a good starting point for the numerical optimisation
in the case that p > 2. A good starting point ensures that we can reduce the
number of times that we evaluate the high-dimensional probability and its
derivatives. Let us first denote f(μ) = 
logP_μ(X_1 > X_2 … > X_p). Taking a first order Taylor
approximation of f about the observed x, we can approximate the
conditional likelihood in (<ref>) with 

    l(μ) ≈  C - 1/2  (x - μ)^T (x - μ) - 
    f(x) - (μ - x)^T ∇ f(x)


Taking derivative with respect to μ and setting the above equation to
0, we can get an approximation to the CCMLE. It works out to be 

    μ̂_0 ≈x - ∇ f(x)


The solution does not mean that a stationary point always exists - remember that
we are solving an approximation to the conditional log-likelihood. For the same
reason, it is also possible that the solution in equation (<ref>)
does not fall within the constrained space Θ. In such cases, we shall use
the orthogonal projection onto Θ as the starting point. Note that the
projection has to be performed numerically - fortunately however, Θ is
convex, and hence we can rely on prior methods from convex optimization theory
in order to perform this step. The R package <cit.> provides a good
implementation for solving this projection problem, using a quadratic
programming technique.
 



§ A COMPARISON BY SIMULATION


In this section, we shall assess the performance of the CCMLE via its Mean
Squared Error, and the bootstrap confidence intervals constructed using it. We
take the opportunity to highlight that the errors and intervals have to respect
the selection procedure.

For instance, suppose that populations A, B and C have true means 3, 2 and 1
respectively. If the corresponding sample means are 2.1, 2.2 and 1.8, then the
population selected as the maximum would be population B. The error in
estimating the mean of the selected population using the sample mean would be 

    2 - 2.2 = -0.2

This is the methodology employed in Section <ref>. The error should
not be computed as 3 - 2.2 = +0.8, since the true mean of the selected
population is in fact 2.

Similarly, when bootstrapping the strata in Section <ref>, the
population selected to have the maximum will not always be population B. It
could be population A or even population C depending on the bootstrap sample
drawn.  Thus it is not accurate to describe it as a confidence interval for
population A (which has the largest mean). It is a confidence interval for
the population selected to have the maximum mean.



 §.§ Mean Squared Error Comparison


In this subsection, we conduct a simulation study to understand the MSE of the
CCMLE, as compared to the MSE of the ordinary MLE. We consider only the cases
when p=2 and p=3 as they are sufficiently informative. 

For the p=2 case, we fix μ_2 = 0, and vary μ_1 from 0 to 5. For each
configuration, we generate 1000 bivariate N(μ, I) random vectors,
and then estimate the mean of the population with the larger sample mean.

The MSE estimate for each configuration has been plotted and smoothed in Figure
<ref>. Notice that the CCMLE performs very well when μ_1 - μ_2
≤ 1.5. Beyond that, it performs approximately 10% worse than the
unconditional MLE until the difference in means becomes quite large. From such
a point onwards, it will be the case that the sample means will be far enough
apart to warrant no shrinkage at all. It seems reasonable to guess that when
populations are close together, the CCMLE will be very beneficial; however, when
the population means are in fact far apart, it is not the right estimator to
use. We shall witness this in the p=3 case as well.



Now let us turn to the situation when p=3. In our experiment, we considered 3
possible values for μ_3:0, 2 and 4. For each μ_3, we varied μ_2
and μ_1 between μ_3 and 5. The ouput for these experiments is shown in 
Figure <ref>. 



Consider the three separate 2-by-3 lattice plots in Figure <ref>. In
each, the level plots in the top row correspond to the ordinary MLE, and the
plots in the second row correspond to the CCMLE. The plots in the left-most
column correspond to the errors when estimating the mean of the population with
the maximum sample mean, and the plots in the right-most column
correspond to the errors for the population with the minimum sample mean.
Note that colors closer to black indicate a good estimator (low MSE) while
colors closer to white indicate a poorer performance (high MSE).

Let us focus first on the case μ_3=0. Observe that there is a diagonal dark
line for the CCMLE in the left-most column. In the configurations on the
diagonal, the means of population 1 and 2 are close to each other, and hence the
CCMLE does well. In the off-diagonals, the populations are better separated and
hence performance goes down. Overall though, it does appear that the MLE and
CCMLE have a similar performance for this setting.


In the second configuration for p=3, we had fixed μ_3 to be 2, while
μ_1 and μ_2 varied from 2 to 5. The level plots for this configuration
dark regions for the CCMLE than the MLE, when compared to the case when
μ_3=0. The reason is that the number of configurations where at least two of
the populations are close together has increased, resulting in overall better
performance from the CCMLE.








Now focus on the displayed plots for μ_3=4, in Figure <ref>.
The difference between the CCMLE and the MLE is more pronounced when we consider
this final configuration, where μ_1 and μ_2 vary between 4 and 5. The
regions in the level plot for the CCMLE are consistently darker than those for
the MLE, due to the close proximity of sample means generated.



 §.§ Confidence Intervals based on the CCMLE



In this subsection, we use the stratified bootstrap to assess the confidence
intervals from the CCMLE procedure. Consider p=3, and the true means to be
μ_1, μ_2, μ_3. We simulate a single sample of size 50 from each
population. We draw from N(μ_i, √(50)) so that the sample mean still has
variance 1. Then we draw 9999 bootstrap samples from within each sample and
repeatedly compute the CCMLE. Each time, we are returned with an estimate of the
mean of the populations selected to be the maximum, middle and minimum. The bias
corrected intervals are then plotted for comparison with the traditional
intervals that do not incorporate the selection process. This procedure is
carried out for 4 different configurations of true means. The output can be seen
in Figure <ref>.

The output is most informative for the configuration where μ_1=10, μ_2 =
9.5 and μ_3 = 9. In this case, the CCMLE point estimates are all equal,
since there is insufficient power to discriminate between them. The CCMLE intervals
are roughly the same, though it is worth pointing out the slight shrinkage, in
opposite directions, of the intervals for the maximum and the minimum. In the
situations where the mean of a group is far from the rest (for instance, in the
bars for “Max” in the plots on the right, and the bars for “Min” in the
bottom left), the CCMLE and the traditional approach provide comparable
intervals. In cases where means are not distinguishable, (for instance, in the
bars for “Mid” and “Min” in the top right), the CCMLE suggests it can provide 
shorter intervals. 





§ DISCUSSION


In this paper, we introduced a new estimator for the means of selected
populations and have started to unveil some of its properties. Although
simulations experiments suggest the estimator is not admissible (see Figures 2
and 3), it performs well, particularly when the population means are
close together. Furthermore, the proposed CCMLE provides simultaneous inference
on selected populations, without a need to pre-specify the number of populations
selected. Another advantage of the procedure is that because it is frequentist
in paradigm, there is no need for prior specification.

Although the focus in this paper has been on selection via ranking of sample
means, conceptually, this approach allows for any other selection criterion to
be used. For instance, if populations were to be selected based on the absolute
values of the sample means, the primary modification would be to the probability
in the denominator of equation (<ref>). Finally, as presented, it is
straightforward to obtain bootstrap confidence intervals for this estimator, as
demonstrated in Section <ref>.

Future work include the study of the asymptotic properties of this
estimator, and computationally efficient methods to approximate the probability
P_μ(X_1>…>X_p) in a general framework.

plain

