# The MicroBooNE Experiment, the NuMI Absorber, and Heavy Neutral Leptons

Kevin J. Kelly and Pedro A.N. Machado  Theory Department, Fermi National
Accelerator Laboratory, Batavia, IL 60510, USA

###### Abstract

Motivated by the recent search for a Higgs-Portal Scalar decaying inside the
MicroBooNE detector, we demonstrate that the same search can be used to
constrain Heavy Neutral Leptons. These are gauge-singlet fermions that
interact with the Standard Model by mixing with neutrinos only and could be
related to the origin of neutrino masses. By recasting the results of the
MicroBooNE Collaboration’s analysis, we show that, for a Heavy Neutral Lepton
that mixes predominantly with muon-flavored neutrinos, previously unexplored
parameter space can be excluded for masses between 30 and 150 MeV.
Additionally, we make our Monte Carlo tools publicly available.

††preprint: FERMILAB-PUB-21-277-T

Introduction – The revelation that neutrinos are massive particles has spurred
great interest in their properties. Perhaps the simplest and most elegant mass
mechanisms for neutrinos are based on the seesaw concept Mohapatra and Valle
(1986); Gell-Mann  _et al._ (1979); Yanagida (1979); Mohapatra and Senjanovic
(1980); Schechter and Valle (1980); Magg and Wetterich (1980); Mohapatra and
Senjanovic (1980, 1981); Lazarides  _et al._ (1981); Ma and Sarkar (1998);
Foot  _et al._ (1989); Ma (1998), in which the scale of neutrino mass is
suppressed with respect to the electroweak scale due to the presence of new
particles and the breaking of lepton number. In particular, several seesaw
models postulate the existence of new leptons, neutral under the standard
model gauge group –the heavy neutral leptons (HNLs).

In this letter, we show that current data from the MicroBooNE experiment
Abratenko  _et al._ (2021) can constrain the existence of HNLs in an
unexplored region of parameter space for masses between ${\sim}20-200$ MeV. In
general, HNLs designate any fermion, neutral under the standard model gauge
group, that mixes with neutrinos. Although HNLs are a common, generic outcome
of many seesaw models, there is no clear indication of what are their masses,
neither how large are their mixings with active neutrino flavors. This
unpredictability has inspired a large effort to search for HNLs in a variety
of environments: kinematics of meson decays Bernardi  _et al._ (1988);
Aguilar-Arevalo  _et al._ (2015, 2018); Sadovsky  _et al._ (2018); Lazzeroni
_et al._ (2017); Cortina Gil  _et al._ (2018); Aguilar-Arevalo  _et al._
(2019); Cortina Gil  _et al._ (2020, 2021); production and decay in beam dump
and neutrino experiments Gorbunov and Shaposhnikov (2007); Artamonov  _et al._
(2015); Alekhin  _et al._ (2016); Ballett  _et al._ (2017, 2020); Abe  _et
al._ (2019); Berryman  _et al._ (2020); de Gouvêa  _et al._ (2021); beta-decay
spectral distortions Schreckenbach  _et al._ (1983); Deutsch  _et al._ (1990);
Hiddemann  _et al._ (1995); Derbin  _et al._ (1997); Holzschuh  _et al._
(1999, 2000); Galeazzi  _et al._ (2001); Belesev  _et al._ (2014); direct
production at high energy colliders Bergsma  _et al._ (1986); Cooper-Sarkar
_et al._ (1985); Abreu  _et al._ (1997); Aaboud  _et al._ (2019); Sirunyan
_et al._ (2019); Aaij  _et al._ (2021); Blondel  _et al._ (2021); and so on
(see, e.g., Refs Atre  _et al._ (2009); de Gouvêa and Kobach (2016); Bryman
and Shrock (2019) for compilations of constraints). In this work, we are
interested in recasting the recent MicroBooNE Higgs-Portal Scalar (HPS)
analysis Abratenko  _et al._ (2021) to constrain HNLs in the 20-200 MeV mass
range, as exhibited by Fig. 1. We include a github repository with simulation
code for generating signal events for both HNL and HPS decays in this setup
Kelly and Machado (2021).

Before going into details of our recasting, it is useful to describe the
aforementioned search by the MicroBooNE experiment, a liquid argon time
projection chamber detector in the Booster Neutrino Beam line Antonello  _et
al._ (2015). In Higgs portal models Patt and Wilczek (2006); Batell  _et al._
(2019); Egana-Ugrinovic  _et al._ (2020); Dev  _et al._ (2020); Archer-Smith
and Zhang (2021); Foroughi-Abari and Ritz (2020), a new, light scalar $S$
mixes with the Higgs, acquiring small couplings to all standard model
fermions. This mixing leads to production of the light scalars via loop
effects in meson decays such as $K^{+}\to\pi^{+}S$, as well as decay of the
light scalar into standard model fermions. The MicroBooNE analysis searched
for $S\to e^{+}e^{-}$ events which would be the dominant signature of the
decay of this light scalar if its mass is below the dimuon production
threshold.

Figure 1: MicroBooNE constraint on heavy neutral lepton parameter space as a
function of its mass and mixing with muon neutrinos. Two assumptions regarding
signal efficiency are shown in black (see text for detail) as well as
comparisons against existing constraints Aguilar-Arevalo  _et al._ (2019);
Bryman and Shrock (2019); de Gouvêa and Kobach (2016); Hayano  _et al._
(1982); Bernardi  _et al._ (1988); Artamonov  _et al._ (2015) (colored regions
as labelled).

The main experimental difficulty here lies in the fact that neutrino
interactions in argon that lead to electromagnetic showers, such as $\pi^{0}$
and photon production or $\nu_{e}$ charged current interactions, are a
background to this search. In view of that, MicroBooNE uses a clever,
unorthodox strategy: to look for events in phase with the NuMI beam which
points $8^{\circ}$ away from MicroBooNE. In NuMI, charged kaons are produced
when protons hit a target, and then travel through a decay pipe. Most $K^{+}$
that do not decay will reach the absorber, stop, and decay at rest. The
absorber is located 100 meters from the MicroBooNE detector and at an angle of
about $125^{\circ}$ with respect to the Booster Neutrino Beam line. The
$8^{\circ}$ off-axis angle between the NuMI beam and MicroBooNE reduces the
neutrino flux from NuMI at MicroBooNE. Moreover, the large angle between the
NuMI absorber and the detector relative to the NuMI beam direction can be used
to further mitigate backgrounds. These two features allows for a competitive
experimental sensitivity to Higgs Portal Scalars (HPS).

MicroBooNE exploits several kinematical features of the $e^{+}e^{-}$ pair to
increase the signal-to-background ratio. These include the opening angle
between the electrons and their angles with respect to the line connecting the
absorber and the detector center111The MicroBooNE detector does not
distinguish electrons from positrons, so we will use “electron” whenever we
refer to the $e^{+}$ or $e^{-}$ experimental signature., and the number of
hits and length of the reconstructed objects in the LArTPC. No attempt to
perform particle identification or energy reconstruction is made.

The key to our recast lies in the fact that HNLs would be produced in the NuMI
absorber and detected at MicroBooNE by processes very similar to those of the
HPS scenario. Stopped kaons could decay isotropically to HNLs via $\nu_{\mu}$
mixing, namely $K^{+}\to\mu^{+}N$. Then, these HNLs could reach MicroBooNE and
decay inside the detector. In particular, the neutral current decay
$N\to\nu_{\mu}e^{+}e^{-}$, via the same mixing, would lead to $e^{+}e^{-}$
pairs in MicroBooNE which would have similar kinematic features to the HPS
scenario. The very same data used to constrain HPS can be readily used to put
a competitive constraint on as-yet unexplored values of HNL masses and
mixings.

Heavy neutral lepton model – As discussed, there is no clear-cut indication of
what would be the HNL masses and mixings, or even how many of them could exist
in nature. In light of that, we take a typical, simplified model approach
which has the advantage of being general. We consider a single HNL which mixes
significantly only with muon neutrinos, that is

$|\nu_{\mu}\rangle=\sum_{i=1}^{3}U_{\mu i}^{*}|\nu_{i}\rangle+U_{\mu
4}^{*}|N\rangle,$ (1)

where $U$ is a generalized version of the usual Pontecorvo-Maki-Nakagawa-
Sakata matrix and $N$ denotes the HNL. Nonzero $U_{\tau 4}$ would not change
the phenomenology we are interested in. Although mixing with electron neutrino
could lead to the same experimental signatures we are focusing on, existing
constraints from $\pi^{\pm}$ decays (for $m_{N}<m_{\pi^{\pm}}-m_{e}$) are
stronger than those on $U_{\mu 4}$ for similar masses, which would make the
MicroBooNE constraint on $U_{e4}$ not competitive in this analysis, as checked
numerically.

A nonzero mixing with muon neutrinos induces kaon branching ratio to $N$

$\text{Br}(K\\!\\!\to\\!\mu N)\simeq\text{Br}(K\\!\\!\to\\!\mu\nu)|U_{\mu
4}|^{2}\rho_{N}\\!\\!\left(\frac{m_{\mu}^{2}}{m_{K}^{2}},\,\frac{m_{N}^{2}}{m_{K}^{2}}\right),$
(2)

where we have approximated $|U_{\mu 4}|^{2}\ll 1$, which is justified by
existing experimental bounds; $m_{\mu}$, $m_{K}$ and $m_{N}$ are the masses of
the muon, the kaon and the HNL; and

$\rho_{N}(x,y)=\frac{[x+y-(x-y)]^{2}\sqrt{1+(x-y)^{2}-2(x+y)}}{x(1-x)^{2}}$

is a phase space factor such that in the limit $m_{\mu,N}\ll m_{K}$ we have
$\rho_{N}\to 1$. The mixing with muon neutrinos also controls the $N$ decay
rate to $e^{+}e^{-}$. This decay depends on the nature of the HNL (Majorana
versus Dirac). Assuming for concreteness that $N$ is a Majorana fermion, its
partial decay width to $e^{+}e^{-}$ pairs is given by Gorbunov and
Shaposhnikov (2007)

$\displaystyle\Gamma(N\to\nu_{\mu}e^{+}e^{-})$
$\displaystyle=2\frac{G_{F}^{2}|U_{\mu
4}|^{2}m_{N}^{5}}{768\pi^{3}}\left(1-4s_{W}^{2}+8s_{W}^{4}\right)$
$\displaystyle\hskip 76.82234pt+\mathcal{O}(m_{e}^{2}/m_{N}^{2}),$ (3)

where terms proportional to the ratio of the electron to HNL mass are
neglected and $s_{W}$ is the sine of the weak mixing angle. We have also kept
explicit the factor of 2 due to the Majorana nature of $N$. This decay, for
$20<m_{N}<200$ MeV corresponds to a ${\sim}10\%$ branching ratio of all $N$
decays. The other $90\%$ is mostly $N\to\nu\nu\nu$, which is unobservable.

Signal Production and Constraint Recasting – The event rate $R_{X}$ of a new-
physics particle $X$, produced in a two-body kaon decay-at-rest, traveling to
MicroBooNE from the NuMI absorber, and decaying within, can be expressed as
the product of the production flux times the probability of decay:

$\displaystyle R_{X}$ $\displaystyle=\Phi_{X}A\,P(X\\!\to
e^{+}\\!e^{-})\varepsilon(m_{X}),$ (4) $\displaystyle\Phi_{X}$
$\displaystyle=\frac{N_{K\mathrm{DAR}}\mathrm{Br}\left(K^{\pm}\to
X\right)}{4\pi D^{2}},$ (5) $\displaystyle P(X\\!\to e^{+}\\!e^{-})$
$\displaystyle=\frac{1}{\gamma c\tau_{X}}\int_{D}^{D+L}e^{-\frac{z}{\gamma
c\tau_{X}}}dz$ $\displaystyle\quad\times\mathrm{Br}(X\\!\to e^{+}\\!e^{-})$
(6)

Here, $\Phi_{X}$ is the flux of $X$ at the MicroBooNE detector, assumed to be
nearly constant over the extent of the detector; $A$ is the cross-sectional
area as viewed by an incoming $X$ from the NuMI absorber; $P(X\\!\to
e^{+}\\!e^{-})$ is the probability of $X$ decaying within the detector volume
to, in our case, $e^{+}e^{-}+\,$anything; and $\varepsilon(m_{X})$ is the
signal reconstruction efficiency, which may depend on the mass of the decaying
particle. $N_{K\mathrm{DAR}}$ is the number of kaons decaying at rest in the
NuMI absorber during the MicroBooNE data collection, $D$ is the
absorber/detector distance (${\sim}100$ m), $L$ is the extent of the detector
for an incoming $X$, $\gamma=E_{X}/m_{X}$ can be determined from the two-body
production of $X$, and $\tau_{X}$ is the proper lifetime of the particle $X$.

If $X$ is long-lived ($c\tau_{X}\gg D,~{}L$) relative to the other relevant
distance scales222This limit is well-satisfied for $\left\lvert U_{\mu
4}\right\rvert^{2}<(m_{N}/100\ \mathrm{MeV})^{-5}$. For the HPS, the
requirement is $\sin^{2}\theta_{S}<10^{-5}\ (m_{S}/100\ \mathrm{MeV})^{-1}$,
satisfied for the parameter space of interest in Ref. Abratenko  _et al._
(2021),

$P(X\\!\to e^{+}\\!e^{-})\simeq\frac{L}{\gamma
c\tau_{X}}\mathrm{Br}\left(X\\!\to
e^{+}\\!e^{-}\right)=\frac{L}{\gamma}\Gamma(X\\!\to e^{+}\\!e^{-}).$ (7)

When multiplied with $A$ in Eq. (4), only the total volume of the detector,
not its specific shape, enters the rate. We can compare the relative rates of
hypothetical HPS and HNL decays within MicroBooNE coming from the NuMI
absorber by calculating $R_{N}/R_{S}$:

$\displaystyle\frac{R_{N}}{R_{S}}=\frac{\Phi_{N}}{\Phi_{S}}\frac{P(N\to\nu
e^{+}e^{-})}{P(S\to
e^{+}e^{-})}\frac{\varepsilon(m_{N})}{\varepsilon(m_{S})},$ (8)
$\displaystyle=\frac{\mathrm{Br}\left(K^{\pm}\to\mu^{\pm}N\right)}{\mathrm{Br}\left(K^{\pm}\to\pi^{\pm}S\right)}\frac{m_{N}E_{S}\Gamma\left(N\to\nu
e^{+}e^{-}\right)}{m_{S}E_{N}\Gamma\left(S\to
e^{+}e^{-}\right)}\frac{\varepsilon(m_{N})}{\varepsilon(m_{S})},$

where $E_{N,S}$ are the energies of the HNL and HPS. This relative rate is
described by the fluxes, probabilities of decaying within the detector, and
signal efficiencies of the respective models. We defer discussion of the
relative efficiencies to the next section.

The first two fractions in Eq. (8) can be determined from the HNL and HPS
model parameters — each is proportional to $\left\lvert U_{\mu
N}\right\rvert^{2}/\sin^{2}\theta_{S}$. In Fig. 2, we present these two
fractions (the relative fluxes in blue and the relative decay probabilities in
red) as a function of the HNL or HPS mass, with benchmark values of
$\left\lvert U_{\mu N}\right\rvert^{2}=10^{-4}$ and
$\sin^{2}\theta_{S}=10^{-6}$.

Figure 2: Relative flux (blue) and decay probabilities (red) between the heavy
neutral lepton $N$ and the Higgs portal scalar $S$ coming from the NuMI
absorber via two-body kaon decay-at-rest and decaying within MicroBooNE. The
product of these two, in addition to relative signal efficiencies, provides
the relative signal event rate in MicroBooNE.

We see that, for these benchmark mixing values, a much larger flux of $N$ is
produced relative to $S$, however, $S$ decays333The sharp mass-dependence of
the red curve in Fig. 2 comes from the relative scaling with masses of
$\Gamma(N\to\nu e^{+}e^{-})\propto m_{N}^{5}$ vs. $\Gamma(S\to
e^{+}e^{-})\propto m_{S}$. much more rapidly than $N$ (effectively,
$c\tau_{S}\gg c\tau_{N}\gg D$), leading to a much smaller probability that a
given $N$ decays into an electron/positron pair in MicroBooNE than a given $S$
does. This combination implies that the HNL and HPS event rates can be
comparable if one assumes similar signal reconstruction efficiencies. Ref.
Abratenko  _et al._ (2021) provides constraints on $\theta_{S}$ and the
efficiency $\varepsilon(m_{S})$ for $m_{S}$ between $2m_{e}$ and $2m_{\mu}$:
with reasonable assumptions about $\varepsilon(m_{N})$, we can recast
MicroBooNE’s HPS constraint onto HNL parameter space.

Event Kinematics & Signal Efficiency – The key challenge in recasting
MicroBooNE’s HPS search onto HNL parameter space is in determining the signal
efficiency $\varepsilon(m_{N})$. We take two approaches: first, we assume that
$\varepsilon(m_{N})=\varepsilon(m_{S})$ and that events from an $N$ of a
particular mass are accepted at the same rate as those from $S$ with the same
mass. There are obvious shortcomings with this approach, however: because
$N$’s decay is three-body, the kinematics of the $e^{+}e^{-}$ pair are
different than those from the two-body $S\to e^{+}e^{-}$ decay. In Ref.
Abratenko  _et al._ (2021), signal identification (and background reduction)
is optimized using boosted decision trees with kinematic quantities,
especially the opening angle of the $e^{+}e^{-}$ pair and their angles
relative to the direction from the NuMI absorber: no particle ID or
calorimetry is used. Therefore, it seems more reasonable that the efficiency
depends more strongly on the opening angle of the $e^{+}e^{-}$ pair than on
the masses of the decaying particles.

To obtain a more realistic efficiency for HNL, we adopt the following second
approach. For each HPS mass, we map its efficiency onto the average cosine of
the $e^{+}e^{-}$ opening angle $\overline{\cos\theta_{e^{+}e^{-}}}$. Then we
assign the HNL efficiency to match the HPS efficiency for the same value of
$\overline{\cos\theta_{e^{+}e^{-}}}$. Simulation code for three-body HNL decay
$N\to\nu e^{+}e^{-}$, including alternate model assumptions such as $N$ being
a Dirac fermion, is available at this URL Kelly and Machado (2021).

To justify this approach, Fig. 3(top) displays several truth-level signal
distributions as a function of the opening angle between the electron/positron
pair $\cos\theta_{e^{+}e^{-}}$ and the angle between the higher-energy
electron and the absorber direction $\theta_{e_{\rm lead}}$. Three
distributions are shown: the two purple/orange ones correspond to HPS with
$m_{S}=100$ and $80$ MeV, as labelled, and the blue/green one corresponds to
an HNL with $m_{N}=100$ MeV. For the HPS ones, we note that there is only one
relevant kinematical quantity in the two-body $S$ decay, the rest-frame angle
of the decay. This produces a strong correlation between these two quantities
in the lab frame. In all three distributions, brighter colors correspond to
where more events are expected.

For comparison, the bottom panel of Fig. 3 displays the same information after
applying a $3\degree$ angular uncertainty to each of the electromagnetic
showers Adams  _et al._ (2020). While we see a considerable overlap for
$m_{N}=m_{S}$, there is even more overlap between $m_{S}=80$ MeV and
$m_{N}=100$ MeV, as three-body decays prefer larger $\cos\theta_{e^{+}e^{-}}$
for the same mass of the mother particle. Note that $m_{S}=80$ MeV and
$m_{N}=100$ MeV exhibit approximately the same value of
$\overline{\cos\theta_{e^{+}e^{-}}}$. This overlap is especially obvious in
the bottom panel in terms of the reconstructed angles.

Figure 3: Distributions of $e^{+}e^{-}$ pairs as a function of the cosine of
their opening angle vs. the angle of the higher-energy particle relative to
the NuMI absorber direction. Top: truth-level information, Bottom:
reconstructed information after applying a 3$\degree$ angular uncertainty on
electron tracks. Blue/green points come from an HNL with $m_{N}=100$ MeV where
the two purple/orange ones come from HPS with $m_{S}=100$ MeV (top) and
$m_{S}=80$ MeV (bottom). In all three, brighter/darker colors correspond to
more signal events.

Finally, since HNLs decay to three daughter particles, the decay products tend
to be less energetic. As the reconstruction of low energy electrons may be
challenging, we also apply two kinematical cuts on top of both HPS and HNL
efficiencies: $E_{e^{\pm}}>10$ MeV, so that the two tracks can be
reconstructed; and
$\left\lvert\cos\theta_{e^{+}e^{-}}\right\rvert<\cos{(10\degree)}$, so that
two distinct tracks can be identified. Fig. 4 shows the HNL efficiency as a
function of the HNL mass for the two benchmark assumptions:
$\varepsilon(m_{S})=\varepsilon(m_{N})$ (solid line) and $\varepsilon(m_{N})$
tracing the average cosine of the $e^{+}e^{-}$ opening angle (dashed line).

Figure 4: Our two assumed signal efficiencies for HNL reconstruction: that the
HNL efficiency matches the HPS efficiency given in Ref. Abratenko  _et al._
(2021) (solid) and that the dependence of the efficiency is based on the
average opening angle of the electron/positron pair (dashed).

Results and discussion – Under these two assumptions for $\varepsilon(m_{N})$,
in tandem with the efficiencies and limits on $\theta_{S}$ provided in Ref.
Abratenko  _et al._ (2021), we can use Eq. (8) to determine the parameter
space of $m_{N}$ vs. $\left\lvert U_{\mu 4}\right\rvert^{2}$ for which a
comparable signal rate is expected and therefore, MicroBooNE can set a
constraint. This is shown in Fig. 1 with the solid line corresponding to
$\varepsilon(m_{N})=\varepsilon(m_{S})$ and the dashed one to
$\varepsilon(m_{N})=\varepsilon(\overline{\cos\theta_{e^{+}e^{-}}})$. Colored
regions depict existing experimental constraints as labelled from Refs.
Aguilar-Arevalo  _et al._ (2019); Bryman and Shrock (2019); de Gouvêa and
Kobach (2016); Hayano  _et al._ (1982); Bernardi  _et al._ (1988); Artamonov
_et al._ (2015). We see that MicroBooNE can constrain new parameter space for
$34$ MeV $<m_{N}<150$ MeV, including by over two orders of magnitude for
$m_{N}\approx 65$ MeV. Even under the more conservative assumption on the HNL
signal efficiency, the MicroBooNE constraint using current data remains world-
leading in the aforementioned HNL mass window.

Conclusions – We have recasted MicroBooNE’s search for Higgs portal scalars
onto the heavy neutral lepton parameter space. We have shown that, even with
conservative assumptions on the HNL signal efficiency, MicroBooNE data still
sets the world leading constraints on HNLs with masses between 34 and 150 MeV.
We encourage the MicroBooNE collaboration to use our publicly available Monte
Carlo simulation code to determine $\varepsilon(m_{N})$ more exactly and carry
out this analysis. Finally, we call attention that a very similar search could
be performed at ICARUS. While its distance to the NuMI absorber is also about
100 meters, its volume is about 5 times larger than that of MicroBooNE.

###### Acknowledgements.

Acknowledgements: We thank Patrick Fox for valuable conversations on this work
and Nikita Blinov and Shirley Li for testing our simulation code. Fermilab is
operated by the Fermi Research Alliance, LLC under contract No. DE-
AC02-07CH11359 with the United States Department of Energy. This project has
received support from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie grant agreement No.
860881-HIDDeN.

## References

  * Mohapatra and Valle (1986) R. N. Mohapatra and J. W. F. Valle, “Neutrino Mass and Baryon Number Nonconservation in Superstring Models,” Phys. Rev. D 34, 1642 (1986).
  * Gell-Mann  _et al._ (1979) Murray Gell-Mann, Pierre Ramond, and Richard Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927, 315–321 (1979), arXiv:1306.4669 [hep-th] .
  * Yanagida (1979) Tsutomu Yanagida, “Horizontal gauge symmetry and masses of neutrinos,” Conf. Proc. C 7902131, 95–99 (1979).
  * Mohapatra and Senjanovic (1980) Rabindra N. Mohapatra and Goran Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44, 912 (1980).
  * Schechter and Valle (1980) J. Schechter and J. W. F. Valle, “Neutrino Masses in SU(2) x U(1) Theories,” Phys. Rev. D 22, 2227 (1980).
  * Magg and Wetterich (1980) M. Magg and C. Wetterich, “Neutrino Mass Problem and Gauge Hierarchy,” Phys. Lett. B 94, 61–64 (1980).
  * Mohapatra and Senjanovic (1981) Rabindra N. Mohapatra and Goran Senjanovic, “Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation,” Phys. Rev. D 23, 165 (1981).
  * Lazarides  _et al._ (1981) George Lazarides, Q. Shafi, and C. Wetterich, “Proton Lifetime and Fermion Masses in an SO(10) Model,” Nucl. Phys. B 181, 287–300 (1981).
  * Ma and Sarkar (1998) Ernest Ma and Utpal Sarkar, “Neutrino masses and leptogenesis with heavy Higgs triplets,” Phys. Rev. Lett. 80, 5716–5719 (1998), arXiv:hep-ph/9802445 .
  * Foot  _et al._ (1989) Robert Foot, H. Lew, X. G. He, and Girish C. Joshi, “Seesaw Neutrino Masses Induced by a Triplet of Leptons,” Z. Phys. C 44, 441 (1989).
  * Ma (1998) Ernest Ma, “Pathways to naturally small neutrino masses,” Phys. Rev. Lett. 81, 1171–1174 (1998), arXiv:hep-ph/9805219 .
  * Abratenko  _et al._ (2021) P. Abratenko _et al._ (MicroBooNE), “Search for a Higgs portal scalar decaying to electron-positron pairs in the MicroBooNE detector,” (2021), arXiv:2106.00568 [hep-ex] .
  * Bernardi  _et al._ (1988) G. Bernardi _et al._ , “FURTHER LIMITS ON HEAVY NEUTRINO COUPLINGS,” Phys. Lett. B 203, 332–334 (1988).
  * Aguilar-Arevalo  _et al._ (2015) A. Aguilar-Arevalo _et al._ (PiENu), “Improved Measurement of the $\pi\to\textrm{e}\nu$ Branching Ratio,” Phys. Rev. Lett. 115, 071801 (2015), arXiv:1506.05845 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2018) A. Aguilar-Arevalo _et al._ (PIENU), “Improved search for heavy neutrinos in the decay $\pi\rightarrow e\nu$,” Phys. Rev. D 97, 072012 (2018), arXiv:1712.03275 [hep-ex] .
  * Sadovsky  _et al._ (2018) A. S. Sadovsky _et al._ (OKA), “Search for heavy neutrino in $K^{+}\rightarrow\mu^{+}\nu_{H}$ decay,” Eur. Phys. J. C 78, 92 (2018), arXiv:1709.01473 [hep-ex] .
  * Lazzeroni  _et al._ (2017) Cristina Lazzeroni _et al._ (NA62), “Search for heavy neutrinos in $K^{+}\rightarrow\mu^{+}\nu_{\mu}$ decays,” Phys. Lett. B 772, 712–718 (2017), arXiv:1705.07510 [hep-ex] .
  * Cortina Gil  _et al._ (2018) Eduardo Cortina Gil _et al._ (NA62), “Search for heavy neutral lepton production in $K^{+}$ decays,” Phys. Lett. B 778, 137–145 (2018), arXiv:1712.00297 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2019) A. Aguilar-Arevalo _et al._ (PIENU), “Search for heavy neutrinos in $\pi\to\mu\nu$ decay,” Phys. Lett. B 798, 134980 (2019), arXiv:1904.03269 [hep-ex] .
  * Cortina Gil  _et al._ (2020) Eduardo Cortina Gil _et al._ (NA62), “Search for heavy neutral lepton production in $K^{+}$ decays to positrons,” Phys. Lett. B 807, 135599 (2020), arXiv:2005.09575 [hep-ex] .
  * Cortina Gil  _et al._ (2021) Eduardo Cortina Gil _et al._ (NA62), “Search for $K^{+}$ decays to a muon and invisible particles,” Phys. Lett. B 816, 136259 (2021), arXiv:2101.12304 [hep-ex] .
  * Gorbunov and Shaposhnikov (2007) Dmitry Gorbunov and Mikhail Shaposhnikov, “How to find neutral leptons of the $\nu$MSM?” JHEP 10, 015 (2007), [Erratum: JHEP 11, 101 (2013)], arXiv:0705.1729 [hep-ph] .
  * Artamonov  _et al._ (2015) A. V. Artamonov _et al._ (E949), “Search for heavy neutrinos in $K^{+}\to\mu^{+}\nu_{H}$ decays,” Phys. Rev. D 91, 052001 (2015), [Erratum: Phys.Rev.D 91, 059903 (2015)], arXiv:1411.3963 [hep-ex] .
  * Alekhin  _et al._ (2016) Sergey Alekhin _et al._ , “A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case,” Rept. Prog. Phys. 79, 124201 (2016), arXiv:1504.04855 [hep-ph] .
  * Ballett  _et al._ (2017) Peter Ballett, Silvia Pascoli, and Mark Ross-Lonergan, “MeV-scale sterile neutrino decays at the Fermilab Short-Baseline Neutrino program,” JHEP 04, 102 (2017), arXiv:1610.08512 [hep-ph] .
  * Ballett  _et al._ (2020) Peter Ballett, Tommaso Boschi, and Silvia Pascoli, “Heavy Neutral Leptons from low-scale seesaws at the DUNE Near Detector,” JHEP 03, 111 (2020), arXiv:1905.00284 [hep-ph] .
  * Abe  _et al._ (2019) K. Abe _et al._ (T2K), “Search for heavy neutrinos with the T2K near detector ND280,” Phys. Rev. D 100, 052006 (2019), arXiv:1902.07598 [hep-ex] .
  * Berryman  _et al._ (2020) Jeffrey M. Berryman, Andre de Gouvea, Patrick J Fox, Boris J Kayser, Kevin J Kelly, and Jennifer L Raaf, “Searches for Decays of New Particles in the DUNE Multi-Purpose Near Detector,” JHEP 02, 174 (2020), arXiv:1912.07622 [hep-ph] .
  * de Gouvêa  _et al._ (2021) André de Gouvêa, Patrick J Fox, Boris J Kayser, and Kevin J. Kelly, “Three-Body Decays of Heavy Dirac and Majorana Fermions,” (2021), arXiv:2104.05719 [hep-ph] .
  * Schreckenbach  _et al._ (1983) K. Schreckenbach, G. Colvin, and F. Von Feilitzsch, “SEARCH FOR MIXING OF HEAVY NEUTRINOS IN THE BETA+ AND BETA- SPECTRA OF THE CU-64 DECAY,” Phys. Lett. B 129, 265–268 (1983).
  * Deutsch  _et al._ (1990) J. Deutsch, M. Lebrun, and R. Prieels, “Searches for admixture of massive neutrinos into the electron flavor,” Nucl. Phys. A 518, 149–155 (1990).
  * Hiddemann  _et al._ (1995) K. H. Hiddemann, H. Daniel, and O. Schwentker, “Limits on neutrino masses from the tritium beta spectrum,” J. Phys. G 21, 639–650 (1995).
  * Derbin  _et al._ (1997) A. V. Derbin, A. I. Egorov, S. V. Bakhlanov, and V. N. Muratova, “Measurement of the Ca-45 beta spectrum in search of deviations from the theoretical shape,” JETP Lett. 66, 88–92 (1997).
  * Holzschuh  _et al._ (1999) E. Holzschuh, W. Kundig, L. Palermo, H. Stussi, and P. Wenk, “Search for heavy neutrinos in the beta spectrum of Ni-63,” Phys. Lett. B 451, 247–255 (1999).
  * Holzschuh  _et al._ (2000) E. Holzschuh, L. Palermo, H. Stussi, and P. Wenk, “The beta-spectrum of S-35 and search for the admixture of heavy neutrinos,” Phys. Lett. B 482, 1–9 (2000).
  * Galeazzi  _et al._ (2001) M. Galeazzi, F. Fontanelli, F. Gatti, and S. Vitale, “Limits on the existence of heavy neutrinos in the range 50-eV - 1000-eV from the study of the Re-187 beta decay,” Phys. Rev. Lett. 86, 1978–1981 (2001).
  * Belesev  _et al._ (2014) A. I. Belesev, A. I. Berlev, E. V. Geraskin, A. A. Golubev, N. A. Likhovid, A. A. Nozik, V. S. Pantuev, V. I. Parfenov, and A. K. Skasyrskaya, “The search for an additional neutrino mass eigenstate in the 2–100 eV region from ‘Troitsk nu-mass’ data: a detailed analysis,” J. Phys. G 41, 015001 (2014), arXiv:1307.5687 [hep-ex] .
  * Bergsma  _et al._ (1986) F. Bergsma _et al._ (CHARM), “A Search for Decays of Heavy Neutrinos in the Mass Range 0.5-GeV to 2.8-GeV,” Phys. Lett. B 166, 473–478 (1986).
  * Cooper-Sarkar  _et al._ (1985) Amanda M. Cooper-Sarkar _et al._ (WA66), “Search for Heavy Neutrino Decays in the BEBC Beam Dump Experiment,” Phys. Lett. B 160, 207–211 (1985).
  * Abreu  _et al._ (1997) P. Abreu _et al._ (DELPHI), “Search for neutral heavy leptons produced in Z decays,” Z. Phys. C 74, 57–71 (1997), [Erratum: Z.Phys.C 75, 580 (1997)].
  * Aaboud  _et al._ (2019) Morad Aaboud _et al._ (ATLAS), “Search for a right-handed gauge boson decaying into a high-momentum heavy neutrino and a charged lepton in $pp$ collisions with the ATLAS detector at $\sqrt{s}=13$ TeV,” Phys. Lett. B 798, 134942 (2019), arXiv:1904.12679 [hep-ex] .
  * Sirunyan  _et al._ (2019) Albert M Sirunyan _et al._ (CMS), “Search for heavy Majorana neutrinos in same-sign dilepton channels in proton-proton collisions at $\sqrt{s}=13$ TeV,” JHEP 01, 122 (2019), arXiv:1806.10905 [hep-ex] .
  * Aaij  _et al._ (2021) Roel Aaij _et al._ (LHCb), “Search for heavy neutral leptons in $W^{+}\to\mu^{+}\mu^{\pm}\text{jet}$ decays,” Eur. Phys. J. C 81, 248 (2021), arXiv:2011.05263 [hep-ex] .
  * Blondel  _et al._ (2021) Alain Blondel, André de Gouvêa, and Boris Kayser, “Z-Boson Decays into Majorana or Dirac (Heavy) Neutrinos,” (2021), arXiv:2105.06576 [hep-ph] .
  * Atre  _et al._ (2009) Anupama Atre, Tao Han, Silvia Pascoli, and Bin Zhang, “The Search for Heavy Majorana Neutrinos,” JHEP 05, 030 (2009), arXiv:0901.3589 [hep-ph] .
  * de Gouvêa and Kobach (2016) André de Gouvêa and Andrew Kobach, “Global Constraints on a Heavy Neutrino,” Phys. Rev. D 93, 033005 (2016), arXiv:1511.00683 [hep-ph] .
  * Bryman and Shrock (2019) D. A. Bryman and R. Shrock, “Constraints on Sterile Neutrinos in the MeV to GeV Mass Range,” Phys. Rev. D 100, 073011 (2019), arXiv:1909.11198 [hep-ph] .
  * Kelly and Machado (2021) Kevin J. Kelly and Pedro A.N. Machado, “muBHNL on GitHub: https://github.com/kjkellyphys/muBHNL,” (2021).
  * Antonello  _et al._ (2015) M. Antonello _et al._ (MicroBooNE, LAr1-ND, ICARUS-WA104), “A Proposal for a Three Detector Short-Baseline Neutrino Oscillation Program in the Fermilab Booster Neutrino Beam,” (2015), arXiv:1503.01520 [physics.ins-det] .
  * Patt and Wilczek (2006) Brian Patt and Frank Wilczek, “Higgs-field portal into hidden sectors,” (2006), arXiv:hep-ph/0605188 .
  * Batell  _et al._ (2019) Brian Batell, Joshua Berger, and Ahmed Ismail, “Probing the Higgs Portal at the Fermilab Short-Baseline Neutrino Experiments,” Phys. Rev. D 100, 115039 (2019), arXiv:1909.11670 [hep-ph] .
  * Egana-Ugrinovic  _et al._ (2020) Daniel Egana-Ugrinovic, Samuel Homiller, and Patrick Meade, “Light Scalars and the Koto Anomaly,” Phys. Rev. Lett. 124, 191801 (2020), arXiv:1911.10203 [hep-ph] .
  * Dev  _et al._ (2020) P. S. Bhupal Dev, Rabindra N. Mohapatra, and Yongchao Zhang, “Constraints on long-lived light scalars with flavor-changing couplings and the KOTO anomaly,” Phys. Rev. D 101, 075014 (2020), arXiv:1911.12334 [hep-ph] .
  * Archer-Smith and Zhang (2021) Paul Archer-Smith and Yue Zhang, “Higgs Portal From The Atmosphere To Hyper-K,” Phys. Lett. B 817, 136309 (2021), arXiv:2005.08980 [hep-ph] .
  * Foroughi-Abari and Ritz (2020) Saeid Foroughi-Abari and Adam Ritz, “LSND Constraints on the Higgs Portal,” Phys. Rev. D 102, 035015 (2020), arXiv:2004.14515 [hep-ph] .
  * Hayano  _et al._ (1982) R. S. Hayano _et al._ , “HEAVY NEUTRINO SEARCH USING K(mu2) DECAY,” Phys. Rev. Lett. 49, 1305 (1982).
  * Adams  _et al._ (2020) C. Adams _et al._ (MicroBooNE), “Reconstruction and Measurement of $\mathcal{O}$(100) MeV Energy Electromagnetic Activity from $\pi^{0}\rightarrow\gamma\gamma$ Decays in the MicroBooNE LArTPC,” JINST 15, P02007 (2020), arXiv:1910.02166 [hep-ex] .


Design-based composite estimation of small proportions in small domains

Andrius Čiginas

Vilnius University

###### Abstract

Traditional direct estimation methods are not efficient for domains of a
survey population with small sample sizes. To estimate the domain proportions,
we combine the direct estimators and the regression-synthetic estimators based
on domain-level auxiliary information. For the case of small true proportions,
we introduce the design-based linear combination that is a robust alternative
to the empirical best linear unbiased predictor (EBLUP) based on the
Fay–Herriot model. We also consider an adaptive procedure optimizing a sample-
size-dependent composite estimator, which depends on a single parameter for
all domains.

We imitate the Lithuanian Labor Force Survey, where we estimate the
proportions of the unemployed and employed in municipalities. We show where
the considered design-based compositions and estimators of their mean square
errors are competitive for EBLUP and its accuracy estimation.

Keywords: small area estimation, area-level model, composite estimator,
sample-size-dependent estimator, Labor Force Survey.

##  1 Introduction

Design-based and model-assisted direct estimators of parameters rely only on
the sample of the estimation domain (area). Therefore, after the sample is
selected, their application for some unplanned domains leads to high variances
of the estimators because of too small sample sizes. In the small area
estimation theory [14], indirect estimators borrow sample information from
neighbor domains through auxiliary information and linking models. These
model-based estimators usually have lower variances than the direct
estimators, but their biases can be relatively large.

To estimate proportions in the domains, one can consider explicit linking
models based on auxiliary data aggregated to the domain level. A popular model
is the Fay–Herriot (FH) model, which is a separate case of linear mixed
models, and the empirical best linear unbiased predictors (EBLUPs) of the
domain means or proportions are derived from it [7]. That small area predictor
is expressed as the linear combination of a regression-synthetic estimator and
the direct estimator. While the former part accounts for a variation reflected
in the auxiliary data, the direct component exploits the unbiasedness
property. Compositions of the synthetic and the direct estimators constitute
an important class of indirect estimators. Before the mixed models,
traditional design-based composite estimators were often used [14, Chapter 3].
However, now it is accepted that the models including random area-specific
effects are more useful. For example, they are more convenient to handle
complex data structures than the traditional estimators with only randomness
induced by the sampling design. Another notable drawback of the latter
estimators is the difficulty in estimating their precision. The problem is
with bias estimation, while it is well elaborated for the estimators like
EBLUP.

We use a conditional analysis to construct the design-based composite
estimator, which is in some sense similar to EBLUP. According to the
construction, it is a robust estimator suitable for small or large domain
proportions. For the comparison, we also consider the sample-size-dependent
(SSD) compositions introduced in [6]. To optimize the latter estimators with
respect to their parameter, we apply the strategy based on the minimization of
estimated average mean square error (MSE) as proposed in [2]. The MSEs of both
the design-based compositions are estimated as suggested in [3].

We compare the estimators and their MSE estimators in the simulation study
using the Lithuanian Labor Force Survey (LFS) data, where fractions of the
unemployed and employed are the proportions of interest estimated in
municipalities. Applications of EBLUPs to LFS unemployment data of other
countries are found, for example, in [1, 8, 12]. SSD compositions, with
subjectively chosen values of the parameter, are used in [6, 18]. The adaptive
selection of values of this parameter is applied to estimate the proportions
of unemployed in [2].

##  2 Basic assumptions and direct estimation

The set ${\cal U}=\\{1,\ldots,N\\}$ consists of the labels of elements of the
survey population. Let $y$ be a binary study variable with the fixed values
$y_{1},\ldots,y_{N}$ assigned to the corresponding elements. To estimate the
proportions in the population and its subsets, the sample $s\subset{\cal U}$
of size $n<N$ is drawn by the sampling design $p(\cdot)$, and
$\pi_{k}=\operatorname{{P}_{p}}\\{k\in s\\}>0$, $k\in{\cal U}$, are inclusion
into the sample probabilities. Here the symbol $\operatorname{{P}_{p}}$, and
hereafter $\operatorname{{E}_{p}}$, $\operatorname{{var}_{p}}$, and
$\operatorname{{MSE}_{p}}$ denote probability, expectation, variance, and MSE
according to $p(\cdot)$, respectively. The characteristic
$\operatorname{{var}_{p}}(\cdot)$ is called the sampling variance or design
variance.

Let ${\cal U}={\cal U}_{1}\cup\cdots\cup{\cal U}_{M}$ be the partition of the
population into the non-overlapping domains, where the domain ${\cal U}_{i}$
contains $N_{i}$ elements. Then the domain sample $s_{i}=s\cap{\cal U}_{i}$ is
of size $n_{i}\leqslant N_{i}$. We aim to estimate the proportions

$\theta_{i}=\frac{1}{N_{i}}\sum_{k\in{\cal U}_{i}}y_{k},\qquad i=1,\ldots,M,$
(1)

where the numbers $N_{i}$ are assumed to be known. If the design $p(\cdot)$
does not ensure the fixed sizes $n_{i}$, then they can be too small to get
suficiently accurate direct estimates $\hat{\theta}_{i}^{\mathrm{d}}$ of (1).

Assume that, for each domain ${\cal U}_{i}$, the auxiliary information is
available as the vector of known characteristics
$\mathbf{z}_{i}=(z_{i1},z_{i2},\ldots,z_{iP})^{\prime}$. This assumption
narrows a choice of direct estimators to the design unbiased Horvitz–Thompson
estimators $\hat{\theta}_{i}^{\mathrm{HT}}=N_{i}^{-1}\sum_{k\in
s_{i}}y_{k}/\pi_{k}$ of $\theta_{i}$ or the weighted sample proportions

$\hat{\theta}_{i}^{\mathrm{H}}=\frac{1}{\widehat{N}_{i}}\sum_{k\in
s_{i}}\frac{y_{k}}{\pi_{k}},\quad\text{where}\quad\widehat{N}_{i}=\sum_{k\in
s_{i}}\frac{1}{\pi_{k}},\qquad i=1,\ldots,M,$ (2)

that are approximately unbiased. The approximate sampling variances of (2) and
their estimators have the expressions [17, p. 185]

$\operatorname{{var}_{p}}(\hat{\theta}_{i}^{\mathrm{H}})\approx\psi_{i}^{\mathrm{H}}=\frac{1}{N_{i}^{2}}\sum_{k\in{\cal
U}_{i}}\sum_{l\in{\cal
U}_{i}}(\pi_{kl}-\pi_{k}\pi_{l})\frac{(y_{k}-\theta_{i})(y_{l}-\theta_{i})}{\pi_{k}\pi_{l}},\qquad
i=1,\ldots,M,$ (3)

and

$\hat{\psi}_{i}^{\mathrm{H}}=\frac{1}{\widehat{N}_{i}^{2}}\sum_{k\in
s_{i}}\sum_{l\in
s_{i}}(1-\pi_{k}\pi_{l}/\pi_{kl})\frac{(y_{k}-\hat{\theta}_{i}^{\mathrm{H}})(y_{l}-\hat{\theta}_{i}^{\mathrm{H}})}{\pi_{k}\pi_{l}},\qquad
i=1,\ldots,M,$ (4)

respectively, where $\pi_{kl}=\operatorname{{P}_{p}}\\{k,l\in s\\}>0$ is the
probability that both of the elements $k$ and $l$ will be included into the
sample.

##  3 EBLUP under the Fay–Herriot model

The direct estimators $\hat{\theta}_{i}^{\mathrm{d}}$ of the domain
proportions can be improved using the FH model [7]. The data for this domain-
level model are the estimates $\hat{\theta}_{i}^{\mathrm{d}}$, their
corresponding estimates $\hat{\psi}_{i}$ of the sampling variances
$\psi_{i}=\operatorname{{var}_{p}}(\hat{\theta}_{i}^{\mathrm{d}})$, and the
covariates $\mathbf{z}_{i}$, $i=1,\ldots,M$. The basic FH model consists of
two parts, see [14, Section 4.2], that are combined into the linear mixed
model

$\hat{\theta}_{i}^{\mathrm{d}}=\mathbf{z}_{i}^{\prime}\boldsymbol{\upbeta}+v_{i}+\varepsilon_{i},\qquad
i=1,\ldots,M,$ (5)

where $\boldsymbol{\upbeta}=(\beta_{1},\ldots,\beta_{P})^{\prime}$ is the
vector of fixed effects, the sampling errors $\varepsilon_{i}$ are assumed
independent with $\operatorname{{E}_{p}}(\varepsilon_{i})=0$ and
$\operatorname{{var}_{p}}(\varepsilon_{i})=\psi_{i}$, and random domain
effects $v_{i}$ are assumed independent of these errors. The latter effects
are supposed to be independent and identically distributed with
$\operatorname{E}(v_{i})=0$ and
$\operatorname{var}(v_{i})=\sigma_{v}^{2}\geqslant 0$ in respect of a
distribution, different from that generated by the design $p(\cdot)$.

Treating the estimates $\hat{\psi}_{i}$ as given numbers, the method of EBLUP
leads to the predictions of domain proportions (1) that are expressed as the
linear combinations [7]

$\hat{\theta}_{i}^{\mathrm{FH}}=\hat{\theta}_{i}^{\mathrm{FH}}(\hat{\psi}_{i})=\hat{\gamma}_{i}\hat{\theta}_{i}^{\mathrm{d}}+(1-\hat{\gamma}_{i})\mathbf{z}_{i}^{\prime}\boldsymbol{\hat{\upbeta}},\quad\text{with}\quad\hat{\gamma}_{i}=\frac{\hat{\sigma}_{v}^{2}}{\hat{\psi}_{i}+\hat{\sigma}_{v}^{2}},\qquad
i=1,\ldots,M,$ (6)

and

$\boldsymbol{\hat{\upbeta}}=\Biggl{(}\sum_{i=1}^{M}\frac{\mathbf{z}_{i}\mathbf{z}_{i}^{\prime}}{\hat{\psi}_{i}+\hat{\sigma}_{v}^{2}}\Biggr{)}^{-1}\sum_{i=1}^{M}\frac{\mathbf{z}_{i}\hat{\theta}_{i}^{\mathrm{d}}}{\hat{\psi}_{i}+\hat{\sigma}_{v}^{2}},$

where $\hat{\sigma}_{v}^{2}$ is an estimator of the variance $\sigma_{v}^{2}$.
One of the ways to estimate $\sigma_{v}^{2}$ is the estimator
$\hat{\sigma}_{v}^{2}$ based on the method of moments, as originally proposed
by [7]. For this estimator, approximately unbiased estimators of MSEs of (6)
were derived in [4]:

$\begin{split}\operatorname{mse}(\hat{\theta}_{i}^{\mathrm{FH}})=\hat{\gamma}_{i}\hat{\psi}_{i}+(1-\hat{\gamma}_{i})^{2}&\Biggl{[}\mathbf{z}_{i}^{\prime}\Biggl{(}\sum_{j=1}^{M}\frac{\mathbf{z}_{j}\mathbf{z}_{j}^{\prime}}{\hat{\psi}_{j}+\hat{\sigma}_{v}^{2}}\Biggr{)}^{-1}\mathbf{z}_{i}+\frac{4M}{\hat{\psi}_{i}+\hat{\sigma}_{v}^{2}}\Biggl{(}\sum_{j=1}^{M}\frac{1}{\hat{\psi}_{j}+\hat{\sigma}_{v}^{2}}\Biggr{)}^{-2}\Biggr{.}\\\
&\Biggl{.}-2\hat{\sigma}_{v}^{2}\Biggl{(}\sum_{j=1}^{M}\hat{\gamma}_{j}\Biggr{)}^{-3}\Biggl{\\{}M\sum_{j=1}^{M}\hat{\gamma}_{j}^{2}-\Biggl{(}\sum_{j=1}^{M}\hat{\gamma}_{j}\Biggr{)}^{2}\Biggr{\\}}\Biggr{]},\qquad
i=1,\ldots,M.\end{split}$ (7)

Predictions (6) and their MSE estimators (7) depend also on the estimators
$\hat{\psi}_{i}$ of the sampling variances $\psi_{i}$ of
$\hat{\theta}_{i}^{\mathrm{d}}$. However, direct estimators
$\hat{\psi}_{i}^{\mathrm{d}}$ of $\psi_{i}$, as, for example, approximately
design unbiased estimators (4) of (3) for $\hat{\theta}_{i}^{\mathrm{H}}$,
have large variances themselves for small sample sizes. Therefore, the direct
estimates $\hat{\psi}_{i}^{\mathrm{d}}$ are smoothed and new more stable
estimates $\hat{\psi}_{i}^{\mathrm{s}}$ are used in (6) and (7). It is called
the generalized variance function (GVF) approach [19]. The specific example of
the GVF method, similar to that used for estimation of census undercounts in
[5], is to assume that $\psi_{i}\approx KN_{i}^{\gamma}$ and estimate the
parameters $K>0$ and $\gamma\in\mathbb{R}$ using the regression model

$\log(\hat{\psi}_{i}^{\mathrm{d}})=\log(K)+\gamma\log(N_{i})+\eta_{i},\qquad
i=1,\ldots,M,$

where errors $\eta_{i}$ are independent and identically distributed. That is,
the smoothed estimates

$\hat{\psi}_{i}^{\mathrm{sD}}=\widehat{K}N_{i}^{\hat{\gamma}},\qquad
i=1,\ldots,M,$ (8)

of $\psi_{i}$ are based on the ordinary least squares estimates of the
regression parameters. Other smoothing examples are pooled variance estimation
[1] and a nonparametric smoothing like in [8]. Despite the smoothing,
estimators (7) tend to underestimate MSEs of (6) because the estimation of the
sampling variances $\psi_{i}$ is ignored in the derivation of (7).

##  4 Design-based composite estimation

###  4.1 Evaluation of optimal compositions and their accuracy estimation

Let us exclude the random effects $v_{i}$ from FH model (5). Then the model
becomes

$\hat{\theta}_{i}^{\mathrm{d}}=\mathbf{z}_{i}^{\prime}\boldsymbol{\upbeta}+\varepsilon_{i},\qquad
i=1,\ldots,M,$ (9)

and, using the estimates $\hat{\psi}_{i}$ of the variances $\psi_{i}$, we
arrive to the regression-synthetic estimators

$\hat{\theta}_{i}^{\mathrm{S}}=\hat{\theta}_{i}^{\mathrm{S}}(\hat{\psi}_{i})=\mathbf{z}_{i}^{\prime}\boldsymbol{\hat{\upbeta}},\qquad
i=1,\ldots,M,$ (10)

of the domain proportions $\theta_{i}$, where

$\boldsymbol{\hat{\upbeta}}=\Biggl{(}\sum_{i=1}^{M}\frac{\mathbf{z}_{i}\mathbf{z}_{i}^{\prime}}{\hat{\psi}_{i}}\Biggr{)}^{-1}\sum_{i=1}^{M}\frac{\mathbf{z}_{i}\hat{\theta}_{i}^{\mathrm{d}}}{\hat{\psi}_{i}}$
(11)

is the generalized least squares estimate of $\boldsymbol{\upbeta}$. Here, as
for EBLUPs, the use of smoothed estimates
$\hat{\psi}_{i}=\hat{\psi}_{i}^{\mathrm{s}}$ instead of
$\hat{\psi}_{i}^{\mathrm{d}}$ stabilizes synthetic estimators (10).

Estimators (10) rely on a synthetic assumption that the parameter
$\boldsymbol{\upbeta}$ is the same across all domains. Therefore, having a
good regression model, their sampling variances are small, compared to that of
the direct estimators $\hat{\theta}_{i}^{\mathrm{d}}$ or even the EBLUPs
$\hat{\theta}_{i}^{\mathrm{FH}}$. However, the design biases of (10) can be
relatively large if the synthetic assumption is not realistic. To find a
trade-off between larger variances of $\hat{\theta}_{i}^{\mathrm{d}}$ and
biases of the synthetic estimators $\hat{\theta}_{i}^{\mathrm{S}}$, we
consider their linear combinations

$\tilde{\theta}_{i}^{\mathrm{C}}=\tilde{\theta}_{i}^{\mathrm{C}}(\lambda_{i})=\lambda_{i}\hat{\theta}_{i}^{\mathrm{d}}+(1-\lambda_{i})\hat{\theta}_{i}^{\mathrm{S}},\qquad
i=1,\ldots,M,$ (12)

with weights $0\leqslant\lambda_{i}\leqslant 1$. Minimizing the function
$\operatorname{{MSE}_{p}}(\tilde{\theta}_{i}^{\mathrm{C}}(\lambda_{i}))$ with
respect to $\lambda_{i}$, the optimal weight for the domain ${\cal U}_{i}$ is
the population parameter [14, Section 3.3]

$\lambda_{i}^{*}=\frac{\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})-C_{i}}{\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{d}})+\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})-2C_{i}}\quad\text{with}\quad
C_{i}=\operatorname{{E}_{p}}(\hat{\theta}_{i}^{\mathrm{d}}-\theta_{i})(\hat{\theta}_{i}^{\mathrm{S}}-\theta_{i}).$
(13)

Assuming that
$|C_{i}|\ll\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})$, the
approximation
$\lambda_{i}^{*}\approx\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})/(\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{d}})+\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}}))$
is applied, but the further difficulty is to evaluate the quantities
$\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})$. A common approach
to this is to use the representation [14, Section 3.2.5]

$\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{S}})=\operatorname{{E}_{p}}(\hat{\theta}_{i}^{\mathrm{S}}-\hat{\theta}_{i}^{\mathrm{d}})^{2}-\operatorname{{var}_{p}}(\hat{\theta}_{i}^{\mathrm{S}}-\hat{\theta}_{i}^{\mathrm{d}})+\operatorname{{var}_{p}}(\hat{\theta}_{i}^{\mathrm{S}}),$
(14)

where $\hat{\theta}_{i}^{\mathrm{d}}$ is assumed to be unbiased, and then to
build an approximately design unbiased estimator

$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{S}})=(\hat{\theta}_{i}^{\mathrm{S}}-\hat{\theta}_{i}^{\mathrm{d}})^{2}-\hat{\sigma}^{2}(\hat{\theta}_{i}^{\mathrm{S}}-\hat{\theta}_{i}^{\mathrm{d}})+\hat{\sigma}^{2}(\hat{\theta}_{i}^{\mathrm{S}})$
(15)

of (14), where $\hat{\sigma}^{2}(\cdot)$ is an estimator of the design
variance $\operatorname{{var}_{p}}(\cdot)$. Unfortunately, estimator (15) can
be very unstable and take negative values for individual small domains.
Therefore, the straightforward estimation of optimal weights (13) is avoided.

To evaluate the optimal coefficients for compositions (12), one can set a
common weight for all domains and then minimize a total MSE with respect to
that weight [13]. A similar approach is to apply James–Stein method [14,
Section 3.4]. One more idea is SSD estimation [6], where estimators of the
weights in (12) are taken to be of the form

$\hat{\lambda}_{i}=\hat{\lambda}_{i}(\delta)=\begin{cases}1&\text{if
$\widehat{N}_{i}/N_{i}\geqslant\delta$,}\\\ \widehat{N}_{i}/(\delta
N_{i})&\text{otherwise.}\end{cases}$ (16)

These weights are dependent on the single subjectively chosen parameter
$\delta$ for all domains with default value $\delta=1$. Similar SSD estimators
were derived in [16] applying a conditional analysis.

Estimation of MSEs of the design-based composite estimators like these is
known as a difficult problem in the literature [14, Chapter 3]. One general
solution is to treat the composition
$\hat{\theta}_{i}^{\mathrm{C}}=\tilde{\theta}_{i}^{\mathrm{C}}(\hat{\lambda}_{i})$
as a synthetic estimator and use the estimator

$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{C}})=(\hat{\theta}_{i}^{\mathrm{C}}-\hat{\theta}_{i}^{\mathrm{d}})^{2}-\hat{\sigma}^{2}(\hat{\theta}_{i}^{\mathrm{C}}-\hat{\theta}_{i}^{\mathrm{d}})+\hat{\sigma}^{2}(\hat{\theta}_{i}^{\mathrm{C}})$
(17)

of $\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{C}})$, see [14, Example
3.3.1] and [2]. However, this estimator has the same drawbacks as (15).
Another general method is to assume that the estimator
$\hat{\theta}_{i}^{\mathrm{C}}$ defined by (12) approximates the optimal
combination
$\hat{\theta}_{i}^{\mathrm{opt}}=\tilde{\theta}_{i}^{\mathrm{C}}(\lambda_{i}^{*})$
quite well and derive the approximation [3]

$\operatorname{{MSE}_{p}}(\hat{\theta}_{i}^{\mathrm{C}})\approx\lambda_{i}(1-\lambda_{i})\psi_{i}+\operatorname{{var}_{p}}(\hat{\theta}_{i}^{\mathrm{C}})$

with the empirical version

$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{C}})=\hat{\lambda}_{i}(1-\hat{\lambda}_{i})\hat{\psi}_{i}+\hat{\sigma}^{2}(\hat{\theta}_{i}^{\mathrm{C}}),$
(18)

where we would set $\hat{\psi}_{i}=\hat{\psi}_{i}^{\mathrm{s}}$ to have robust
MSE estimators. Estimator (18) takes only non-negative values.

###  4.2 Composition based on a ratio of variances

The sampling variance $\psi_{i}$ is approximately proportional to the product
$\theta_{i}(1-\theta_{i})$. That is, one can use the approximation

$\psi_{i}\approx\frac{D_{i}\theta_{i}(1-\theta_{i})}{n_{i}},$ (19)

where $D_{i}$ is the design effect reflecting the sample efficiency of the
complex sampling design [9]. Then, inserting $\hat{\theta}_{i}^{\mathrm{d}}$
and an appropriate estimator $\widehat{D}_{i}$ of $D_{i}$ into (19), we would
approximate the direct estimator $\hat{\psi}_{i}^{\mathrm{d}}$ of $\psi_{i}$.
Let us first suppose that the domain proportions $\theta_{i}$ are small, say
$\theta_{i}<0.1$. In that case, it is even more complicated to get reliable
direct estimates and estimates of their accuracy [10]. For example, the direct
estimator of the proportion of the unemployed can take zero value even for a
sample of moderate size in the municipality.

Consider two candidate estimators $\hat{\psi}_{i}^{\mathrm{d}}$ and
$\hat{\psi}_{i}^{\mathrm{s}}$ of $\psi_{i}$ used in regression-synthetic
estimator (10). Assume that we got too small estimate
$\hat{\theta}_{i}^{\mathrm{d}}$ of $\theta_{i}$ for the specific sample $s$.
The direct estimate $\hat{\psi}_{i}^{\mathrm{d}}$ then underestimates the
sampling variance $\psi_{i}$. Therefore, the inequality
$\hat{\psi}_{i}^{\mathrm{s}}>\hat{\psi}_{i}^{\mathrm{d}}$ should often hold,
that is, the smoothed variance $\hat{\psi}_{i}^{\mathrm{s}}$ could be a better
choice than $\hat{\psi}_{i}^{\mathrm{d}}$. Now suppose that
$\hat{\theta}_{i}^{\mathrm{d}}$ overestimated the parameter $\theta_{i}$. Then
$\hat{\psi}_{i}^{\mathrm{d}}$ overestimates $\psi_{i}$ as well, and the
inequality $\hat{\psi}_{i}^{\mathrm{s}}<\hat{\psi}_{i}^{\mathrm{d}}$ should
hold if $\hat{\theta}_{i}^{\mathrm{d}}$ is an outlier. That larger estimate
$\hat{\psi}_{i}^{\mathrm{d}}$ can be employed to down-weight the outlying
observation $\hat{\theta}_{i}^{\mathrm{d}}$ used in (11) thus robustifying
synthetic estimators (10). From these considerations, we derive the combined
estimators

$\hat{\psi}_{i}^{\mathrm{c}}=\max\\{\hat{\psi}_{i}^{\mathrm{s}},\hat{\psi}_{i}^{\mathrm{d}}\\},\qquad
i=1,\ldots,M,$

of the sampling variances $\psi_{i}$ that should improve the regression-
synthetic estimation. Next, in line with the same ideas, we define the design-
based composite estimators

$\hat{\theta}_{i}^{\mathrm{C}}=\hat{\lambda}_{i}\hat{\theta}_{i}^{\mathrm{d}}+(1-\hat{\lambda}_{i})\hat{\theta}_{i}^{\mathrm{S}}(\hat{\psi}_{i}^{\mathrm{c}})\quad\text{with}\quad\hat{\lambda}_{i}=\frac{\min\\{\hat{\psi}_{i}^{\mathrm{s}},\hat{\psi}_{i}^{\mathrm{d}}\\}}{\hat{\psi}_{i}^{\mathrm{c}}},\qquad
i=1,\ldots,M,$ (20)

of domain proportions (1). If the estimate $\hat{\theta}_{i}^{\mathrm{d}}$ is
an outlier by its small or large value, then relatively more weight is
attached to the synthetic part of composition (20). The composition is a
shrinkage estimator because it shrinks the direct estimator toward the
synthetic one.

We apply the same arguments to create (20) if the parameters $\theta_{i}$ are
not small, but then the inequalities
$\max\\{\theta_{i},\hat{\theta}_{i}^{\mathrm{d}}\\}<1/2$ or
$\min\\{\theta_{i},\hat{\theta}_{i}^{\mathrm{d}}\\}>1/2$ must be satisfied. If
these inequalities are not valid, the composite estimator is still applicable,
but it can be less efficient. The worst scenario here would be a large
difference $\theta_{i}-\hat{\theta}_{i}^{\mathrm{d}}$ and the relation
$\theta_{i}\approx 1-\hat{\theta}_{i}^{\mathrm{d}}$ but those events are rare.

To estimate MSE of composition (20), we suggest to apply general estimator
(18). We study the accuracy of both these estimators in Section 5.

###  4.3 Sample-size-dependent estimation

A choice of the parameter $\delta$ in (16) varies from survey to survey. That
is, the values 2/3 and 1 are good for LFS in [6], the authors of [18] try the
larger points 1.5 and 2 for their data, and optimal values of $\delta$ are
even higher in [2]. Therefore, to select the value of the parameter for the
composition
$\tilde{\theta}_{i}^{\mathrm{C}}(\delta)=\tilde{\theta}_{i}^{\mathrm{C}}(\hat{\lambda}_{i}(\delta))$
defined by (12), we minimize numerically the sample based function [2]

$r(\delta)=\frac{1}{M}\sum_{i=1}^{M}\operatorname{mse_{u}}(\tilde{\theta}_{i}^{\mathrm{C}}(\delta))$
(21)

with respect to $\delta$. This function is the average of individual MSE
estimators (17) over domains and therefore it is stable unlike the individual
ones. So we get the adaptive composite estimators

$\hat{\theta}_{i}^{\mathrm{SSD}}=\tilde{\theta}_{i}^{\mathrm{C}}(\hat{\delta}^{*})\quad\text{with}\quad\hat{\delta}^{*}=\operatorname*{arg\,min}_{\delta>0}r(\delta),\qquad
i=1,\ldots,M,$ (22)

of the domain proportions. We apply estimators (18) to evaluate MSEs of these
compositions.

##  5 Simulations using the Labor Force Survey data

The main LFS variable is the categorical one that indicates an individual’s
participation in the labor market. This variable is decomposed into three
binary variables: is the person unemployed, is employed, and is not in the
labor force. We are going to estimate the proportions of the former two
variables in the municipalities of Lithuania. To imitate the real survey, we
construct the artificial population from the sample data of the fourth quarter
of 2018 as follows: we remove municipalities with too small fractions of
observed unemployed persons and then replicate the data of each individual the
number of times equal to the rounded survey weight. The size of that
population ${\cal U}$ is $N=1396763$, and it contains $M=30$ municipalities.
In LFS, the sample of households is drawn without replacement with
probabilities proportional to the number of their members, and then the
selected households are surveyed entirely. We use the same sampling design to
draw $R=10^{3}$ independent samples of households of size $n^{\prime}=3700$,
and it yields the samples of persons of sizes close to $n=7667$. Then, for the
$k$th individual that belongs to the $l$th household of size $h_{l}$, we apply
the approximation $\pi_{k}\approx h_{l}n^{\prime}/N$, $k\in{\cal U}$.

We compare the direct estimator
$\hat{\theta}_{i}^{\mathrm{d}}=\hat{\theta}_{i}^{\mathrm{H}}$ from (2),
regression-synthetic estimator (10), EBLUP (6) calculated using the package
sae for R [11], and two design-based compositions (20) and (22). Moreover, we
compare the accuracy of MSE estimator (7) for (6) with that of two MSE
estimators (17) and (18) applied to both compositions (20) and (22). We
consider also the optimal combination $\hat{\theta}_{i}^{\mathrm{opt}}$ that
uses (13) in (12), and its MSE estimator calculated by (18).

To model the direct estimates of the proportions of interest by (5) and (9),
we use the municipality characteristics
$\mathbf{z}_{i}=(1,z_{i2},z_{i3},z_{i4},z_{i5},z_{i6})^{\prime}$, where
$z_{i2}$ is the proportion of registered unemployed individuals derived from
the administrative Lithuanian Labor Exchange data, $z_{i3}$ is the proportion
of persons who, according to the register of the State Social Insurance Fund
Board, paid the social contribution one month before they participated in the
survey, $z_{i4}$ is the proportion of males, and $z_{i5}$ and $z_{i6}$ are the
proportions of individuals from age intervals 26–40 and 41–55, respectively.

Since the sampling fractions are small in the municipalities, we take
$\pi_{kl}\approx\pi_{k}\pi_{l}$, $k\neq l$, and so approximate direct
estimators (4) of sampling variances (3) by

$\hat{\psi}_{i}^{\mathrm{H}}\approx\hat{\psi}_{i}^{\mathrm{d}}=\frac{1}{\widehat{N}_{i}^{2}}\sum_{k\in
s_{i}}w_{k}(w_{k}-1)(y_{k}-\hat{\theta}_{i}^{\mathrm{d}})^{2},\qquad
i=1,\ldots,M,$

where we write $w_{k}=1/\pi_{k}$. Then we smooth these
$\hat{\psi}_{i}^{\mathrm{d}}$ to obtain
$\hat{\psi}_{i}=\hat{\psi}_{i}^{\mathrm{sD}}$ according to (8) and use the
smoothed estimates for (6), (7), (10), (18), and in the synthetic parts of
(22) and $\hat{\theta}_{i}^{\mathrm{opt}}$.

We apply the bootstrap method of [15] to evaluate the estimators of the design
variances in (17), (18), and (21). Denote by $\hat{\theta}_{i}$ any estimator
for which we need to estimate the design variance. The bootstrap procedure
works as follows: (i) Draw a simple random sample of $m=n^{\prime}-1$
households with replacement from $n^{\prime}$ sample households. Let
$m_{l}^{*}$ be the number of times the $l$th sample household is selected, and
then $\sum_{l=1}^{n^{\prime}}m_{l}^{*}=m$. Define the bootstrap weights
$w_{l}^{*}=n^{\prime}m_{l}^{*}w_{l}/m$, $l=1,\ldots,n^{\prime}$. Calculate the
bootstrap estimate $\hat{\theta}_{i}^{*}$ using the weights $w_{l}^{*}$ in the
formula for $\hat{\theta}_{i}$. (ii) Repeat step (i) $B$ times independently
to obtain the estimates $\hat{\theta}_{i}^{*(b)}$, $b=1,\ldots,B$. Then

$\hat{\sigma}^{2}(\hat{\theta}_{i})=\frac{1}{B}\sum_{b=1}^{B}(\hat{\theta}_{i}^{*(b)}-\bar{\theta}_{i}^{*})^{2},\quad\text{where}\quad\bar{\theta}_{i}^{*}=\frac{1}{B}\sum_{b=1}^{B}\hat{\theta}_{i}^{*(b)},$

is the bootstrap estimator of the design variance
$\operatorname{{var}_{p}}(\hat{\theta}_{i})$. We take $B=200$.

We evaluate all estimators for each of $R$ samples and calculate
approximations to their root mean squared errors (RMSEs) and absolute biases
(ABs). That is, we use the accuracy measures

$\operatorname{RMSE}(\hat{\mu}_{i})=\Biggl{(}\frac{1}{R}\sum_{r=1}^{R}(\hat{\mu}_{i}^{(r)}-\mu_{i})^{2}\Biggr{)}^{1/2}\quad\text{and}\quad\operatorname{AB}(\hat{\mu}_{i})=\Biggl{|}\frac{1}{R}\sum_{r=1}^{R}\hat{\mu}_{i}^{(r)}-\mu_{i}\Biggr{|},$
(23)

where $\hat{\mu}_{i}^{(r)}$ is a realization of the specific estimator
$\hat{\mu}_{i}$ of the parameter $\mu_{i}$, based on the $r$th sample. We
classify the municipalities by the expected domain sample size into three
classes of equal size, and calculate the average of RMSEs as well as ABs over
domains of each class. We also present the averages of (23) over all
municipalities as common indicators of accuracy.

The results for the proportions of the unemployed and the employed are
presented in Tables 1 and 2, respectively. Let us use the superscripts of
estimators to discuss the output. In both the tables, any indirect estimator
of the proportions improves the direct one in the sense of RMSE, and
theoretical composition opt is the best estimator. Among the indirect
estimators, regression-synthetic estimator S has much larger design biases
than compositions FH, C, and SSD. In Table 1, the averages of RMSEs over all
domains of design-based composite estimators C and SSD are smaller than that
of EBLUP FH. It is not valid for estimator C in Table 2 because the
proportions of the employed are distributed near the point $1/2$ if to look at
the five-number summary $(0.379,0.585,0.634,0.668,0.766)$ for the true
proportions.

MSE estimation (18) for design-based compositions C and SSD evidently improves
estimation (17), and yields similar or even better results than MSE estimator
(7) for FH. The best MSE estimation using (18) is obtained for optimal
composition opt. Composite estimators C and SSD only approximate the optimal
one and, therefore, their MSE estimators have larger errors. On the other
hand, these errors are acceptable if to compare them with the results for FH.

The same experiment but with the twice smaller sample size $n^{\prime}=1850$
leads to similar conclusions. In this case, design-based composition C
improves EBLUP FH more for small proportions.

##  6 Conclusions

The construction of composite estimator (20) is based on the monotonicity of
the variance of the direct estimator as the function of the proportion.
Approximation (19) is the monotone function in two separate parts of the
interval $[0,1]$. Therefore, the composition loses its efficiency for the
proportions close to turning point $1/2$, where the monotonicity changes.

In general, the sampling variance of any direct estimator of the domain mean
is not the monotone function of the target parameter. On the other hand, some
GVF models from [19, p. 274] suggest that this function might be treated as an
approximately monotonic one. Therefore, if we can find the GVF model that fits
the data well and is the monotonic function, then estimator (20) could be also
applied to the domain means with this fitted model used instead of smoothed
sampling variances (8).

The simulation study shows that the design-based compositions might be an
alternative to the classical EBLUP estimating proportions in small domains.
Adaptive composite estimator SSD works well for both unemployment and
employment cases, while simpler composition (20) is efficient for the
unemployment fractions that are small proportions.

Design-based estimators and estimators of MSE under the design-based approach
are desirable in practice [12]. That design MSE estimator (18) works well in
our simulations, and its formula is simple compared to that of model MSE
estimator (7) for EBLUP.

## References

  * [1] H.J. Boonstra, J.A. van den Brakel, B. Buelens, S. Krieg, M. Smeets,  Towards small area estimation at Statistics Netherlands,  Metron, 66(1):21–49, 2008.
  * [2] A. Čiginas,  Adaptive composite estimation in small domains,  Nonlinear Anal. Model. Control, 25(3):341–357, 2020.
  * [3] A. Čiginas,  Design-based composite estimation rediscovered,  arXiv:2108.05052 [stat.ME], 2021.
  * [4] G.S. Datta, J.N.K. Rao, D.D. Smith,  On measuring the variability of small area estimators under a basic area level model,  Biometrika, 92(1):183–196, 2005.
  * [5] P. Dick,  Modelling net undercoverage in the 1991 Canadian census,  Surv. Methodol., 21(1):45–54, 1995.
  * [6] J.D. Drew, M.P. Singh, G.H. Choudhry,  Evaluation of small area estimation techniques for the Canadian Labour Force Survey,  Surv. Methodol., 8:17–47, 1982.
  * [7] R.E. Fay, R.A. Herriot,  Estimates of income for small places: an application of James-Stein procedures to census data,  J. Amer. Statist. Assoc., 74(366):269–277, 1979.
  * [8] W. González-Manteiga, M.J. Lombardía, I. Molina, D. Morales, L. Santamaría,  Small area estimation under Fay–Herriot models with non-parametric estimation of heteroscedasticity,  Stat. Model., 10(2):215–239, 2010.
  * [9] L. Kish,  Methods for design effects,  J. Off. Stat., 11(1):55–77, 1995.
  * [10] E.L. Korn, B.I. Graubard,  Confidence intervals for proportions with small expected number of positive counts estimated from survey data,  Surv. Methodol., 24(2):193–201, 1998.
  * [11] I. Molina, Y. Marhuenda,  sae: An R package for small area estimation,  R J., 7(1):81–98, 2015,  available from: https://journal.r-project.org/archive/2015/RJ-2015-007/RJ-2015-007.pdf.
  * [12] I. Molina, E. Strzalkowska-Kominiak,  Estimation of proportions in small areas: application to the labour force using the Swiss Census Structural Survey,  J. Roy. Statist. Soc. Ser. A, 183(1):281–310, 2020.
  * [13] N.J. Purcell, L. Kish,  Estimation for small domains,  Biometrics, 35:365–384, 1979.
  * [14] J.N.K. Rao, I. Molina,  Small Area Estimation,  John Wiley, New Jersey, 2 edition, 2015.
  * [15] J.N.K. Rao, C.F.J. Wu, K. Yue,  Some recent work on resampling methods for complex surveys,  Surv. Methodol., 18(2):209–217, 1992.
  * [16] C.-E. Särndal, M.A. Hidiroglou,  Small domain estimation: a conditional analysis,  J. Amer. Statist. Assoc., 84(405):266–275, 1989.
  * [17] C.-E. Särndal, B. Swensson, J. Wretman,  Model Assisted Survey Sampling,  Springer-Verlag, New York, 1992.
  * [18] M.D. Ugarte, T. Goicoa, A.F. Militino, M. Sagaseta-López,  Estimating unemployment in very small areas,  SORT, 33(1):49–70, 2009.
  * [19] K.M. Wolter,  Introduction to Variance Estimation,  Springer-Verlag, New York, 2 edition, 2007.

Table 1: Average RMSEs and ABs of estimators for the unemployed proportions in domain size classes as $n\approx 7667$. The domain is small if its expected sample size $\bar{n}_{i}=\operatorname{{E}_{p}}(n_{i})<116$, is medium for $116\leqslant\bar{n}_{i}<159$, and is large as $\bar{n}_{i}\geqslant 159$. | Average RMSE ($\times 10^{2}$) | Average AB ($\times 10^{2}$)
---|---|---
Estimator | Domain size class by $\bar{n}_{i}$ | Domain size class by $\bar{n}_{i}$
|  any |  small |  medium |  large |  any |  small |  medium |  large
$\hat{\theta}_{i}^{\mathrm{d}}$ |  2.4793 |  3.8540 |  2.4578 |  1.1259 |  0.0636 |  0.1200 |  0.0485 |  0.0223
$\hat{\theta}_{i}^{\mathrm{S}}$ |  1.8174 |  2.8950 |  1.5632 |  0.9940 |  1.3461 |  2.3656 |  1.0677 |  0.6050
$\hat{\theta}_{i}^{\mathrm{FH}}$ |  1.7857 |  2.6707 |  1.7156 |  0.9707 |  0.7349 |  1.4738 |  0.5496 |  0.1811
$\hat{\theta}_{i}^{\mathrm{C}}$ |  1.7511 |  2.6798 |  1.6838 |  0.8897 |  0.7951 |  1.4777 |  0.6130 |  0.2946
$\hat{\theta}_{i}^{\mathrm{SSD}}$ |  1.7529 |  2.7228 |  1.6162 |  0.9196 |  0.8649 |  1.4974 |  0.6928 |  0.4045
$\hat{\theta}_{i}^{\mathrm{opt}}$ |  1.4712 |  2.3804 |  1.2486 |  0.7846 |  0.7301 |  1.3978 |  0.5206 |  0.2720
$\operatorname{mse}(\hat{\theta}_{i}^{\mathrm{FH}})$ |  0.0223 |  0.0445 |  0.0173 |  0.0051 |  0.0180 |  0.0373 |  0.0128 |  0.0039
$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{C}})$ |  0.0708 |  0.1540 |  0.0491 |  0.0094 |  0.0263 |  0.0532 |  0.0215 |  0.0041
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{C}})$ |  0.0173 |  0.0371 |  0.0119 |  0.0030 |  0.0135 |  0.0296 |  0.0087 |  0.0021
$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{SSD}})$ |  0.0593 |  0.1164 |  0.0494 |  0.0120 |  0.0115 |  0.0290 |  0.0051 |  0.0005
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{SSD}})$ |  0.0257 |  0.0497 |  0.0210 |  0.0063 |  0.0172 |  0.0314 |  0.0153 |  0.0050
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{opt}})$ |  0.0098 |  0.0206 |  0.0064 |  0.0023 |  0.0050 |  0.0110 |  0.0027 |  0.0012
Table 2: Average RMSEs and ABs of estimators for the employed proportions in domain size classes as $n\approx 7667$. The domain is small if its expected sample size $\bar{n}_{i}=\operatorname{{E}_{p}}(n_{i})<116$, is medium for $116\leqslant\bar{n}_{i}<159$, and is large as $\bar{n}_{i}\geqslant 159$. | Average RMSE ($\times 10^{2}$) | Average AB ($\times 10^{2}$)
---|---|---
Estimator | Domain size class by $\bar{n}_{i}$ | Domain size class by $\bar{n}_{i}$
|  any |  small |  medium |  large |  any |  small |  medium |  large
$\hat{\theta}_{i}^{\mathrm{d}}$ |  4.7718 |  6.9201 |  4.7104 |  2.6848 |  0.1516 |  0.2577 |  0.1395 |  0.0575
$\hat{\theta}_{i}^{\mathrm{S}}$ |  3.4061 |  4.9905 |  3.2215 |  2.0061 |  2.6481 |  4.1247 |  2.5006 |  1.3188
$\hat{\theta}_{i}^{\mathrm{FH}}$ |  3.3054 |  4.6679 |  3.1768 |  2.0716 |  1.7276 |  2.8992 |  1.6535 |  0.6302
$\hat{\theta}_{i}^{\mathrm{C}}$ |  4.2265 |  6.0532 |  4.1539 |  2.4724 |  0.4024 |  0.6893 |  0.4213 |  0.0967
$\hat{\theta}_{i}^{\mathrm{SSD}}$ |  3.2747 |  4.7502 |  3.1425 |  1.9314 |  1.6996 |  2.6130 |  1.6237 |  0.8622
$\hat{\theta}_{i}^{\mathrm{opt}}$ |  2.8602 |  4.1800 |  2.6261 |  1.7746 |  1.6026 |  2.5092 |  1.4717 |  0.8269
$\operatorname{mse}(\hat{\theta}_{i}^{\mathrm{FH}})$ |  0.0724 |  0.1281 |  0.0670 |  0.0221 |  0.0561 |  0.1070 |  0.0469 |  0.0143
$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{C}})$ |  0.0666 |  0.1371 |  0.0481 |  0.0144 |  0.0297 |  0.0590 |  0.0216 |  0.0086
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{C}})$ |  0.0535 |  0.1093 |  0.0406 |  0.0107 |  0.0102 |  0.0170 |  0.0120 |  0.0016
$\operatorname{mse_{u}}(\hat{\theta}_{i}^{\mathrm{SSD}})$ |  0.1992 |  0.3628 |  0.1784 |  0.0563 |  0.0297 |  0.0674 |  0.0201 |  0.0014
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{SSD}})$ |  0.0655 |  0.1091 |  0.0697 |  0.0178 |  0.0442 |  0.0715 |  0.0491 |  0.0119
$\operatorname{mse_{b}}(\hat{\theta}_{i}^{\mathrm{opt}})$ |  0.0181 |  0.0374 |  0.0110 |  0.0059 |  0.0082 |  0.0155 |  0.0049 |  0.0043