# Enhancing direct exoplanet spectroscopy with apodizing and beam shaping
optics

Benjamin Calvin California Institute of Technology, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  Nemanja Jovanovic California Institute of Technology,
1200 E. California Blvd., Pasadena, CA 91125, USA  Garreth Ruane California
Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA  Jacklyn Pezzato  California Institute of
Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA  Jennah Colborn
California Institute of Technology, 1200 E. California Blvd., Pasadena, CA
91125, USA  Daniel Echeverri  California Institute of Technology, 1200 E.
California Blvd., Pasadena, CA 91125, USA  Tobias Schofield  California
Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
Michael Porter  California Institute of Technology, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  J. Kent Wallace  Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109,
USA  Jacques-Robert Delorme California Institute of Technology, 1200 E.
California Blvd., Pasadena, CA 91125, USA  W. M. Keck Observatory, Kamuela, HI
96743, USA  Dimitri Mawet  California Institute of Technology, 1200 E.
California Blvd., Pasadena, CA 91125, USA  Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109,
USA

###### Abstract

Direct exoplanet spectroscopy aims to measure the spectrum of an exoplanet
while simultaneously minimizing the light collected from its host star.
Isolating the planet light from the starlight improves the signal-to-noise
ratio (S/N) per spectral channel when noise due to the star dominates, which
may enable new studies of the exoplanet atmosphere with unprecedented detail
at high spectral resolution ($>$30,000). However, the optimal instrument
design depends on the flux level from the planet and star compared to the
noise due to other sources, such as detector noise and thermal background.
Here we present the design, fabrication, and laboratory demonstration of
specially-designed optics to improve the S/N in two potential regimes in
direct exoplanet spectroscopy with adaptive optics instruments. The first is a
pair of beam-shaping lenses that increase the planet signal by improving the
coupling efficiency into a single-mode fiber at the known position of the
planet. The second is a grayscale apodizer that reduces the diffracted
starlight for planets at small angular separations from their host star. The
former especially increases S/N when dominated by detector noise or thermal
background, while the latter helps reduce stellar noise. We show good
agreement between the theoretical and experimental point spread functions in
each case and predict the exposure time reduction ($\sim 33\%$) that each set
of optics provides in simulated observations of 51 Eridani b using the Keck
Planet Imager and Characterizer instrument at W.M. Keck Observatory.

astronomical instrumentation: high angular resolution, techniques: high
angular resolution, exoplanets

††journal: PASP††facilities: KeckII

##  1 Introduction

With over 4000 exoplanets confirmed to date, detection has given way to the
era of characterization, critical to understanding the properties of these
systems. In this vein, high contrast imagers (HCI), which isolate the light
from the planet, offer numerous advantages over indirect techniques that rely
on the signal from the host star alone. Typically, HCI use advanced wavefront
control techniques combined with coronagraphs to extinguish the star light and
minimize contamination to the planet signal (Macintosh et al., 2014; Jovanovic
et al., 2015; Males et al., 2018; Beuzit et al., 2019). These systems often
exploit low to medium resolving power (R$\sim$10-1,000), integral field
spectrographs to characterize the targets (e.g. the CHARIS instrument (Groff
et al., 2016)). Besides providing spectral information about the target,
dispersing the light offers two advantages. Firstly, the distance between a
speckle and the optical axis increases linearly with wavelength, while a
planet remains a fixed distance away, which can be used to differentiate
between the two. Indeed, this technique is referred to as spectral
differential imaging (SDI)(Sparks & Ford, 2002). Secondly, the spectrum of a
planet has a different profile to that of a star, which enables the observer
to constrain the nature of the orbiting object. This technique has been
exploited to detect or characterize both planets (Barman et al., 2015) as well
as disks (Currie et al., 2017).

Recently, the field has shifted focus towards much greater resolving power
(R$>$30,000), to really take advantage of SDI in the regime where spectral
lines in the target begin to be resolved (Snellen et al., 2015). The
technique, dubbed High dispersion Coronagraphy (HDC), optimally combines high
contrast imaging techniques such as adaptive optics/wavefront control plus
coronagraphy to high resolution spectroscopy (Wang et al., 2017; Mawet et al.,
2017). One approach that has been explored is to use an optical fiber to route
the light of the known planet from the focal planet to the spectrograph. A
single-mode fiber (SMF) is the ideal transport vehicle owing to the fact it
has a field-of-view (FOV) that can be matched to the 1 $\lambda/D$ width of
the point spread function (PSF), enabling efficient coupling for the planet
light and suppression of unwanted star light. Its small core size only allows
light to be guided in a single mode providing spatial filtering, which further
suppresses starlight speckles from coupling in and more importantly delivers
an ultra stable output beam profile, highly desired by many instruments
(Schwab et al., 2012; Woillez et al., 2003; Crass et al., 2020). In addition,
its narrow FOV reduces the amount of sky background injected into the
spectrograph, which can be several orders of magnitude greater in the case of
a seeing-limited spectrograph.

HDC provides the ability to do species-by-species molecular characterization
(e.g. oxygen, water, carbon dioxide, methane), thermal (vertical) atmospheric
structure, planetary spin measurements (length of day), and potentially
Doppler imaging of atmospheric (clouds) and/or surface features (continents
versus oceans) (Wang et al., 2017). As such several projects have been
initiated to realize this new technique from the ground including: the Keck
Planet Imager and Characterizer (KPIC) (Jovanovic et al., 2019), which
combines Keck AO and NIRSPEC, the Rigorous Exoplanetary Atmosphere
Characterization with High dispersion coronography instrument (REACH)
(Jovanovic et al., 2017a), which combines SCExAO and IRD and High-Resolution
Imaging and Spectroscopy of Exoplanets (HiRISE), which combines SPHERE and
CRIRES+(Vigan et al., 2018). The phase I version of KPIC and the REACH
instrument are both transitioning from commissioning to early science at the
time of writing of this article and offer complimentary wavelength coverage
across the near-IR (NIR) on Maunakea (REACH operates from y-H and KPIC
operates in K and L bands), while HiRISE is still in the development stage.

Key to being able to disperse the faint planet signal across many pixels and
maintain signal-to-noise (S/N), as is the case in a high resolution
spectrograph, is maximizing the planet throughput. This can be done by
improving the adaptive optics (AO) correction to boost the Strehl ratio and
hence coupling efficiency, reshaping the PSF to better match the fiber mode
(Jovanovic et al., 2017b), or minimizing the number of optics in the system.
In addition, it is also important to minimize the stellar leakage into the
planet fiber to reduce the photon noise contribution from the star. Improved
AO correction can also help here, but specific focal plane wavefront control
techniques such as speckle nulling (Bottom et al., 2016; Sayson et al., 2019)
or the use of a pupil plane apodizing mask (Zhang et al., 2018), are better
suited to reducing speckle noise. Although these technologies achieve
different goals, they are used to reduce total exposure time, which is
critical when characterizing faint exoplanets with high spectral resolution on
large telescopes.

The technology or approach of choice depends on the parameters of the system
being observed, such as planet-to-star flux ratio, angular separation, and the
waveband of observation. In this paper, we present the development of two
technologies that can help boost the performance of HDC on KPIC and similar
instruments, namely: Phase Induced Amplitude Apodization (PIAA) optics for re-
shaping the beam and boosting coupling and a grayscale microdot apodizer (MDA)
for suppressing diffraction features at small angular separations. We aim to
deploy these as part of the phase II upgrade of the KPIC instrument in late
2021. In Section 2 we outline the two technologies and present the design and
simulated performance of such devices. In Section 3 we present the
experimental setup and procedures to evaluate the two technologies and in
Section 4 we summarize and compare the results to simulations. Section 5
compares the two technologies in the context of a hypothetical observing
scenario to highlight the benefits of each in their respective domains and
Section 6 rounds out the paper with some concluding thoughts.

##  2 Apodization Optics

This section provides an overview of the two technologies that we have
developed to enhance the S/N of direct spectroscopic, or HDC, observations of
exoplanets: PIAA (see Fig. 1) and MDA (see Fig. 2) optics. Each has been
employed previously in apodized coronagraphs in the context of high-contrast
imaging (Watson et al., 1991; Nisenson & Papaliolios, 2001; Kasdin et al.,
2003; Soummer et al., 2003), where the PIAA (Guyon, 2003) or MDA (Dorrer &
Zuegel, 2007; Martinez et al., 2009a) modifies the shape of the beam at, or
near, the re-imaged telescope pupil to reduce diffracted starlight at small
angular separations where exoplanets may be directly imaged. These pupil-plane
optics do not require any focal plane masks to work for our purposes with
KPIC. Our application differs from conventional high contrast imaging in that
we aim to maximize the S/N in the measured spectrum of an exoplanet using a
diffraction-limited spectrograph that is fed by a SMF. As such, the PIAA and
MDA optics that we developed for KPIC are optimized to maximize the coupling
efficiency for planet light, $\eta_{p}$, and minimize the fraction of the
starlight that is coupled into the SMF, $\eta_{s}$.

Figure 1: a) A ray trace through the PIAA optics. Here, the 50 mm represents
the inner spacing between the powered surfaces of the PIAA optics. b) An image
of the first lens, which refracts the beam to reshape it. c) An image of the
second lens in the pair, which re-collimates the light after beam shaping. d)
The radial sag profile for the first lens. e) The radial sag profile for the
second lens. Both sag profiles are azimuthally invariant. Figure 2: (Left)
Design for the MDA. It consists of a 500$\times$500 grid of 25 $\mu$m chrome
squares forming a 12.3 mm annulus of transmission. The white areas show the
region of high transmission, while the black parts are covered with chrome.
For scale, the outer black edge of the MDA is 15 mm. (Right) A white-light
photograph of a device manufactured for KPIC placed on textured background
showing the transmission function in the visible regime appears similar to the
design.

The coupling efficiency of a coherent, scalar field, $E(\mathbf{r})$, into a
fiber mode $\Psi(\mathbf{r})$ is given by:

$\eta=\frac{\left|\int
E^{*}(\mathbf{r})\Psi(\mathbf{r})dA\right|^{2}}{\int\left|E(\mathbf{r})\right|^{2}dA\int\left|\Psi(\mathbf{r})\right|^{2}dA},$
(1)

where $\mathbf{r}$ is the coordinate in the plane transverse to the beam at
the fiber. We use Eqn. 1, often referred to as the overlap integral, to
compute $\eta_{p}$ and $\eta_{s}$ by plugging in the field at the fiber,
$E(\mathbf{r})$, due to the planet and star, respectively. When stellar photon
noise dominates, $S/N\propto\eta_{p}/\sqrt{\eta_{s}}$, while in most other
cases $S/N\propto\eta_{p}$. In essence, the PIAA optics are designed for the
latter scenario, whereas the MDA is optimized for the former.

The PIAA optics consist of a pair of beam shaping lenses designed to alter the
effective distribution of light at the pupil, in order to reshape the PSF in
the focal plane. In the original implementation, PIAA lenses were used to
reduce the Airy rings in coronagraphs (Guyon, 2003). Later, PIAA optics were
used to optimize the coupling to a SMF (Jovanovic et al., 2017b).

Figure 3: Theoretical, monochromatic beam shapes and PSFs for the Keck
telescope (Top row) without apodization, (Middle row) with the PIAA lenses,
and (Bottom row) with the MDA. (Column 1) Beam intensity before focusing onto
the SMF. (Column 2-4) Log base-10 PSF intensity for (Column 2) an on-axis
source as well as sources with an angular separation of (Column 3)
5$\lambda/D$ and (Column 4) 10$\lambda/D$ from the optical axis (white cross).
All are normalized to the peak intensity in the non-apodized case and using
the corresponding central design wavelength (3.8 $\mu$m for the PIAA and
2.2$\mu$m for the MDA).

The KPIC PIAA optics consist of two aspheric lenses fabricated from CaF2
because of its high transmission across the wavelength regime that KPIC
operates ($K$ and $L$ bands; 2.0-4.2 $\mu$m). AR coatings for this spectral
region were optimized and deposited on both sides of each optic to minimize
Fresnel reflections. Figure 1(a) shows a ray diagram of the PIAA lenses, which
are designed to re-distribute an annular intensity distribution, similar to
the centrally-obscured Keck pupil, into a quasi-Gaussian profile. As a result,
the PSF is a better match to the fundamental mode (LP01) of a SMF and,
therefore, the coupling efficiency is improved.

The lenses each consist of one flat surface (facing outwards) and one aspheric
surface (facing inwards). As the light propagates from the first aspheric
surface to the second 50 mm away, the collimated beam is remapped and
recollimated. In this way the PIAA optic pair can be inserted into a
collimated beam without modifying the downstream focusing optics. The sag
profiles of the two lenses were designed by solving the differential equation
in Guyon (2003), which leverages Snell’s law to achieve the desired remapping.
Figure 1(b) and (c) show the fabricated lens pair while (d) and (e) show the
cross-section of the rotationally-symmetric sag profiles. We chose the sag
profiles to achieve the desired remapping function at the central wavelength
in the $L^{\prime}$ band ($\lambda_{0}$ = 3.8 $\mu$m) because we expect the
$S/N$ to be limited by thermal background in the longest wavelength filters.
However, the PIAA can be used over the full KPIC range with slightly degraded
coupling efficiency in the shorter wavelength filters.

One limitation of the PIAA optics is that the remapped intensity depends on
the position of the source with respect to the optical axis. While the lenses
are designed to provide a quasi-Gaussian beam shape for an on-axis source, the
exiting wavefront deviates from this and becomes heavily aberrated when the
source is off-axis by even a few resolution elements (i.e. $\lambda/D$, where
$\lambda$ is the wavelength and $D$ is the beam diameter at the pupil). Figure
3 shows the beam shape and corresponding PSFs without any alteration (Fig. 3,
top row) and with (Fig. 3, middle row) the PIAA optics. The PIAA modifies the
PSF favorably in the case of an on-axis source by concentrating more light
into the PSF core. However, the off-axis PSFs suffer from a strong coma
aberration with a magnitude that increases with the the angular separation
(Fig. 3 shows offsets of 5 and 10 $\lambda/D$). Due to the rotational symmetry
of the PIAA design, the direction of the coma aberration will always point
radially away from the core.

To utilize PIAA optics for HDC, the planet must be aligned to the optical axis
of the PIAA lenses to receive the coupling boost. The star would then be off
axis, and for the separations depicted in the figure, would diffract a small
amount of light onto the location of the planet marked by the white cross.
This is where its important to note that Fig. 3 displays the intensity of the
PSF, which does not necessarily indicate the amount of light that will couple
into a SMF. So although there may be star light of non-zero intensity at the
location of the planet fiber, we need to compute the coupling of that speckle
pattern to the SMF to understand the stellar leakage term. Figure 4 shows an
azimuthally averaged line profile taken from a 2D coupling map computed for
the PIAA, MDA and also with no optic in the beam (labelled Keck in the
figure). It can be seen that the off-axis coupling response with the PIAA is
very similar to that with a clear pupil, emphasizing that the PIAA does not
increase stellar leakage.

Where the PIAA makes use of transparent aspheric lenses, the MDA uses an array
of opaque microdots to shape the pupil illumination to reduce the brightness
of diffraction features around the PSF (Martinez et al., 2009a). Since the MDA
can be placed at, or near, the reimaged telescope pupil, the PSF is spatially-
invariant; i.e. the PSF is approximately the same for the planet and star (see
Fig. 3, bottom row). The MDA for KPIC uses a semi-random pattern of 25 $\mu$m
squares in a 200 nm thick chrome layer on an AR-coated CaF2 substrate (seen in
Fig. 2). The size of the microdots was chosen to maintain a $\sim$10:1 ratio
between the mask feature size and the wavelength of light to avoid anomalous
diffraction effects (Martinez et al., 2009b; Zhang et al., 2018). This design
is lossy because the chrome dots block light, creating a grayscale beam
intensity distribution. The mask is optimized by maximizing the $S/N$ for a
planet spectrum assuming the measurement is dominated by stellar photon noise.
For computational convenience, rather than maximizing the
$\eta_{p}/\sqrt{\eta_{s}}$ ratio, we minimize $\eta_{s}/\eta_{p}^{2}$, which
is proportional to the required exposure time to achieve a given $S/N$ in a
stellar photon noise limited regime.

To optimize the effective transmittance of the MDA, we assumed a polynomial
radial apodization function with field amplitude

$A(r)=\sum_{j=0}^{N}c_{j}(r/a)^{j},$ (2)

where $c_{j}$ are constants, $r$ is the radial coordinate, and $a$ is the beam
radius. In our simulations, we multiplied $A(r)$ by the expected Keck pupil
function, which is not rotationally-symmetric. However, using a one-
dimensional apodization function simplifies the optical system by not
requiring a rotational alignment between the apodizer and the beam. For each
set of $c_{j}$ coefficients, we computed the fraction of on-axis planet light
that couples into the SMF, normalized to the total energy incident on the
apodizer. This slightly modified definition of $\eta_{p}$ includes both losses
in coupling efficiency and the semi-transparent design of grayscale mask. The
effective $\eta_{s}$ was defined as the fraction of starlight that couples
into the on-axis SMF when the star is imaged off-axis (see Fig. 4), averaged
over planet-star separations of 3-15 $\lambda/D$. Using these metrics, we
found the optimal apodization function by minimizing $\eta_{s}/\eta_{p}^{2}$
with a simplex algorithm (see Fig. 3, bottom left). We repeated the
optimization for various values of $N$ and opted to use $N=5$ since adding
more terms did not significantly improve performance.

Figure 4: Azimuthal average of the coupling efficiency of an off-axis PSF
shown in Fig. 3 at the central design wavelength and assuming no other
aberrations. A coupling of 100% means that all of the light of the PSF couples
into the on-axis fiber. We can see that the PIAA does not significantly affect
the off-axis coupling of the clear Keck pupil. We can also see that the MDA
decreases the off-axis coupling from 3-12 $\lambda$/D, excluding two localized
regions of 4.5-6 and 8.5-9.

We converted the continuous apodization function into a binary microdot
pattern using the Floyd-Steinberg error diffusion algorithm (Floyd &
Steinberg, 1976; Dorrer & Zuegel, 2007; Zhang et al., 2018) taking into
account the relative transmission of the 200 nm thick chrome layer and the
exposed substrate. Figure 2 shows the designed pattern and a photo of the
fabricated MDA, which has approximately 500 microdots across the expected 12.3
mm beam size. The full coated region is 15 mm in diameter. Given that the
$S/N$ is more likely to be stellar photon noise limited at the shorter
wavelengths, we optimized the MDA for the central wavelength of $K$ band
($\lambda_{0}$ = 2.2 $\mu$m), but the MDA can technically be used with any
KPIC filter. At the longest wavelengths, where thermal background tends to
dominate, the MDA is not likely to improve performance because of the reduced
throughput and potentially high emissivity of the opaque regions.

Figure 5: Azimuthal average of the on-axis PSFs shown in Fig. 3 at the central
design wavelength and assuming no aberrations. The intensities in all cases
are normalized to the peak of the Keck PSF. The top panel shows the intensity
at the PSF core on a linear scale and the bottom panel shows the diffracted
intensity further from the source on a log scale.

Figure 5 shows azimuthally averaged line profiles of each on-axis PSF,
normalized to the non-apodized case. The PIAA increases the flux in the core
of the PSF, which leads to improved fiber coupling. On the other hand, the MDA
reduces the diffraction from the star at angular separations beyond 3
$\lambda/D$, but at the cost of adding the semi-transparent optic and lowering
the coupling efficiency. These technologies provide the ability to enhance the
performance of KPIC in two different regimes. However, each has a significant
downside. While the PIAA increases the throughput for an on-axis source, the
source and PIAA must be carefully aligned to the optical axis of the system;
the off-axis PSF becomes heavily aberrated. The MDA suppresses diffracted
starlight and thereby reduces stellar photon noise in the measurement of the
planet spectrum, but at the cost of overall throughput.

Figure 6: Image of the experimental lab bench in the PSF imaging mode. The red
line indicates the beam path for the PSF imaging mode. The orange line
indicates the beam path for fiber coupling mode after the objects circled in
orange have been inserted into the beam path. Inset in the bottom left corner
is the pupil mask emulating the Keck telescope pupil.

##  3 Experimental characterization

In this section we describe the experimental setup and tests carried out to
characterize the fabricated devices.

The key properties to examine for the two technologies are the PSF shape, the
throughput, and the coupling efficiency. To specify, we use throughput in this
paper to refer to the throughput of the individual optics and not the overall
effect that the optic has on the system’s throughput. The throughput of the
individual optics is less than 100% as a result of Fresnel reflections,
absorption, and reflections from the micro-dots in the case of the MDA.
Testing the throughput will help ensure that there were no errors in the
manufacturing of the optics and that the AR coatings are performing according
to spec. In addition, we will measure what effect these technologies have on
the light coupled to a SMF on/off axis and compare this to simulations.

###  3.1 Experimental Setup

Figure 6 shows the experimental setup with the PIAA mechanics in the beam. The
setup to characterize and measure the MDA is equivalent, but with the KPIC
coronagraph module at the “Test Optic” location.

This testbed can be separated into 4 main components. First, there is the
“telescope simulator,” which collimates light injected by a SMF. The SMF can
be swapped between a 2 $\mu$m laser and a blackbody light source that is
single mode across the 2-5 $\mu$m wavelength region. The collimating optic is
an AR coated CaF2 lens, with f=200 mm. Its distance from the light source can
be adjusted with a translation stage to ensure that the light is collimated
regardless of the wavelength. The beam then passes through a pupil mask, which
imitates the Keck pupil. Given the long focal length of the collimating lens,
the beam can be approximated as having a flat top illumination after passing
through the mask.

The second component to the testbed is the test optic. These are located in
the beam immediately following the telescope simulator in collimated space.
Both the PIAA and MDA are mounted on stages that can repeatably move the
optics in/out of the beam path. This allows for comparative measurements to be
taken quickly between the apodized beam and the non-apodized beam.

Following the test optics is the third component, the coupling arm. This
begins with a field-steering pickoff mirror (FSM1) that directs light toward
the off-axis parabola (OAP) used to couple the light. FSM1 is mounted on a
magnetic plate that can be removed and repeatably replaced to give quick
access to the two main functionalities of the testbed. This turning mirror is
complemented by a second, fixed steering mirror (FSM2). The two in combination
can translate the collimated beam and adjust its angle of incidence with
respect to the OAP (inside of the OAP block). This allows for coupling to be
efficiently optimized into the SMF. The OAP is a protected gold coated mirror
with an off-axis focal length of 36.6 mm, which injects the light into a ZBLAN
fiber (Le Verre Fluoré, MFD = $7.5~{}\mu$m at a $2.0~{}\mu$m wavelength, core
diameter = $6.5~{}\mu$m, cladding diameter = $125~{}\mu$m). With the pupil
mask diameter of 12.3 mm, this creates an $F\\#=3.0$, which closely matches
the measured ZBLAN fibers numerical aperture of 0.175 $\pm$ 0.01.

Lastly, the setup also contains an imaging arm. The imaging arm has two modes,
with the first taking the collimated output of the test optic and passing it
through a f=1000 mm AR coated CaF2 focusing lens (FL1) to create an image.
This arm is folded for convenience by a 2” silver mirror (FSM3) toward the
camera. The second imaging mode is re-imaging the output of the coupling fiber
onto the camera. The mounted fiber output is incident on a f=5.7 mm AR coated
aspheric lens (FL2), which collects all of the light from the fiber and re-
images it onto the camera. As seen in Fig. 6, this reimaging assembly can be
inserted and removed when switching from coupling mode to imaging mode.

An InSb camera (Merlin, Indigo) is used to carry out all imaging. It is
sensitive from 1-5 $\mu$m, and needs to be liquid nitrogen cooled. It has 240
x 318 pixels and a pixel pitch of 30 $\mu$m. The PSF formed by the $F\\#=81.3$
beam incident from FL1 is sampled by approximately 6 pixels per FWHM on this
detector. Immediately in front of the camera, we place a spectral filter, used
to both limit the background incident on the detector and to select the
measurement band when using the blackbody light source. When we use the 2
$\mu$m laser and do Ks measurements, a Ks-band filter (Asahi-Spectra) is used.
For L band measurements, we use a 3-5 micron bandpass filter.

###  3.2 Procedures

Before the camera could be used to make quantitative measurements, it required
us to linearize its response, maximize its dynamic range and do a one-point
non-uniformity correction (NUC) flat fielding. The camera’s exposure time and
video offset were coarsely adjusted to match the expected flux levels for the
experiments. Then, they were carefully tuned and the response of the camera to
varying flux levels was compared to the reading of a power meter. This
resulted in a linear camera response with the maximum possible dynamic range.
The next step was to perform a one-point NUC to ensure a uniform response
across the detector. This was done with the camera entrance window blocked by
a black piece of metal, such that a uniform thermal background from the plate
filled the camera’s entire field-of-view, while running the NUC correction
algorithm built into the camera.

####  3.2.1 PSF Imaging

We recognize that the shape of the PSF is critical and that the faint wings
can only be revealed at a high SNR. With the beam of interest projected onto
the camera, the exposure time and video offset were adjusted both to bring the
background counts down to around 50-100 ADUs and to maximize the signal. The 1
point NUC was applied as outlined above. The light source was blocked and a
cube of 100 frames that constitutes the background were collected. The light
source was unblocked and a cube of 100 frames with signal was collected.
Finally, the light source was blocked once again and another cube with 100
background frames was collected. The two sets of background frames, bracketing
the signal frames, were used to construct a background that was as similar to
the background when the signal was acquired as possible.

####  3.2.2 Throughput

To make quantitative measurements of the flux using the camera, aperture
photometry/radiometry was conducted. To obtain the throughput of a given
optical element (e.g. PIAA or MDA), the flux from the PSF generated when the
optic was in the beam was normalized by the flux of the PSF when the optic was
out of the beam. As such, it was important to keep the exposure time and other
camera settings constant between the two measurements to assure an accurate
comparison.

We take the same pattern of 100 background frames, 100 signal frames, and then
another 100 background frames both when the optics were in and out of the
beam.

####  3.2.3 Coupling Efficiency

To determine the coupling efficiency to the ZBLAN optical fiber, we again
compare the flux between two states: 1) light transmitted through the SMF re-
imaged onto the detector and 2) the PSF that was directly incident on the
fiber before being imaged onto the detector. These two states have non-common
optics and so the throughputs of those optics had to be accounted for during
the analysis.

To photometrically measure the coupling efficiency in a given band, we:

  1. 1.

Adjust the light source power, exposure time, and video offset to maximize the
signal from the re-imaged coupling fiber,

  2. 2.

Image, with backgrounds, the re-imaged fiber,

  3. 3.

Carefully remove the re-imaging system and FSM1 to reconfigure into PSF
imaging mode,

  4. 4.

Image, with backgrounds, the PSF without adjusting any settings.

After normalizing out the differential losses in each arm, we can take the
ratio between the extracted fluxes and determine the coupling efficiency in
isolation.

###  3.3 Data Acquisition and Analysis

There are multiple things we can do with the imaged PSF to analyze and
quantify it’s properties in comparison to the model. For example, one can
compute the overlap integral (Eqn 1) between the model and the experimental
PSFs to assess their similarity. For all other forms of analysis (i.e.
throughput or coupling), we must first extract the fluxes from the detector by
carefully removing the background. The PSF is localized to, at most, a region
with a radius of 25 pixels. Meanwhile the thermal background is a non-zero
bias across the frame, that can be assumed to be uniform following the 1 point
NUC. When using the single-mode blackbody source, the signal is small, so this
background subtraction must be done carefully. To accurately only extract the
flux from the PSF, we follow a two step process: 1) subtract a background
frame collected close in time to the data and 2) model out the remainder of
the background.

To achieve step 2, we take measurements of the total enclosed aperture flux at
increasing aperture radii. Beyond the extent of the PSF, the encircled flux
increases or decreases linearly (for a uniform background) with the total
enclosed pixels, $F_{internal}=F_{PSF}+\alpha\cdot n_{pix}$, where $\alpha$ is
the average background flux. Thus, in order to get $F_{PSF}$, we fit the
linear part of the curve and solve to determine the y-intercept, which is the
background free flux. The uncertainty in this measurement is the relevant
value from the covariance matrix of the fit.

With the flux’s extracted, the throughput can next be calculated. It is the
ratio of the fluxes: $\xi=F_{in}/F_{out}$

The coupling efficiency is mathematically very similar to the throughput, the
ratio of the flux coupled into the mode of the SMF and the flux incident on
the fiber face. However, there is more involved with getting to this step. As
was stated in the procedure, there are many non-common optics between the
fiber re-imaging block and the focusing arm of the bench. Once those are
accounted for, the coupling efficiency is as well: $\eta=F_{coup}/F_{in}$

The key to qualifying the experimental measurements is to be able to compare
them to accurate models. To ensure we get the highest degree of agreement, we
replicate the parameters of Fig. 6 in a simulator. There are uncertainties in
some of the physical parameters in the coupling arm of the setup: the NA of
the coupling fiber and the effective focal length of the coupling OAP are
known only to $\sim 5\%$. These uncertainties were accounted for in
simulations and resulted in a range of expected values.

##  4 Results

Here we summarize the measurements of the parameters outlined in the methods
section and compare them to simulations.

###  4.1 PSF profiles

The critical feature for both technologies is the structure of the PSFs and
specifically the faint structure in the wings. Only the 2 $\mu$m laser had
enough power to reveal the faint structure on the MIR camera. Figure 7 and 8,
present the PSFs obtained from the Indigo camera side-by-side with the
corresponding simulations and the residual difference between the two. This
enables qualitative comparison between the two, which is important to building
confidence in both the manufacturing and design process.

Figure 7: PSFs for the PIAA shown in both linear (top) and a base-10
logarithmic (bottom) space with the same spatial extent. (Left) Experimental,
(Center) Simulated, (Right) Residual of difference between the model
subtracted from the experimental PSF. The linear scaling residuals emphasize
the high degree of agreement between the model and experiment, while the
logarithmic scaling shows the structural differences between the two. All PSFs
are shown at $2~{}\mu$m.  Figure 8: PSFs for the MDA shown in both linear
(top) and a base-10 logarithmic (bottom) space with the same spatial extent.
(Left) Experimental, (Center) Simulated, (Right) Residual of difference
between the model subtracted from the experimental PSF. The linear scaling
residuals emphasize the high degree of agreement between the model and
experiment, while the logarithmic scaling shows the structural differences
between the two. All PSFs are shown at $2~{}\mu$m.

It can be seen from both figures, that there is a strong similarity between
the experimental and simulated PSFs. The right hand panel in both figures
shows that the residuals of the difference between the model and the
experimental PSFs are very low, with a peak around the $10^{-2}$ or 1%
intensity level. Further, we can also use the overlap integral as a
quantitative comparison tool for the PSFs. To account for misalignment’s in
the peaks of the two PSFs, we scan the model PSF across the experimental PSF
in 2D and compute an overlap integral map as shown in Fig. 9. From the overlap
integral maps, we can see that both technologies have a peak greater than
99.8%, indicating a high level of consistency between the experimental and
simulated PSFs, validating the design and manufacturing process.

Figure 9: Overlap integral maps between experimental and modelled PSFs for
(Left) PIAA and (Right) MDA. The figures show a high level of consistency
between the imaged PSFs and the simulated ones.

We can also compare the azimuthally averaged line profiles of the PSFs. By
over plotting the line profile of the experimental PSF with the simulated PSF,
we can identify locations where the flux is unexpectedly high or low. In Fig.
10, we show these line profiles of the two technologies. The shaded regions
indicate the standard deviation in the azimuthal direction of the imaged PSF.
We can see that both the PIAA and MDA closely follow the simulated curves,
once again emphasizing the similarity between the design and the characterized
devices.

Figure 10: Azimuthally averaged line profile of the a) PIAA PSF and b) MDA
PSF. The shaded region accounts for the uncertainty in the measurements.

###  4.2 Throughput

The throughput is defined as $\xi=F_{out}/F_{in}$, where $F_{in}$ is the flux
at the input to the test optic and $F_{out}$ is the flux at the output of the
test optic. The simulations were carried out in discrete monochromatic bands,
while the experimental measurements, aside from the 2 $\mu$m laser, were taken
using a broadband black-body light source combined with bandpass filters. The
simulated and measured throughputs are shown in Fig. 11, where the horizontal
error bars represent the bandwidth of the lab filters.

As can be seen from the figure, both technologies met the expected throughput
to within uncertainty. The PIAA is designed to reshape the pupil in a lossless
fashion as described in Section 2. The throughput is well above 90% in both K
and L bands as expected, and any systematic losses can be easily explained
through residual reflections across the four AR coated surfaces.

We can see that the MDA had a throughput consistent with simulations of $\sim
50\%$ in K-band. However, we see it begin to roll off in the experimental
results in the L band. This is due to the AR coating on the MDA not being
optimized for L band, because the losses associated with the MDA reduce the
L-band flux to the point where the thermal background dominates. Therefore,
the MDA was not expected to be scientifically useful in L band and hence the
coating was not optimized in this range.

Figure 11: A comparison of throughput curves between simulations and our
experimental measurements for the PIAA and MDA optics.

###  4.3 Coupling Efficiency

The simulations for the coupling efficiency were also produced in discrete
monochromatic bands, and compared to the polychromatic results obtained with a
filtered black-body light source in the laboratory. The coupling efficiencies
for the two technologies are shown in Fig. 12. The expected coupling for the
non-apodized pupil is overlaid on these figures for reference.

It can be seen that the coupling for the PIAA optics matches our simulation
well. We can also see that the coupling for the MDA matches the simulation
quite well, except for the measurement with the 2 $\mu$m laser. We are unsure
why this data point is off, but believe it has to do with the increased
coherence of the 2 $\mu$m laser and interference effects generated by the MDA.
To ensure that this was not merely an effect of the camera, the measurement
was repeated by placing a Thorlabs S148c power meter near the final focal
plane and making the same flux comparison as the camera. The power meter
confirmed the higher than expected coupling when using the MDA and the 2
$\mu$m laser beam. While this is not well understood, it is not an undesirable
result in the context of KPIC.

Figure 12: A comparison of coupling efficiency curves between the simulation
and our experimental measurements for both test optics.

From the figure it is clear that the PIAA boosts the coupling from $\sim 60\%$
in K and L bands without any additional optics to $\sim 70\%$ in K and $84\%$
in L band. In comparison, the MDA reduces the coupling in K band to $\sim
48\%$ and in L band to $45\%$. It should be made clear that these results are
for the Keck like pupil mask in our experiments, which has oversized spiders,
and at the Keck telescope, the coupling will be $3-4\%$ greater in all cases.

Looking at the impact of the optics at the system level by combining the
throughput and coupling efficiency, the PIAA has a total system throughput of
$66\%$ in K band and $77\%$ in L band while the MDA is closer to $\sim 24\%$
in both K and L bands. The impact of these efficiencies will be discussed in
the following section.

##  5 Discussion

It is clear from the results that the PIAA and MDA have different
characteristics that will be useful in different contexts, so we will address
them individually.

###  5.1 Beam Shaping Optics

The PIAA are lossless beam-shaping optics. We can see from the Ks and L band
measurements that there is an increase in coupling efficiency, from $\sim
60\%$ without the PIAA to $\sim 75\%$ with PIAA in K band and $\sim 85\%$ in L
band, where the lenses were optimized to operate. This boost in coupling is
only present for the on-axis source as PIAA optics suffer strong aberrations
off-axis as shown in Section 2 above. The lenses have a small loss ($\sim
8\%$), which is equal for both the star and the planet. For KPIC to utilize
these optics effectively, the known planet will need to be aligned with the
optical axis of the PIAA optics, i.e, on-axis. This will result in a
throughput increase for the planet light (improved coupling), and a reduction
in the starlight (losses from optics and aberrations for off-axis star).
Although the PIAA optics may scatter some of the star light onto the location
of the planet fiber, for small angular separations between the star/planet as
shown in Fig. 3, the coupling to the fiber will be extremely low. In this way,
the PIAA optics can reduce the exposure time to reach a given SNR, which will
also allow fainter objects to be observed.

###  5.2 Micro-dot Apodizer

The MDA’s effect is more subtle. The role of the MDA is to decrease the amount
of contaminating light from an off-axis source that gets coupled into the SMF.
As can be seen from the line profiles (Fig. 5), the PSF from the apodizer has
nearly an order of magnitude less flux between 3 and 12 $\lambda$/D compared
to the PSF from a non-apodized pupil. So, in situations where the light from
the on-axis host star is drowning out the signal from the off-axis planet,
using the MDA can help isolate the planet light. One drawback of the MDA is,
to achieve such starlight suppression, the apodizer must decrease the
throughput of the on-axis (star) and off-axis (planet) objects to $\sim 50\%$.
Furthermore, the coupling is $48\%$, lower than the native PSF and
significantly lower than when using the PIAA optics. Therefore, one must
carefully consider when to use this technology.

###  5.3 Exposure Time Calculation

Figure 13: Exposure time contours for S/N=10 detection of 51 Eri b using the
different technologies in K-band and L-band in units of hours. Left, we can
see in the K-band how the MDA decreases the planet’s throughput, but the
resultant increase in star contrast creates a net improvement in detection
exposure time. Right, we can see how the increased throughput of the PIAA
helps to decrease the detection exposure time.

In this section we look at the impact of the two technologies on the exposure
time of a simulated observing scenario of the planetary system 51 Eridani b,
with KPIC, in both K and L band. Although the star is typically on-axis during
acquisition, a tip/tilt mirror in the injection unit is adjusted in the final
steps of acquisition to put the known planet on-axis and align it with the
SMF. In this scenario, the exposure time ($\tau_{exp}$) of an observation to
reach a given SNR is proportional to the coupling efficiency of the stars PSF
($\eta_{s}$) and approximately inversely proportional to the square of the
coupling efficiency of the planet PSF ($\eta_{p}$):
$\tau_{exp}\propto\eta_{s}/\eta_{p}^{2}$

To see how each technology effects the exposure time, the properties of KPIC
were included in a realistic exposure time calculator (ETC). $\eta_{s}$ and
$\eta_{p}$ were scanned across a grid of reasonable values for the case of
observing 51 Eridani b. Figure 13 shows the contour maps of exposure time vs
raw contrast and throughput to achieve a cross-correlation function (CCF)
S/N=10. On these figures, discrete points can be seen for the case of the
native PSF (red dots), the PIAA (white dots) and the MDA (gold dots) PSFs. In
addition, the cyan dots highlight the phase I level of performance of KPIC on-
sky before deploying these technologies in phase II.

It can be seen that in the case of 51 Eridani b, the PIAA optics offer little
improvement in K-band compared to the expected phase II performance, while the
MDA can reduce the exposure time by $\sim 33\%$ from 3 to 2 hrs. Despite the
fact the MDA decreases the planet’s throughput, it also suppresses the star
light, improving contrast and the net effect is an improvement in the
detection exposure time. However, in L-band, which is dominated by thermal
background, the PIAA optics offer the advantages in exposure time, reducing
the exposure time of the expected phase II system by $\sim 33\%$ from 3 to 2
hrs, while the MDA should be avoided. Factors of $33\%$ are significant when
it comes to observing time on 8-10 m class telescopes.

This demonstrates the relative strengths and weaknesses of the two
technologies, and their applicability to a given observing scenarios depends
on the waveband of observation, the flux ratio between the star/planet and
their separation.

##  6 Conclusion

We have demonstrated the design, fabrication, characterization and accurate
simulation of both PIAA and MDA optics in the context of HDC. The optics were
optimized to operate in the K and L bands, which overlaps with the operating
spectral bands of the KPIC instrument. The PSFs, throughput and coupling
efficiencies all matched expectations from simulations, providing a high level
of confidence in both the simulation tools and manufacturing capabilities. In
a simulated observing scenario of 51 Eridani b, it was determined that the MDA
would offer a reduction in exposure time of $\sim 33\%$ in K-band, while the
PIAA would offer a similar level of reduction for L-band observations. Such
reductions in exposure time are important on 8-10 m class telescopes and are
necessary to be able to target even fainter objects in future. The application
of a particular technology to a given observation will depend on waveband of
choice, the flux ratio of the star/planet and their angular separation, and
needs to be determined on a case-by-case basis.

This work was supported by the Heising-Simons Foundation through grants
#2019-1312, #2017-318 and #2015-129. Support was provided by the Simons
Foundation grant ”Planetary Context of Habitability and Exobiology” G. Ruane
was supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship
under award AST-1602444. W. M. Keck Observatory is operated as a scientific
partnership among the California Institute of Technology, the University of
California, and the National Aeronautics and Space Administration (NASA). The
Observatory was made possible by the generous financial support of the W. M.
Keck Foundation. The authors wish to recognize and acknowledge the very
significant cultural role and reverence that the summit of Maunakea has always
had within the indigenous Hawaiian community. We are most fortunate to have
the opportunity to conduct observations from this mountain. Part of this work
was carried out at the Jet Propulsion Laboratory, California Institute of
Technology, under contract with NASA.

## References

  * Barman et al. (2015) Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61, doi: 10.1088/0004-637X/804/1/61
  * Beuzit et al. (2019) Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155, doi: 10.1051/0004-6361/201935251
  * Bottom et al. (2016) Bottom, M., Femenia, B., Huby, E., et al. 2016, in Proc. SPIE, Vol. 9909, Adaptive Optics Systems V, 990955
  * Crass et al. (2020) Crass, J., Bechter, A., Sands, B., et al. 2020, MNRAS, doi: 10.1093/mnras/staa3355
  * Currie et al. (2017) Currie, T., Guyon, O., Tamura, M., et al. 2017, ApJ, 836, L15, doi: 10.3847/2041-8213/836/1/L15
  * Dorrer & Zuegel (2007) Dorrer, C., & Zuegel, J. D. 2007, J. Opt. Soc. Am. B, 24, 1268, doi: 10.1364/JOSAB.24.001268
  * Floyd & Steinberg (1976) Floyd, R. W., & Steinberg, L. 1976, Proc. Soc. Inf. Disp., 17, 75
  * Groff et al. (2016) Groff, T. D., Chilcote, J., Kasdin, N. J., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, 99080O, doi: 10.1117/12.2233447
  * Guyon (2003) Guyon, O. 2003, A&A, 404, 379, doi: 10.1051/0004-6361:20030457
  * Jovanovic et al. (2015) Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890, doi: 10.1086/682989
  * Jovanovic et al. (2017a) Jovanovic, N., Guyon, O., Kotani, T., et al. 2017a, arXiv e-prints, arXiv:1712.07762.  https://arxiv.org/abs/1712.07762
  * Jovanovic et al. (2017b) Jovanovic, N., Schwab, C., Guyon, O., et al. 2017b, A&A, 604, A122, doi: 10.1051/0004-6361/201630351
  * Jovanovic et al. (2019) Jovanovic, N., Delorme, J.-R., Bond, C. Z., et al. 2019, in Techniques and Instrumentation for Detection of Exoplanets IX, ed. S. B. Shaklan, Vol. 11117, International Society for Optics and Photonics (SPIE), doi: 10.1117/12.2529330
  * Kasdin et al. (2003) Kasdin, N. J., Vanderbei, R. J., Spergel, D. N., & Littman, M. G. 2003, Astrophys. J., 582, 1147
  * Macintosh et al. (2014) Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the National Academy of Science, 111, 12661, doi: 10.1073/pnas.1304215111
  * Males et al. (2018) Males, J. R., Close, L. M., Miller, K., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10703, Adaptive Optics Systems VI, 1070309, doi: 10.1117/12.2312992
  * Martinez et al. (2009a) Martinez, P., Dorrer, C., Aller Carpentier, E., et al. 2009a, Astron. Astrophys., 495, 363, doi: 10.1051/0004-6361:200810918
  * Martinez et al. (2009b) Martinez, P., Dorrer, C., Kasper, M., Boccaletti, A., & Dohlen, K. 2009b, Astron. Astrophys., 500, 1281, doi: 10.1051/0004-6361/200911824
  * Mawet et al. (2017) Mawet, D., Ruane, G., Xuan, W., et al. 2017, ApJ, 838, 92, doi: 10.3847/1538-4357/aa647f
  * Nisenson & Papaliolios (2001) Nisenson, P., & Papaliolios, C. 2001, Astrophys. J. Lett., 548, L201.  http://stacks.iop.org/1538-4357/548/i=2/a=L201
  * Sayson et al. (2019) Sayson, J. L., Ruane, G., Mawet, D., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 019004, doi: 10.1117/1.JATIS.5.1.019004
  * Schwab et al. (2012) Schwab, C., Leon-Saval, S. G., Betters, C. H., Bland-Hawthorn, J., & Mahadevan, S. 2012, arXiv e-prints, arXiv:1212.4867.  https://arxiv.org/abs/1212.4867
  * Snellen et al. (2015) Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59, doi: 10.1051/0004-6361/201425018
  * Soummer et al. (2003) Soummer, R., Aime, C., & Falloon, P. E. 2003, Astron. Astrophys., 397, 1161, doi: 10.1051/0004-6361:20021573
  * Sparks & Ford (2002) Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543, doi: 10.1086/342401
  * Vigan et al. (2018) Vigan, A., Otten, G. P. P. L., Muslimov, E., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, 1070236, doi: 10.1117/12.2313681
  * Wang et al. (2017) Wang, J., Mawet, D., Ruane, G., Hu, R., & Benneke, B. 2017, AJ, 153, 183, doi: 10.3847/1538-3881/aa6474
  * Watson et al. (1991) Watson, S. M., Mills, J. P., Gaiser, S. L., & Diner, D. J. 1991, Appl. Opt., 30, 3253, doi: 10.1364/AO.30.003253
  * Woillez et al. (2003) Woillez, J., Guerin, J., Perrin, G. S., et al. 2003, in Interferometry for Optical Astronomy II, ed. W. A. Traub, Vol. 4838, International Society for Optics and Photonics (SPIE), 1341 – 1343, doi: 10.1117/12.459362
  * Zhang et al. (2018) Zhang, M., Ruane, G., Delorme, J.-R., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, 106985X, doi: 10.1117/12.2312831


# Federated Learning Algorithms for Generalized Mixed-effects Model (GLMM) on
Horizontally Partitioned Data from Distributed Sources

Wentao Li  Jiayi Tong  Md. Monowar Anjum  Noman Mohammed  Yong Chen  Xiaoqian
Jiang

###### Abstract

Objectives: This paper developed federated solutions based on two
approximation algorithms to achieve federated generalized linear mixed effect
models (GLMM). The paper also proposed a solution for numerical errors and
singularity issues. And compared the developed model’s outcomes with each
other, as well as that from the standard R package (‘lme4’).

Methods: The log-likelihood function of GLMM is approximated by two numerical
methods (Laplace approximation and Gaussian Hermite approximation), which
supports federated decomposition of GLMM to bring computation to data. To
solve the numerical errors and singularity issues, the loss-less estimation of
log-sum-exponential trick and the adaptive regularization strategy was used to
tackle the problems caused by federated settings.

Results: Our proposed method can handle GLMM to accommodate hierarchical data
with multiple non-independent levels of observations in a federated setting.
The experiment results demonstrate comparable (Laplace) and superior
(Gaussian-Hermite) performances with simulated and real-world data.

Conclusion: We modified and compared federated GLMMs with different
approximations, which can support researchers in analyzing versatile
biomedical data to accommodate mixed effects and address non-independence due
to hierarchical structures (i.e., institutes, region, country, etc.).

###### keywords:

GLMM , Federated learning , Mixed effects , Laplace approximation , Gauss-
Hermite approximation

††journal: Journal of Biomedical Informatics

[inst1]organization=School of Biomedical Informatics,
UTHealth,addressline=7000 Fannin St, city=Houston, postcode=77030, state=TX,
country=USA

[inst2]organization=Department of Biostatistics, Epidemiology and Informatics,
University of Pennsylvania,addressline=3400 Civic Center Boulevard,
city=Philadelphia, postcode=19104, state=PA, country=USA

[inst3]organization=Department of Computer Science, University of
Manitoba,city=Winnipeg, country=Canada

##  1 Introduction

There is an increasing surge of interest in analyzing biomedical data to
improve health. Biostatisticians and machine learning researchers are keen to
access personal health information for a deeper understanding of diagnostics,
disease development, and potential preventive or treatment options [1].

In the US, healthcare and clinical data are often collected by local
institutions. For many situations, combining these datasets would increase
statistical power in hypothesis testing and provide better means to
investigate regional differences and subpopulation bias (e.g., due to
differences in disease prevalence or social determinants). However, such an
information harmonization process needs to respect the privacy of individuals,
as healthcare data contain sensitive information about personal
characteristics and health conditions. As a minimum requirement, HIPAA (Health
Insurance Portability and Accountability Act)[2] specifies PHIs (protected
health information) and regulations to de-identify the sensitive information
(i.e., safe harbor mechanism). But HIPAA compliance does not mean full
protection of the data, as several studies demonstrated re-identifiability of
HIPAA de-identified data [3, 4, 5]. Ethical healthcare data sharing and
analysis should also respect the “minimum necessary” principle to reduce the
unnecessary risk of potential data leakage, which might increase the
likelihood of information leakage.

Figure 1: Schema of federated learning model in multiple geographically
distributed healthcare institutions. The local institutions periodically
exchange intermediate statistics and update the convergence situation of the
global model.

The recent development of federated learning, which intends to build a shared
global model without moving local data from their host institutions (Fig. 1),
shows good promise in addressing the challenge in data sharing mentioned
above. Despite the exciting progress, there is still an important limitation
as existing models cannot effectively handle mixed-effects (i.e., both fixed
and random effects), which is very important to analyzing non-independent,
multilevel/hierarchical, longitudinal, or correlated data. Also, due to the
sampling errors (i.e., smaller sample size in local sites), variances from
these local statistics are larger than those of the global model. These
issues, if not addressed appropriately, would lead to failure in global
optimization. The goal of this paper is to improve existing techniques and
provide practical solutions with open-source implementation and to allow
ordinary biomedical/healthcare researchers to build federated mixed effect
learning models for their studies.

##  2 Related Work

Federated learning for healthcare data analysis is not a new topic, and there
have been many previous studies. However, most of the existing methods assume
the observations are independent and identically distributed[6, 7]. In the
presence of non-independence due to hierarchical structures (e.g., due to
institutional or regional differences), existing federated models have strong
limitations in ignoring the regional differences. The generalized linear mixed
model (GLMM), which takes the heterogeneous factors into consideration, is
more amenable to accommodate the heterogeneity across healthcare systems.
There have been very few studies in this area and one relevant work is a
privacy-preserving Bayesian GLMM model[8], which proposed an Expectation-
Maximization (EM) algorithm to fit the model collaboratively on horizontally
partitioned data. The convergence process is relatively slow (due to the
Metropolis–Hastings sampling in the E-step) and it is also not very stable
(likely to be trapped in local optima[9] in high-dimensional data). In the
experiment, a loose threshold (i.e., 0.08) was used as a convergence condition
[8] while typical federated learning algorithms [10] in healthcare use much
stringent convergence threshold (i.e., $10^{-6}$).

Another related work to fit GLMM in a federated manner is the distributed
penalized quasi-likelihood (dPQL) algorithm[11]. This algorithm reduces the
computational complexity by considering the target function of penalized
quasi-likelihood, which is motivated from Laplacian approximation. The model
has communication efficiency over the EM approach and can converge in a few
shots. However, the target function PQL can have first order asymptotic bias
[12] due to the Laplacian approximation of the integrated likelihood.

There is an alternative strategy, Gauss-Hermite (GH), which supports high-
order approximation. It is computationally more intensive and requires special
techniques to handle the numerical instability of the logSumExp operation (due
to the overflow issue when the dimensionality grows in the sum of the
exponential terms). We will explain both models in this manuscript and compare
their performance on simulated and real-world data.

##  3 Methods

In this section, we will discuss the statistic model along with challenges to
be tackled. A high-level schema of the method is shown in algorithm 1.

###  3.1 Notation

Before we introduce the formation of GLMM, let us define some notations.

$i$ | Index of sites | $l_{i}$ | Log-likelihood function for site $i$
---|---|---|---
$j$ | Index of patients in a specific site | $\bm{\beta}$ | Parameters of fixed effect
$k$ | Index of Hermite polynomial | $\mu_{i}$ | Parameters of random effect in site $i$
$K$ | Order of Hermite polynomial | $\tau$ | Hyper-parameters
$m$ | Number of sites | $\bm{\theta}$ | Parameter space $(\bm{\beta},\tau)$
$n_{i}$ | Number of patients in site $i$ | $X_{ij}$ |  | A vector represents the data of $j$-th
---
patient in $i$-th site
$\mathcal{L}_{i}$ | Likelihood function for site $i$ | $y_{ij}$ | The outcome of patient $j$ from site $i$
$\lambda$ | The parameter of regularization term | $p$ | Number of variables

###  3.2 Fitting GLMM with quasi-likelihood

Let us provide the formation of the GLMM. Define $\mathbb{P}$ is the
distribution of interest and depending on patient-level data $X_{ij},\
y_{ij}$. Define $\phi$ as the distribution of random effects. We can compose
the joint distribution as following

$\prod_{j=1}^{n_{i}}\mathbb{P}(\bm{\theta}|X_{ij},y_{ij})\phi(\mu_{i};\tau)$

Now we have the log-likelihood function of the joint distribution:

$\log\\{\mathcal{L}(\bm{\theta})\\}=\sum_{i=1}^{m}\log\left\\{\int_{\mu_{i}}\left[\prod_{j=1}^{n_{i}}\mathbb{P}(\bm{\theta}|X_{ij},y_{ij})\right]\phi(\mu_{i};\tau)\text{d}\mu_{i}\right\\}$
(1)

From the log-likelihood function Eq. (1), one can see that it does not support
direct linear decomposition. In order to support federated learning, we will
leverage approximation strategies to make the objective linearly decomposable
with simple summary statistics.

We will compare Laplace approximation and Gauss-Hermite approximation in the
following sections.

###  3.3 Laplace (LA) approximation

With the help of Laplace approximation, the integration from Eq.(1) can be
approximated by an exponential family expression.

$\int_{\mu_{i}}f_{\theta}(\mu_{i})\text{d}\mu_{i}=\int_{\mu_{i}}e^{\log{f_{\bm{\theta}}(\mu_{i})}}\text{d}\mu_{i}\triangleq\int_{\mu_{i}}e^{g(\mu_{i},\theta)}\text{d}\mu_{i}$
(2)

After the deduction (A.1), the intractable problem is solved and the objective
is to maximize the following formula with respect to $\theta$, where $g$ is an
exponential family function defined above (Eq.(2))

$\sum_{i=1}^{n_{i}}\left(g(\hat{\mu}_{i},\bm{\theta})-\dfrac{n_{i}}{2}\log\left(g_{\mu\mu}(\hat{\mu}_{i},\bm{\theta})\right))\right)$

, for which the terms are linearly decomposable from local sites. Site $i$
needs to calculate the following aggregated data:

  * •

$p\times p$ matrix:

$\dfrac{\hat{\omega}_{\beta\beta}\hat{\omega}-\hat{\omega}_{\beta}\hat{\omega}_{\beta}}{\hat{\omega}^{2}}+\hat{\mu}_{\beta\beta}g_{\mu}+\hat{\mu}_{\beta}(\hat{\mu}_{\beta}g_{\mu\mu}+g_{\mu\beta})+\hat{\mu}_{\beta}g_{\mu\beta}+g_{\beta\beta}$
(3)

  * •

$p$ \- dim vector:

$\dfrac{\hat{\omega}_{\beta}}{\hat{\omega}}+\hat{\omega}^{2}g_{\mu\beta}(\hat{\mu}_{i})g_{\mu}+g_{\beta}$
(4)

  * •

scalar of random effect: $\hat{\mu}_{i}$ and first order derivative of $\tau$
by

$\dfrac{\hat{\omega}_{\tau}}{\hat{\omega}}+\hat{\omega}^{2}g_{\mu\tau}(\hat{\mu}_{i})g_{\mu}+g_{\tau}$

where $\hat{\omega}=\sqrt{-\dfrac{1}{g_{\mu\mu}(\hat{\mu}_{i0})}}$

###  3.4 Gauss-Hermite (GH) approximation

Gauss-Hermite approximation [13] implements Hermite interpolation concerning
Eq. (2).  And after the deduction in A.2, notice that when the order of
Hermite polynomial $K=1$, the objective function is identical to the method
with Laplace approximation. Because GH is more generalizable, we will describe
the distributed federated learning model on the GLMM problem with the
formation of Gauss-Hermite approximation Eq. (9). For each site $i$, the
followings need to calculate and transmit:

  * •

$p\times p$ matrix:

$\dfrac{\hat{\omega}_{\beta\beta}\hat{\omega}-\hat{\omega}_{\beta}\hat{\omega}_{\beta}}{\hat{\omega}^{2}}+\dfrac{1}{\sum_{k=1}^{K}f_{k}}\sum_{k=1}^{K}\dfrac{\partial}{\partial\bm{\beta}}(f_{k_{\mu}}\hat{\mu}_{\beta}+f_{k_{\omega}}\hat{\omega}_{\beta}+f_{k_{\beta}})-\dfrac{1}{(\sum_{k=1}^{K}f_{k})^{2}}\|\sum_{k=1}^{l}(f_{k_{\mu}}\hat{\mu}_{\beta}+f_{k_{\omega}}\hat{\omega}_{\beta}+f_{k_{\beta}})\|^{2}_{2}$
(5)

  * •

$p$ \- dim vector:

$\dfrac{\hat{\omega}_{\beta}}{\hat{\omega}}+\dfrac{1}{\sum_{k=1}^{K}f_{k}}\sum_{k=1}^{K}(f_{k_{\mu}}\hat{\mu}_{\beta}+f_{k_{\omega}}\hat{\omega}_{\beta}+f_{k_{\beta}})$
(6)

  * •

scalar of random effect: $\hat{\mu}_{i}$ and first order derivative of $\tau$
by

$\dfrac{\hat{\omega}_{\tau}}{\hat{\omega}}+\dfrac{1}{\sum_{k=1}^{K}f_{k}}\sum_{k=1}^{K}(f_{k_{\mu}}\hat{\mu}_{\tau}+f_{k_{\omega}}\hat{\omega}_{\tau}+f_{k_{\tau}})$

###  3.5 Training Penalization GLMM with GH approximation

The convergence of the approximation of the likelihood function may be
compromised due to over-fitting. Also, for those spatially correlated data,
the convergence of them may lead to a complex model. Hence, $L2$
regularization is added to the local log-likelihood function of Gauss-Hermite
approximation form, and as shown below

$l_{i}=\log\mathcal{L}_{i}=\log\left(\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\theta})+x_{k}^{2}\right\\}\right)-\lambda\|\bm{\beta}\|^{2}_{2}$
(7)

note that when $K=1$, it is represented as regularized Laplace approximation
to the problem. To evaluate and find the optimum $\lambda$, we steadily
increased the value of $\lambda$ in range $[0,10]$ by $1$. Set
$\lambda_{\text{opt}}$ as the optimized regularization term with largest
$\sum_{i}^{m}l_{i}$. And choose $\hat{\bm{\beta}}_{\text{opt}}$ as the
optimized estimator for $\bm{\beta}$.

Due to the limited computation digits, computers are not able to calculate the
correct results of the local log-likelihood function $l_{i}$ of the Gauss-
Hermite approximation form as stated above. Such problem is also known as the
Log-Sum-Exponential problem and can be solved by shifting the center of the
exponential sum for easier computation,

$\log\sum_{k=1}^{K}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\theta})+x_{k}^{2}\right\\}=a+\log\sum_{k=1}^{K}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\theta})+x_{k}^{2}-a\right\\}$

where $a$ is an arbitrary number.

Thus, the global problem of maximizing $\sum_{i}^{m}l_{i}$ can be divided into
several local maximization problems (7). Each local site $i$ will update the
regression intermediates, and they will be combined to update the iteration
status. Specifically, in each iteration of the federated GLMM algorithm, the
following statistics are exchanged from each site to contribute aggregated
data for the global model

LA | GH
---|---
number of variables $p$ | number of variables $p$
$p\times p$ matrix (Eq. 3) |  $p\times p$ matrix (Eq. 5)
$p$ \- dim vector (Eq. 4) |  $p$ \- dim vector (Eq. 6)
$p$ \- dim vector $\bm{\beta}$ |  $p$ \- dim vector $\bm{\beta}$
scalar $\lambda$ | scalar $\lambda$
scalar $\hat{\mu}_{i}$ | scalar $\hat{\mu}_{i}$
scalar first order derivative $\tau$ | scalar first order derivative $\tau$
| scalar $K$

Detailed derivatives with the logistic regression setting of the optimization
are presented in the appendix A.3.

Data: Local data $X_{i}$ from site $i$

Result: Global model with coefficients $\hat{\bm{\beta}}_{\text{global}}$,
test P-values, upper and lower bound

Initialization: coefficients $\bm{\beta}$, random effects $\mu_{i}$ and
$\tau$, regularization term $\lambda$;

for  _$\lambda=0$ to $10$_ do

while  _$\Delta\mu\neq 0$_ do

Maximized $\mu_{i}$ with respect to $\bm{\theta}$

while  _$\Delta\bm{\theta}\neq 0$_ do

1\. Approximate the log-likelihood functions $l_{\text{approx}}$ with Laplace
or Gauss-Hermite;

2\. Calculate intermediate statistics. For Laplace, Eq. 3, 4; for Gauss-
Hermite, Eq. 5, 6;

3\. Send the intermediate statistics to center server, and then the center
server will aggregate them;

4\. Update $\bm{\theta}$ and $l_{\text{approx}}({\bm{\theta}},\lambda)$ in
center server and send back to each client $i$;

end while

end while

end for

Return: The largest $l_{\text{approx}}(\hat{\bm{\theta}},\hat{\lambda})$, the
coefficients $\hat{\bm{\beta}}_{\text{global}}=\hat{\bm{\beta}}$,
hyperparameter $\hat{\tau}$, and the regularization term $\hat{\lambda}$

Algorithm 1 Distributed GLMM with approximation methods

##  4 Experiments

Our algorithm is developed in Python with packages pandas, numpy, scipy, and
the benchmark algorithm is glmer function in R package ‘lme4’.

###  4.1 Benchmarking the methods using synthetic data

To test the performance of our proposed methods, we first designed a stress
test based on a group of synthetic data, which include 8 different settings
(Tab.1), and each set contains 20 datasets. In each dataset, it consists of 4
categorical variables with value in $\\{0,1\\}$; 6 categorical variables with
value in range $[-1,1.5]\in\mathbb{R}$; 1 outcome variable with value in
$\\{0,1\\}$; Site ID, represents the id of which site the entry belongs to;
Site sample size, represents the number of samples in this specific setting;
Log-odds ratio for each sample; Number of true positive, true negative, true
positive, false positive, false negative.

To evaluate which method can reach better performance, we proposed the
following evaluation measurements: discrimination of the estimated
coefficients $\hat{\bm{\beta}}$, the test power of each coefficient, and the
precision and recall of the number of significant coefficients.

Table 1: The summary of data in each setting. Setting |  | Number of
---
sites
| Sample size
---
in each site
variance
1 | 2 | 500 | small
2 | 2 | 500 | large
3 | 10 | 500 | small
4 | 10 | 500 | large
5 | 2 | 30 | small
6 | 2 | 30 | large
7 | 10 | 30 | small
8 | 10 | 30 | large

The valuation experiments were conducted among federated GLMM with Laplace
approximation, federated GLMM with Gauss-Hermite approximation, and
centralized GLMM (all of the data stored in single host) in the R package. And
the stress test will be run in 160 different datasets in 8 different settings
as mentioned in Tab.1. All of the data in different settings were randomly
separated into training sets and validation sets with a ratio of 7:3. And we
trained the federated learning model on training data sets, then by slowly
increasing the regularization term $\lambda$, we chose the optimum model with
the best Akaike information criterion and Bayesian information criterion
performance on the validation sets. All testing was performed on 2017 iMac
with 16 GB memory, CPU (4.2 GHz Quad-Core Intel Core i7), macOS Big Sur
version 11.6, Python 3.8, and R version 3.5.0.

Although we tested the data sets with the state-of-art benchmark algorithm for
centralized GLMM in R, the regression is not perfect for the ground truth
coefficients we used to generate the data (Fig.2). So, it is also important to
have the P-values of variables into consideration when interpreting the model.
Thus, We made comparisons among centralized GLMM, Laplace method, and Gauss-
Hermite method concerning the p-values of coefficients. Tables in the appendix
captured the performance of different methods. Fig.3 shows the precision and
recall results of centralized, Laplace, and Gauss-Hermite methods. Noted that
we set our Gauss-Hermite approximation to 2-degree. See tables in Appendix
(Tab.5, 6, 7).

Figure 2: The difference from coefficients to the true parameters that are
used to generate data. (Left) The distributed GLMM with Laplace approximation;
(Middle) The distributed GLMM with 2-degree Gauss-Hermite approximation.
Reminds that $X_{1}$ is the intercept; (Right) The benchmark of centralized
GLMM in R package.

Figure 3: The precision and recall among centralized, Laplace, and Gauss-Hermite method under significance level $\alpha=0.05$. (Left) The precision of the test compared to the true value. (Right) The recall of the test compared to the true value.  Figure 4: The curve of test power among centralized, Laplace, and Gauss-Hermite methods. (Left) The power of the test of the Laplace method. (Middle) The power of the test of the 2-degree Gauss-Hermite method. (Right) The power of the test of the Centralized method. Power was calculated as the two-sided t-test on p-values among different methods. Table 2: The convergence rates on approximation methods LA and GH. (Both LA and GH held the same convergence threshold $10^{-3}$. The mean values and standard deviations (in parentheses) were given) | LA | GH
---|---|---
Setting | Steps | Runtime (s) | Steps | Runtime (s)
1 | 22.875 (21.623) | 47.953 (20.513) | 34.850 (9.213) | 104.460 (10.614)
2 | 21.500 (21.977) | 40.947 (36.466) | 35.000 (8.711) | 100.940 (19.940)
3 | 29.867 (31.719) | 108.931 (65.486) | 34.900 (6.138) | 1259.285 (231.956)
4 | 27.846 (24.034) | 84.343 (76.502) | 36.650 (6.310) | 1342.695 (250.603)
5 | 59.722 (42.057) | 10.631 (3.945) | 33.750 (10.146) | 12.568 (2.116)
6 | 67.188 (48.994) | 10.499 (4.054) | 31.400 (11.081) | 11.430 (3.064)
7 | 96.286 (53.635) | 96.501 (38.632) | 37.450 (3.818) | 369.165 (41.998)
8 | 116.083 (46.479) | 91.304 (62.410) | 37.150 (4.295) | 309.693 (36.621)

The simulation results showed the federated Gauss-Hermite approximation
performed better than the method based on Laplace approximation on every
variable. Also, the federated Gauss-Hermite method achieved higher test power
(Fig.4). When considering the convergence rates between the two approximation
methods, both showed less convergence efficiency in Setting 7 and 8 (Tab.2).
The result Indicates that more local sites and smaller sample sizes will make
the federated GLMM more inefficient to converge. Also, GH approximation method
will required more computation time compared with LA approximation. In sum,
one-degree increase of the approximation function in LA with our developed GH
method, GH outperformed LA methods for federated GLMM implementation.

###  4.2 Mixed-effects logistic regression on mortality for patients with
COVID-19

We analyzed the data of COVID-19 electronic health records collected by Optum®
from February 2020 to January 28, 2021, from a network of healthcare
providers. The dataset has been de-identified and based on HIPAA statistical
de-identification rules and managed by Optum® customer data user agreement. In
this database, there are 56,898 unique positive tested COVID-19 patients.
After removing the patients with missing data, the final cohort contains 4,531
patients who died and the rest population (41,781) survived. The database
contains a regional variable with five levels (Midwest, Northwest, South,
West, Others/unknown) to provide privacy-preserving area information to
indicate where the samples were collected.

We have conducted a GLMM model (considering region-distinct random effect)
using this dataset with the following predictors: age, gender, race,
ethnicity, Chronic obstructive pulmonary disease (COPD), Congestive heart
failure (CHF), Chronic kidney disease (CKD), Multiple sclerosis (MS),
Rheumatoid arthritis (RA), LU (other lung diseases), High blood pressure
(HTN), ischemic heart disease (IHD), diabetes (DIAB), Asthma (ASTH), obesity
(Obese). Our proposed method with GH approximation performed the best with
both the smallest Akaike information criterion (AIC) and Bayesian information
criterion (BIC) according to the table of the goodness of fit (Tab.3). And the
performance of different methods can be shown in (Tab.4).

Table 3: Statistics of goodness of fit among different methods | Log-likelihood | AIC | BIC
---|---|---|---
R | 13562.9 | 27165.9 | 27340.8
LA | -13695.0 | 27428.0 | 27594.1
GH | -11.8 | 61.6 | 227.7
Figure 5: The ROC curve with Area Under Curve (AUC) among centralized, Laplace, and Gauss-Hermite methods. The orange ROC curve is the centralized method without regularization and the Laplace approximation(i.e., R implementation in the ‘lme4’ package, which does not have an option for including regularization). AUC values are also included, a higher AUC value implicates better performance of the model. The green ROC curve is the 2-degree Gauss-Hermite method with regularization. Table 4: Statistics of performances among different methods(95% CIs were generated by Wilson Score interval) |  | Precision | Recall | F1-score | AUC | threshold
---|---|---|---|---|---|---
Centralized method with LA | Value | 0.1507 | 0.6204 | 0.2425 | 0.6789 | 0.0900
Lower bound (0.95) | 0.1474 | 0.6160 | 0.2386 | 0.6700 |
Upper bound (0.95) | 0.1539 | 0.6248 | 0.2464 | 0.6878 |
GH with regularization | Value | 0.1705 | 0.6546 | 0.2705 | 0.7178 | 0.0108
Lower bound (0.95) | 0.1670 | 0.6503 | 0.2664 | 0.7091 |
Upper bound (0.95) | 0.1739 | 0.6589 | 0.2745 | 0.7265 |

We also compared the ROC curves (Fig.5) between our proposed GH method and
centralized method to check their performance. And the result showed that GH
approximation (AUC=0.72) outperforms the centralized method without
regularization (AUC=0.68). Indicating GH-based GLMM method has better
classification performance than the GLMM based on LA approximation. In our
proposed model, it showed variables: Unknown race, Chronic kidney disease
(CKD), Multiple sclerosis (MS), and other lung diseases (LU) are not
significant to the mortality of COVID-19. The result of the regression is in
the Appendix (Tab.8, 9, 10).

##  5 Discussion

We developed solutions to address the limited digit problem (i.e., overflow
issue of fixed-length object types due to extremely large numbers in local
estimation) using an alternative loss-less estimation of log-sum-exponential
term, and the singularity issue (involved in Newton optimization) with an
adaptive regularization strategy to avoid inverting low-rank matrices without
imposing too much unnecessary smoothness. We further compared two federated
GLMM algorithms with our developed federated solutions (LA vs. GH) and
demonstrated the performance of the federated GLMM based on the GH method
surpassed the method based on LA in terms of the accuracy of estimation, power
of tests, and AUC. Although the GH method is requiring slightly more
computations than the LA method, it is still acceptable for more accurate
results. For example, in the prediction of COVID-19 mortality rates, the
accuracy of prediction will be more reliable, as we have shown in the previous
section. During the optimization iterations, we noticed that some sites have
already achieved convergence in very few steps. If those sites stop
communicating with the central server, they can be released from extra
computations. We would investigate more efficient algorithms based on such a
strategy of ‘lazy regression’ for minimizing communication for federated
learning models.

Another limitation of the proposed federated GLMM model is not yet
differentially private and iterative summary statistics exchange can lead to
incremental information disclosure, which might increase the re-identification
risk over time. There are several strategies to improve the model based on
secure operations like homomorphic encryption and differential privacy, which
we have previously studied in GLM models [14]. Finally, in practice, there can
be extra heterogeneity that cannot be explained by random intercepts only, it
is of interest to further develop our algorithms toward GLMM that allows
multiple random effects including random coefficients in the regression
models.

##  6 Acknowledgement

XJ is CPRIT Scholar in Cancer Research [RR180012]; and he was supported in
part by Christopher Sarofim Family Professorship, UT Stars award, UTHealth
startup, the National Institute of Health (NIH) under award number
[R01AG066749, R01GM114612, U01TR002062]; and the National Science Foundation
(NSF) RAPID [#2027790]. NM is supported by NSERC Discovery Grants
[RGPIN-2015-04147]. JT and YC’s research is supported in part by NIH under
award number [R01AI130460 and R01LM012607] and Patient-Centered Outcomes
Research Institute (PCORI) Project Program Award (ME-2019C3-18315). All
statements in this report, including its findings and conclusions, are solely
those of the authors and do not necessarily represent the views of the
Patient-Centered Outcomes Research Institute (PCORI), its Board of Governors
or Methodology Committee.

## References

  * [1] A. Telenti, X. Jiang, Treating medical data as a durable asset, Nature Genetics 52 (10) (2020) 1005–1010.
  * [2] Hipaa privacy rule, 45 c.f.r. s 164 (2014).
  * [3] L. Bonomi, X. Jiang, Linking temporal medical records using non-protected health information data, Statistical methods in medical research 27 (11) (2018) 3304–3324.
  * [4] V. Janmey, P. L. Elkin, Re-identification risk in hipaa de-identified datasets: The mva attack, in: AMIA Annual Symposium Proceedings, Vol. 2018, American Medical Informatics Association, 2018, p. 1329.
  * [5] L. Sweeney, J. S. Yoo, L. Perovich, K. E. Boronow, P. Brown, J. G. Brody, Re-identification risks in hipaa safe harbor data: A study of data from one environmental health study, Technology science 2017 (2017).
  * [6] J. A. Nelder, R. W. Wedderburn, Generalized linear models, Journal of the Royal Statistical Society: Series A (General) 135 (3) (1972) 370–384.
  * [7] N. R. Draper, H. Smith, Applied regression analysis, Vol. 326, John Wiley & Sons, 1998.
  * [8] R. Zhu, C. Jiang, X. Wang, S. Wang, H. Zheng, H. Tang, Privacy-preserving construction of generalized linear mixed model for biomedical computation, Bioinformatics 36 (Supplement_1) (2020) i128–i135.
  * [9] L. Xu, M. Jordan, On convergence properties of the EM algorithm for Gaussian mixtures, Neural computation 8 (1) (1996) 129–151.
  * [10] Y. Wu, X. Jiang, J. Kim, L. Ohno-Machado, G rid binary lo gistic re gression (glore): building shared models without sharing data, Journal of the American Medical Informatics Association 19 (5) (2012) 758–764.
  * [11] C. Luo, M. N. Islam, N. E. Sheils, J. Buresh, Y. Chen, dpql: a lossless distributed algorithm for generalized linear mixed model with application to privacy-preserving hospital profiling, medRxiv (2021).
  * [12] X. Lin, N. E. Breslow, Bias correction in generalized linear mixed models with multiple components of dispersion, Journal of the American Statistical Association 91 (435) (1996) 1007–1016.
  * [13] Q. Liu, D. A. Pierce, A note on gauss—hermite quadrature, Biometrika 81 (3) (1994) 624–629.
  * [14] M. Kim, J. Lee, L. Ohno-Machado, X. Jiang, Secure and differentially private logistic regression for horizontally distributed data, IEEE Transactions on Information Forensics and Security 15 (2019) 695–710.
  * [15] G. Broström, H. Holmberg, Generalized linear models with clustered data: Fixed and random effects models, Computational Statistics & Data Analysis 55 (12) (2011) 3123–3134.

##  Appendix A Supplementary proofs

###  A.1 Laplace approximation

Let us explain the Laplace approximation. Denote that

$f_{\theta}(\mu_{i})\triangleq\prod_{j=1}^{n_{i}}\mathbb{P}(\theta|X_{ij},y_{ij})\phi(\mu_{i};\tau)$

and one can see the log term inside the log-likelihood function that

$\int_{\mu_{i}}f_{\theta}(\mu_{i})\text{d}\mu_{i}=\int_{\mu_{i}}e^{\log{f_{\theta}(\mu_{i})}}\text{d}\mu_{i}\triangleq\int_{\mu_{i}}e^{g(\mu_{i},\theta)}\text{d}\mu_{i}$
(8)

Applying a Taylor expansion on $g(\mu_{i},\theta)$, and we choose
$\hat{\mu}_{i}$ that maximized $g(\mu_{i},\theta)$. See that $\hat{\mu}_{i}$
satisfies $g_{\mu}(\hat{\mu}_{i},\theta)=0$ and
$g_{\mu\mu}(\hat{\mu}_{i},\theta)<0$, we have

$g(\mu_{i},\theta)=g(\hat{\mu}_{i},\theta)-\dfrac{1}{2}(\hat{\mu}_{i}-\mu_{i})^{2}\left(-g_{\mu\mu}(\hat{\mu}_{i},\theta)\right)+o(\mu_{i}^{2})$

, which is plugged into Eq. (2). With Laplace approximation [15], one can see
that

$\displaystyle\int_{\mu_{i}}e^{g(\mu_{i},\theta)}\text{d}\mu_{i}\approx
exp\\{{g(\hat{\mu}_{i},\theta)}\\}\left[2\pi\cdot-\dfrac{1}{g_{\mu\mu}(\hat{\mu}_{i},\theta)}\right]^{n_{i}/2}$

###  A.2 Gauss Hermite approximation

The Hermite polynomial $H_{k}(x)$ and weight $h_{k}$ are defined as
followings,

$H_{k}(x)\triangleq(-1)^{k}e^{x^{2}}\dfrac{d^{k}}{dx^{k}}e^{-x^{2}}\quad\quad
h_{k}\triangleq\dfrac{2^{k-1}k!\sqrt{\pi}}{k^{2}\left[H_{k-1}(x_{k})\right]^{2}}$

where $x_{k}$ are the roots of $H_{k}(x)=0$.

Thus, with the Gauss-Hermite approximation, Eq. (2) can be approximated by

$\int_{\mu_{i}}e^{g(\mu_{i})}\text{d}\mu_{i}\approx\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k})+x_{k}^{2}\right\\},\quad\hat{\omega}=\sqrt{-\dfrac{1}{g^{\prime\prime}(\hat{\mu}_{i})}}$
(9)

notice that when $K=1$, it is a Laplace approximation.

Our final objective function for GH is

$\mathcal{L}_{i}=\mathcal{L}_{i}(\bm{\beta},\mu_{i};X_{i\cdot},y_{i\cdot})=\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\beta})+x_{k}^{2}\right\\}$

Denote that

$\displaystyle
f_{k}=h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\beta})+x_{k}^{2}\right\\}$
$\displaystyle f_{k_{\beta}}=f_{k}g_{\beta}(\hat{\mu}_{i})$ $\displaystyle
f_{k_{\mu}}=f_{k}g_{\mu}(\hat{\mu}_{i})$ $\displaystyle
f_{k_{\omega}}=f_{k}g_{\mu}(\hat{\mu}_{i})\sqrt{2\pi}x_{k}$

###  A.3 Optimization

A logistic regression model with random effects is developed under the form of

$\log\\{\mathcal{L}(\theta)\\}=\sum_{i=1}^{m}\log\left\\{\int_{\mu_{i}}\left[\prod_{j=1}^{n_{i}}\mathbb{P}(\theta|X_{ij},y_{ij})\right]\phi(\mu_{i};\tau)\text{d}\mu_{i}\right\\}$

and the distribution $\mathbb{P}$ follows density of logit and $\phi$ is a
univariate normal, see that

$\prod_{j=1}^{n_{i}}\mathbb{P}(\theta|X_{ij},y_{ij})=\prod_{j=1}^{n_{i}}\pi_{ij}^{y_{ij}}(1-\pi_{ij})^{(1-y_{ij})}$
$\phi(\mu_{i};\theta)=\dfrac{1}{\sqrt{2\pi}\tau}\exp(-\mu_{i}^{2}/2\tau^{2})$

where $\pi_{ij}$ is a Sigmoid function of $\mu_{i}$ and defined as

$\pi_{ij}=\dfrac{\exp(X_{ij}^{\top}\bm{\beta}+\mu_{i})}{1+\exp(X_{ij}^{\top}\bm{\beta}+\mu_{i})}$

with the Gauss-Hermite approximation set up, the objective function can be
approximated as

$\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\beta})+x_{k}^{2}\right\\}$

where

$\displaystyle g(\mu_{i};\bm{\beta})$
$\displaystyle=\log\prod_{j=1}^{n_{i}}\mathbb{P}(\theta|X_{ij},y_{ij})\phi(\mu_{i};\tau)$
$\displaystyle=\sum_{j=1}^{n_{i}}\left[\log\mathbb{P}(\theta|X_{ij},y_{ij})\right]+\log\phi(\mu_{i};\tau)$
$\displaystyle=\sum_{j=1}^{n_{i}}\left[y_{ij}\log\pi_{ij}+(1-y_{ij})\log(1-\pi_{ij})\right]+\log\phi(\mu_{i};\tau)$

####  A.3.1 Step 1: Maximize $g(\mu_{i})$

To maximize $g(\mu_{i})$, we need to get the derivatives

$\displaystyle\dfrac{\partial g}{\partial\mu_{i}}$
$\displaystyle=\sum_{j=1}^{n_{i}}\left[y_{ij}\dfrac{1}{\pi_{ij}}\dfrac{\partial\pi_{ij}}{\partial\mu_{i}}-(1-y_{ij})\dfrac{1}{1-\pi_{ij}}\dfrac{\partial\pi_{ij}}{\partial\mu_{i}}\right]+\dfrac{1}{\phi}\dfrac{\partial\phi}{\partial\mu_{i}}$
$\displaystyle=\sum_{j=1}^{n_{i}}(y_{ij}-\pi_{ij})-\dfrac{\mu_{i}}{\tau^{2}}$

where

$\dfrac{\partial\phi}{\partial\mu_{i}}=(\sqrt{2\pi}\tau)^{-1}\exp(-\mu_{i}^{2}/2\tau^{2})\cdot(-\mu_{i}/\tau^{2})$

and

$\displaystyle\dfrac{\partial^{2}g}{\partial\mu_{i}^{2}}$
$\displaystyle=-\sum_{j=1}^{n_{i}}\dfrac{\partial\pi_{ij}}{\partial\mu_{i}}-\dfrac{1}{\tau^{2}}<0$

where

$\dfrac{\partial\pi_{ij}}{\partial\mu_{i}}=\dfrac{\exp\left(X_{ij}^{\top}\bm{\beta}+\mu_{i}\right)}{\left[1+\exp\left(X_{ij}^{\top}\bm{\beta}+\mu_{i}\right)\right]^{2}}$

see that it is a convex problem, using newton’s method, we can derive
$\hat{\mu}_{i}=\arg\max_{\mu_{i}}g(\mu_{i})$ that is global optimum.

####  A.3.2 Step 2: Maximization preparation of $\beta$ in LOCAL

See the derivative

$\dfrac{\partial\pi_{ij}}{\partial\bm{\beta}}=\dfrac{X_{ij}\exp\left(X_{ij}^{\top}\bm{\beta}+\mu_{i}\right)}{\left[1+\exp\left(X_{ij}^{\top}\bm{\beta}+\mu_{i}\right)\right]^{2}}$

and denote
$f_{k}(\hat{\mu}_{i};\bm{\beta}):=h_{k}\exp\left\\{g(\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k};\bm{\beta})+x_{k}^{2}\right\\}$,
then

$\displaystyle\dfrac{\partial\mathcal{L}_{i}}{\partial\bm{\beta}}$
$\displaystyle=\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}\left\\{f_{k}(\hat{\mu}_{i};\bm{\beta})h_{k}\left.\dfrac{\partial
g(\mu_{i};\bm{\beta})}{\partial\bm{\beta}}\right|_{\mu_{i}=\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k}}\right\\}$
(10)
$\displaystyle=\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}\left\\{f_{k}(\hat{\mu}_{i};\bm{\beta})h_{k}\sum_{j=1}^{n_{i}}(X_{ij}y_{ij}-X_{ij}\pi_{ij})\right\\}$
(11)

and the second derivative

$\displaystyle\dfrac{\partial^{2}\mathcal{L}_{i}}{\partial\bm{\beta}^{2}}=\sqrt{2\pi}\hat{\omega}\sum_{k=1}^{K}$
$\displaystyle\left\\{f_{k}(\hat{\mu}_{i};\bm{\beta})h_{k}^{2}\sum_{j=1}^{n_{i}}(X_{ij}y_{ij}-X_{ij}\pi_{ij})\left[\sum_{j=1}^{n_{i}}(X_{ij}y_{ij}-X_{ij}\pi_{ij})\right]^{\top}\right.$
(12)
$\displaystyle\quad+\left.f_{k}(\hat{\mu}_{i};\bm{\beta})h_{k}\sum_{j=1}^{n_{i}}\left(-X_{ij}\dfrac{\partial\pi_{ij}}{\partial\bm{\beta}}\right)\right\\}$
(13)

Notice that $\mu_{i}$ in (3), (4) and (5) are replaced by
$\hat{\mu}_{i}+\sqrt{2\pi}\hat{\omega}x_{k}$ where $\hat{\mu}_{i}$ is the
maximand of function $g(\cdot)$ with respect to $\mu_{i}$.

####  A.3.3 Step 3: Maximization of $\beta$ in GLOBAL

Reminds that $\mathcal{L}=\sum_{i=1}^{m}\log\mathcal{L}_{i}$, then another
Newton’s method is applied in global log-likelihood function,

$\dfrac{\partial\mathcal{L}}{\partial\bm{\beta}}=\sum_{i=1}^{m}\dfrac{\mathcal{L}_{i}^{\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}\quad\quad\dfrac{\partial^{2}\mathcal{L}}{\partial\bm{\beta}^{2}}=\sum_{i=1}^{m}\left[\dfrac{\mathcal{L}_{i}^{\prime\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}-\left(\dfrac{\mathcal{L}_{i}^{\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}\right)^{2}\right]$

Now, focus on
$\bm{\beta}^{(n+1)}=\bm{\beta}^{(n)}-\dfrac{\mathcal{L}^{\prime}(\bm{\beta}^{(n)})}{\mathcal{L}^{\prime\prime}(\bm{\beta}^{(n)})}$,
deduce that

$\dfrac{\mathcal{L}^{\prime}(\bm{\beta}^{(n)})}{\mathcal{L}^{\prime\prime}(\bm{\beta}^{(n)})}=\dfrac{\sum_{i=1}^{m}\dfrac{\mathcal{L}_{i}^{\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}}{\sum_{i=1}^{m}\dfrac{\mathcal{L}_{i}^{\prime\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}-\sum_{i=1}^{m}\left(\dfrac{\mathcal{L}_{i}^{\prime}(\bm{\beta})}{\mathcal{L}_{i}(\bm{\beta})}\right)^{2}}$

###  A.4 Synthetic data generation

There are 8 settings of data sets generated from the process, and each setting
can be summarized in table 1.

Fig.6 shows the distribution of each setting about different variables.

Figure 6: The box-plots of variables in Setting 1 to Setting 8.

We set true sensitivity and specificity as $sen=0.6\text{ and }sp=0.9$ and
$\bm{\beta}=(-1.5,0.1,-0.5,-0.3,0.4,-0.2,-0.25,0.35,-0.1,0.5)$. Also define
$X_{1}=\mathbb{1}_{N}$ as the intercept, and $X_{2},X_{3},X_{4}$ are generated
with Bernoulli distribution (2) with probability $p=0.1,0.3,0.5$ respectively.

$f(X_{i};p)=\left\\{\begin{array}[]{ll}p&\text{ if }X_{i}=1\\\ q=1-p&\text{ if
}X_{i}=0\end{array}\right.$ (14)

then $X_{5},X_{6},X_{7}$ are generated from normal distributions
$\mathcal{N}(0,0.5),\mathcal{N}(0,1),\mathcal{N}(0,1.5)$ respectively. Lastly,
$X_{8},X_{9},X_{10}$ are generate from uniform distributions
$\mathcal{U}(-0.5,0.5)$, $\mathcal{U}(-0.7,0.7)$, $\mathcal{U}(-1,1)$
respectively. We also generate the random effect $\mu$ using trivariate normal
distribution

$\mathcal{N}_{3}\left(\left(\begin{array}[]{ccc}0\\\ 0\\\ 0\\\
\end{array}\right),\Sigma\right),\quad\Sigma=I_{3}$ (15)

and with the settings, we can deduce the log-odds ratio with following formula

$\log(\pi)=f(X\bm{\beta}+\mu+\epsilon)$ (16)

where $f$ is the sigmoid function defined as

$f(x)=\dfrac{e^{x}}{1+e^{x}}$ (17)

Now, we generate the outcomes $y$ for each sample with Bernoulli distribution
(2) where the log-odds ratio served as the probability $p$. Also, the
sensitivity and specificity can be calculated with Binomial distribution with
probability $sen+\mu_{2}$ and $sp+\mu_{3}$.

##  Appendix B Supplementary tables

Table 5: The performance of centralized GLMM on R package with a significance threshold $\alpha=0.05$. TNR refers to the True negative rate. | Precision | Recall | TNR | Accuracy
---|---|---|---|---
X1 | 0.7162 | 1.0000 | 0.0000 | 0.7162
X2 | 0.1000 | 0.9231 | 0.2000 | 0.2635
X3 | 0.5917 | 0.9726 | 0.3467 | 0.6554
X4 | 0.4400 | 0.9649 | 0.2308 | 0.5135
X5 | 0.5703 | 0.9865 | 0.2568 | 0.6216
X6 | 0.4797 | 1.0000 | 0.0000 | 0.4797
X7 | 0.7162 | 1.0000 | 0.0000 | 0.7162
X8 | 0.3167 | 0.9268 | 0.2336 | 0.4257
X9 | 0.1417 | 0.8500 | 0.1953 | 0.2838
Table 6: The performance of distributed GLMM with Laplace transformation with a significance threshold $\alpha=0.05$. TNR refers to the True negative rate. | Precision | Recall | TNR | Accuracy
---|---|---|---|---
X1 | 0.7063 | 1.0000 | 0.0000 | 0.7063
X2 | 0.0833 | 1.0000 | 0.2667 | 0.3125
X3 | 0.5750 | 0.9857 | 0.4333 | 0.6750
X4 | 0.4609 | 1.0000 | 0.3168 | 0.5688
X5 | 0.5303 | 0.9859 | 0.3034 | 0.6063
X6 | 0.4375 | 1.0000 | 0.0000 | 0.4375
X7 | 0.6563 | 1.0000 | 0.0000 | 0.6563
X8 | 0.3000 | 1.0000 | 0.3226 | 0.4750
X9 | 0.1500 | 1.0000 | 0.2817 | 0.3625
X10 | 0.5316 | 1.0000 | 0.0263 | 0.5375
Table 7: The performance of distributed GLMM with 2-degree Gauss-Hermite transformation with significance threshold $\alpha=0.05$. TNR refers to True negative rate. | Precision | Recall | TNR | Accuracy
---|---|---|---|---
X1 | 0.5455 | 1.0000 | 0.0000 | 0.5455
X2 | 0.3158 | 1.0000 | 0.3390 | 0.4935
X3 | 0.6404 | 1.0000 | 0.4938 | 0.7338
X4 | 0.6066 | 1.0000 | 0.4000 | 0.6883
X5 | 0.5952 | 1.0000 | 0.3544 | 0.6688
X6 | 0.5260 | 1.0000 | 0.0000 | 0.5260
X7 | 0.6948 | 1.0000 | 0.0000 | 0.6948
X8 | 0.5000 | 0.9828 | 0.4063 | 0.6234
X9 | 0.4386 | 1.0000 | 0.3846 | 0.5844
X10 | 0.5658 | 1.0000 | 0.0294 | 0.5714
Table 8: The result of centralized GLMM in R package | Coef | Std.Err | z | P-value | [0.025 | 0.975]
---|---|---|---|---|---|---
(Intercept) | -5.882 | 0.133 | -44.294 | 0.000 | -6.142 | -5.621
age | 0.043 | 0.001 | 30.918 | 0.000 | 0.040 | 0.045
Gen_M | 0.370 | 0.033 | 11.052 | 0.000 | 0.304 | 0.435
race_Asian | 0.365 | 0.122 | 2.995 | 0.003 | 0.126 | 0.604
race_Caucasian | 0.157 | 0.049 | 3.185 | 0.001 | 0.061 | 0.254
race_Other.Unknown | 0.385 | 0.070 | 5.530 | 0.000 | 0.249 | 0.521
ethnicity_Not.Hispanic | -0.181 | 0.060 | -3.030 | 0.002 | -0.299 | -0.064
ethnicity_Unknown | -0.088 | 0.069 | -1.274 | 0.203 | -0.224 | 0.048
COPD_Y | 0.096 | 0.038 | 2.556 | 0.011 | 0.022 | 0.169
CHF_Y | 0.153 | 0.040 | 3.796 | 0.000 | 0.074 | 0.232
CKD_Y | 0.029 | 0.044 | 0.674 | 0.500 | -0.056 | 0.115
MS_Y | -0.090 | 0.124 | -0.725 | 0.469 | -0.332 | 0.153
RA_Y | 0.144 | 0.070 | 2.048 | 0.041 | 0.006 | 0.281
LU_Y | -0.001 | 0.204 | -0.004 | 0.997 | -0.400 | 0.398
HTN_Y | 0.113 | 0.053 | 2.121 | 0.034 | 0.009 | 0.217
IHD_Y | 0.334 | 0.037 | 9.077 | 0.000 | 0.262 | 0.406
DIAB_Y | 0.149 | 0.035 | 4.214 | 0.000 | 0.080 | 0.218
ASTH_Y | -0.170 | 0.054 | -3.143 | 0.002 | -0.276 | -0.064
Obese_Y | 0.211 | 0.043 | 4.893 | 0.000 | 0.126 | 0.295
Table 9: The result of federated GLMM with GH method | Coef | Std.Err | z | P-value | [0.025 | 0.975]
---|---|---|---|---|---|---
(Intercept) | -3.476 | 0.066 | -53.064 | 0.000 | -3.605 | -3.348
age | 0.043 | 0.001 | 55.794 | 0.000 | 0.041 | 0.044
Gen_M | 0.370 | 0.021 | 17.847 | 0.000 | 0.329 | 0.410
race_Asian | 0.364 | 0.075 | 4.847 | 0.000 | 0.217 | 0.511
race_Caucasian | 0.156 | 0.028 | 5.538 | 0.000 | 0.101 | 0.211
race_Other.Unknown | 0.383 | 0.041 | 9.365 | 0.000 | 0.303 | 0.463
ethnicity_Not.Hispanic | -0.178 | 0.036 | -4.936 | 0.000 | -0.249 | -0.107
ethnicity_Unknown | -0.091 | 0.042 | -2.154 | 0.031 | -0.174 | -0.008
COPD_Y | 0.096 | 0.026 | 3.700 | 0.000 | 0.045 | 0.146
CHF_Y | 0.154 | 0.029 | 5.205 | 0.000 | 0.096 | 0.211
CKD_Y | 0.028 | 0.031 | 0.901 | 0.367 | -0.033 | 0.090
MS_Y | -0.093 | 0.078 | -1.205 | 0.228 | -0.245 | 0.059
RA_Y | 0.144 | 0.049 | 2.959 | 0.003 | 0.048 | 0.239
LU_Y | 0.001 | 0.119 | 0.008 | 0.994 | -0.233 | 0.235
HTN_Y | 0.113 | 0.028 | 3.988 | 0.000 | 0.057 | 0.169
IHD_Y | 0.334 | 0.024 | 14.162 | 0.000 | 0.288 | 0.381
DIAB_Y | 0.149 | 0.022 | 6.677 | 0.000 | 0.105 | 0.192
ASTH_Y | -0.169 | 0.032 | -5.242 | 0.000 | -0.233 | -0.106
Obese_Y | 0.211 | 0.028 | 7.483 | 0.000 | 0.156 | 0.267
Table 10: The result of federated GLMM with LA method | Coef | Std.Err | z | P-value | [0.025 | 0.975]
---|---|---|---|---|---|---
(Intercept) | -3.162 | 0.064 | -49.437 | 0.000 | -3.288 | -3.037
age | 0.041 | 0.001 | 54.895 | 0.000 | 0.040 | 0.043
Gen_M | 0.359 | 0.021 | 17.374 | 0.000 | 0.318 | 0.399
race_Asian | 0.315 | 0.075 | 4.214 | 0.000 | 0.168 | 0.461
race_Caucasian | 0.130 | 0.028 | 4.675 | 0.000 | 0.076 | 0.185
race_Other.Unknown | 0.330 | 0.041 | 8.138 | 0.000 | 0.251 | 0.410
ethnicity_Not.Hispanic | -0.211 | 0.036 | -5.886 | 0.000 | -0.281 | -0.140
ethnicity_Unknown | -0.114 | 0.042 | -2.707 | 0.007 | -0.197 | -0.031
COPD_Y | 0.098 | 0.026 | 3.750 | 0.000 | 0.047 | 0.149
CHF_Y | 0.156 | 0.030 | 5.234 | 0.000 | 0.098 | 0.215
CKD_Y | 0.029 | 0.032 | 0.928 | 0.353 | -0.033 | 0.092
MS_Y | -0.092 | 0.077 | -1.192 | 0.233 | -0.244 | 0.059
RA_Y | 0.138 | 0.049 | 2.829 | 0.005 | 0.042 | 0.234
LU_Y | -0.007 | 0.118 | -0.061 | 0.952 | -0.239 | 0.225
HTN_Y | 0.104 | 0.028 | 3.722 | 0.000 | 0.049 | 0.159
IHD_Y | 0.337 | 0.024 | 14.208 | 0.000 | 0.290 | 0.383
DIAB_Y | 0.144 | 0.022 | 6.461 | 0.000 | 0.100 | 0.187
ASTH_Y | -0.175 | 0.032 | -5.445 | 0.000 | -0.238 | -0.112
Obese_Y | 0.182 | 0.028 | 6.450 | 0.000 | 0.127 | 0.237