# HI–shielding of ${\bf H_{2}}$ in UV–irradiated protogalaxies: suppression of
the photodissociation rate

Meredith Neyer1,2 ID and Jemma Wolcott-Green1 ID
1Department of Physics, University of California Santa Barbara, MC 9530, Santa
Barbara, CA 93106, USA
2Department of Physics and Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the impact of neutral hydrogen absorption on ${\rm H_{2}}$
photodissociation in protogalactic haloes exposed to soft-UV radiation. Lyman-
series absorption can significantly deplete dissociating photons as line
overlap with the ${\rm H_{2}}$ Lyman-Werner bands occurs for neutral column
densities exceeding $10^{22}~{}{\rm cm^{-2}}$, but this effect has not been
previously included in studies of protogalactic haloes. We use high–resolution
three–dimensional hydrodynamic simulations to investigate this “HI–shielding”
in three metal–free atomic cooling haloes collapsing at redshift $z\sim
10-20$. We use cloudy modeling to update a previous fitting formula for
HI–shielding which is a better model for shielding of non–ground state ${\rm
H_{2}}$ rovibrational populations and implement the new fit in our
simulations. We find that the inclusion of HI–shielding increases the
“critical flux” for suppression of ${\rm H_{2}}$ cooling in these haloes by
$\sim 60-100$ per cent. The larger critical flux has implications in
particular for the predicted numbers of candidate haloes in which“direct
collapse” could seed massive ($\sim 10^{5}~{}{\rm M_{\odot}}$) black holes at
$z\sim 15$.

###### keywords:

cosmology: theory – early Universe – galaxies: formation – molecular processes
– stars: Population III

††pagerange: HI–shielding of ${\bf H_{2}}$ in UV–irradiated protogalaxies:
suppression of the photodissociation rate–References††pubyear: 2022

##  1 Introduction

Molecular hydrogen, ${\rm H_{2}}$, has been extensively studied in the context
of the first generation of stars and galaxies, in which it plays a crucial
role as the primary coolant of primordial gas below $\sim 10^{4}$K (for a
review, see Bromm & Yoshida, 2011). Prior to the production and dispersion of
metals by the supernovae, the thermodynamic evolution of pristine primordial
gas depends sensitively on the ${\rm H_{2}}$ abundance and therefore on the
photodissociation of ${\rm H_{2}}$, which occurs in the presence of soft UV
photons in the “Lyman–Werner” (LW) bands (11.1-13.6 eV).

Depletion of ${\rm H_{2}}$ by LW radiation has been shown to raise the minimum
mass of protogalactic haloes in which gas is able to condense and cool, thus
delaying star formation in smaller “minihaloes” (Haiman et al., 1997, 2000;
Machacek et al., 2001; Yoshida et al., 2003; Mesinger et al., 2006; Wise &
Abel, 2007; O’Shea & Norman, 2008; Kulkarni et al., 2020; Schauer+21). In more
massive haloes, with virial temperatures $\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{4}$K, cooling by
neutral hydrogen allows gas to condense in haloes even in the absence of
significant ${\rm H_{2}}$ cooling, rendering these “atomic cooling haloes”
(ACHs) less vulnerable to feedback from a cosmological background LW radiation
(e.g. Oh & Haiman, 2002). Typically, the column densities of ${\rm H_{2}}$ in
ACHs grow large enough that ${\rm H_{2}}$ becomes self–shielding: systematic
depletion of LW band photons in the outer layers of the halo depresses
photodissociation of ${\rm H_{2}}$ in the core, allowing the gas to cool to
temperatures of a few hundred Kelvin. However, sufficiently strong LW
radiation fields have been shown suppress the ${\rm H_{2}}$ abundance and
thereby to prevent gas in ACHs from cooling below the virial temperature of
the halo (see Inayoshi et al., 2020, and references therein). This threshold
LW flux strength is commonly referred to as the critical flux or “$J_{\rm
crit}$” and has been typically found in hydrodynamic simulations to be in the
range $10^{3-4}$ in the customary units $10^{-21}~{}{\rm
erg~{}s^{-1}~{}cm^{-2}~{}Hz^{-1}~{}sr^{-1}}$. While this is orders of
magnitude larger than the expected cosmological background (e.g. Dijkstra et
al., 2008), a collapsing halo near a particularly bright neighboring galaxy
with recently-formed Pop III stars may be exposed to a such a flux (Visbal et
al., 2014; Regan et al., 2017); in this ”synchronized collapse” scenario, if
the two collapse within a short period of time – of order a few Myr – the
second halo to cross the atomic cooling threshold may have ${\rm
H_{2}}$–cooling entirely suppressed.

The presence of a super–critical flux has implications for the formation of
massive seed black holes, $\sim 10^{(4-5)}M_{\odot}$; rapid accretion in ACHs
that remain near the virial temperature; these “heavy seeds” could help
explain the existence of the earliest supermassive black holes, observed to
have masses $\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{9}{\rm M_{\odot}}$ at redshifts
$z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}6$ and as high as $\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}7.5$ (Fan et al.,
2001; Fan et al., 2003; Morganson et al., 2012; Mazzucchelli et al., 2017;
Wang et al., 2019; Yang et al., 2019; Wang et al., 2021) These heavy seeds are
commonly referred to as “direct collapse” black holes, though they’re thought
to form via an intermediary supermassive star phase (e.g. Haemmerlé et al.,
2018).

Extensive work has been done to constrain the value of $J_{\rm crit}$ using
hydrodynamic simulations of ACHs, which relies on detailed modeling of the
${\rm H_{2}}$ chemistry. Since the fraction of haloes exposed to a
super–critical UV flux depends sensitively on $J_{\rm crit}$, even small
changes in the photodissociation rate significantly alters the predicted
prevalence of direct collapse halo candidates (Dijkstra et al., 2008; Ahn et
al., 2009; Agarwal et al., 2012; Dijkstra et al., 2014; Chon et al., 2016).

###  1.1 The ${\bf H_{2}}$ photodissociation rate

Self–shielding by ${\rm H_{2}}$ occurs as the LW bands become optically thick
at column densities $\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{13}~{}{\rm cm}^{-2}$, suppressing
the photodissociation rate (e.g. Draine & Bertoldi, 1996). The optically–thick
rate in general depends on the column density, gas temperature, rovibrational
populations of ${\rm H_{2}}$(Wolcott-Green & Haiman, 2019), and details of the
incident spectrum (Agarwal & Khochfar, 2014; Sugimura et al., 2014; Wolcott-
Green et al., 2017), and is prohibitively computational expensive to calculate
on–the–fly in simulations, due to the large number of LW transitions.
Simulations most often therefore implement a fitting formula to model the
optically–thick rate and rely on local estimates of the column density
(Wolcott-Green et al. 2011, but see Hartwig et al. 2015).

In addition to self–shielding, absorption of LW photons by neutral hydrogen
Lyman series resonances can decrease the rate of ${\rm
H_{2}}$–photodissociation. Processing of the cosmological UV background by HI
in the pre–reionization IGM has been studied in detail (Haiman et al., 1997,
2000); however, the effects of HI absorption within protogalactic halos has
not previous been included in 3D simulations. Using one–zone models, Wolcott-
Green & Haiman (2011) found that this HI–shielding of ${\rm H_{2}}$ can
significantly decrease the LW photodissociation rate when the column density
exceeds $N_{\rm HI}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{22}~{}\rm{cm^{-3}}$ and provided
an analytic fit for the suppression factor $f_{\rm shield,HI}$.

In a study of the escape fraction of LW out of ACHs, Schauer et al. (2015)
used the WH11 fitting formula to quantify the effect of HI absorption of LW
photons emitted stars within the halo. They found the LW escape fraction was
reduced by a factor of three, owing to the large neutral column density. In a
later study of more massive halos, $10^{7-8}M_{\odot}$, Schauer et al. (2017)
found a significantly smaller effect, with escape fractions reduced by up to
$\sim 29$ per cent due to HI absorption, possibly due to significantly more
ionization by stellar clusters in the haloes resulting in lower neutral column
densities. Nevertheless, these results point to the possible importance of
HI–shielding of ${\rm H_{2}}$ in primordial ACHs irradiated by an external LW
field.

In this study, we use the three–dimensional hydrodynamic simulation code enzo
to test the effect of absorption by HI on $J_{\rm crit}$ in three
UV–irradiated protogalaxies collapsing at z$\sim 10$. We use a modified
version of the WH11 fitting formula for HI–shielding of ${\rm H_{2}}$ that we
updated to better fit non–ground state ${\rm H_{2}}$ rovibrational
populations, which become important at the temperatures and densities of
gravitationally collapsing ACHs. Our modified fitting formula, obtained using
data for the ${\rm H_{2}}$ rovibrational populations from cloudy, is accurate
to within $\sim 30$ per cent at $T=500-8000$K, $n=10^{0-5}~{}{\rm cm^{-3}}$,
and column densities $N_{\rm HI}\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}10^{24}~{}cm^{-2}$,
and can be easily implemented in chemical models for future studies.

The rest of this paper is organized as follows. In § 2 we provide details of
our simulations and calculations of the HI–shielding; we discuss our results
in § 3 and conclude with a summary in § 4.

##  2 Numerical Modeling

Figure 1: The effect of shielding by HI on the H2 photodissociation rate is
parameterized by a shielding factor $f_{\rm HI}=k_{\rm diss}(n,T,N_{\rm
H_{2}},N_{\rm HI})/k_{\rm diss}(n,T,N_{\rm H_{2}})$. The blue dashed lines
show the HI–shielding factor from the full calculation, $f_{\rm exact}$, with
cloudy–derived rovibrational populations for ${\rm H_{2}}$. The magenta solid
lines show our fit. The black dot-dashed line in each lower panel shows the
ratio of our fit to the exact calculation. In order to isolate the HI fitting
formula accuracy from that for self–shielding, we calculate $f_{\rm{fit}}$ as
= $f_{\rm{HI,fit}}\times f_{\rm{H_{2},exact}}$. All are shown for
$\log(n/{\rm{cm^{-3}}})=3$, and $\log(N_{\rm{H_{2}}}/{\rm{cm^{-2}}})=16$.
Figure 2: Spherically–averaged radial profiles for Halo B at the collapse
redshift showing density (upper left), temperature (upper right), electron
fraction (lower left), and ${\rm H_{2}}$ fraction (lower right). Lyman-Werner
fluxes are $J_{21}=20,000$ (sub–critical; blue solid lines) and
$J_{21}=32,000$ (super–critical; magenta dashed lines). The radial distance is
measured in physical units. Figure 3: Phase plots showing Halo B at the
collapse redshift with sub–critical flux (left) and super–critical flux
(right). Figure 4: Histogram of HI column densities along lines of sight from
the center of Halo B just before cooling occurs (sub–critical flux). Column
densities exceed the $N_{\rm{HI}}=10^{22}\ \rm{cm}^{-2}$ threshold, above
which HI–shielding of ${\rm H_{2}}$ becomes significant.

###  2.1 Simulations

We run simulations of three atomic cooling halos using enzo, a publicly-
available three-dimensional adaptive mesh refinement (AMR) hydrodynamic code
(Bryan et al., 2014). Initial conditions for a box $1h^{-1}$ Mpc on a side and
$128^{3}$ root grid were generated using music (Hahn & Abel, 2011). In order
to select haloes for higher resolution “zoom–in” simulations, we performed an
initial dark-matter only enzo run from $z=99$ to $z=10$. We used the rockstar
halo finder package (Behroozi et al., 2013) to identify a halo with mass above
the atomic cooling threshold at $z=10$; we then added three additional levels
of refinement using nested grids which enclose the Lagrangian volume for the
halo of interest, yielding an effective $1024^{3}$ resolution and dark matter
particle mass $\sim 85{\rm M_{\odot}}$.

Each halo is then run with $+$ DM “zoom-in” simulations initialized at $z=99$
to a maximum refinement level of 18, which results in the highest resolution
regions having a minimum cell size of $.0298h^{-1}$pc. Additional refinement
is added each time the baryon or dark matter mass exceeds four times that of
the most refined cell. We also imposed that the local Jeans length is resolved
by at least 16 cells to prevent spurious fragmentation (Truelove et al.,
1997).

We utilize the 9–species non–equilibrium primordial chemistry network within
enzo to model the chemical evolution of the gas. The cooling function from
Galli & Palla (1998) is implemented to model the radiative cooling by ${\rm
H_{2}}$. Several of the reaction rate calculations have been modified in the
enzo chemistry code as described in Wolcott-Green et al. 2021 (see their
Appendix A for details). For ${\rm H_{2}}$ self–shielding, we use the local
column density from the “Sobolev–like” method described in WGHB11 and their
fitting formula for the optically–thick rate.111This fit has since been
updated by Wolcott-Green & Haiman (2019) to account for non-ground state
rovibrational populations; however, using the updated fit would not affect our
conclusions, since we are interested here in the change in $J_{\rm crit}$ due
to HI–shielding, rather than the precise value of the critical flux. We assume
a blackbody incident radiation field at temperature $T=10^{5}$ K.

In order to determine the impact of HI–shielding, we run realizations of each
halo to determine the value of $J_{\rm crit}$ first with ${\rm H_{2}}$ self-
shielding only, using the Newton-Raphson method to find $J_{\rm crit}$, and
then subsequently run each with HI–shielding included, using our new fitting
formula § 2.2.

We use the publicly-available package YT (YT10) for simulation data analysis
and visualization222yt-project.org. Throughout, we adopt the cosmological
parameters from the Planck 2018 collaboration (Planck Collaboration et al.,
2018), $\Omega_{\rm m}=0.315$, $\Omega_{\Lambda}=0.685$, $\Omega_{b}=0.0493$,
$h=0.674$, $\sigma_{8}=0.811$, and $n=0.965$.

###  2.2 HI–shielding of ${\bf H_{2}}$

In order to find the exact optically–thick photodissociation rate, we use the
method described in detail in WGH19 and summarized briefly here. The rate
calculation includes contributions from LW transitions originating in the 301
bound rovibrational levels of the electronic ground state. We use the spectral
synthesis code cloudy (Ferland et al., 2017) to model the ${\rm H_{2}}$
rovibrational populations at $T=(500-8000)$K, densities $T=10^{(0-5)}~{}{\rm
cm^{-3}}$, $N_{\rm HI}=10^{(20-25)}~{}{\rm cm^{-}2}$, and $N_{\rm
H2}=10^{(14-17)}~{}{\rm cm^{-}2}$. The fractional populations are then input
in the rate calculation for each density and temperature combination.

In order to determine the impact of HI–shielding, the rate is calculated with
and without HI Lyman series absorption for each ${\rm n,T,N_{H2}}$. We define
the dimensionless HI shield factor as:

$f_{\rm sh,HI}=\frac{k_{\rm diss}(n{\rm,T,N_{HI},N_{H2}})}{k_{\rm diss}(n,{\rm
T,N_{H2}})}.$ (1)

In order to develop our fit, we began with the form used in WH11,

$f_{\rm sh,HI}=\frac{\chi}{(1+x)^{\delta}}\exp{(-\alpha x)}$ (2)

which was used in that study for photodissociation of ${\rm H_{2}}$ in the
ground rovibrational state only; we modified the parameters using the downhill
simplex method amoeba provided in Numerical Recipes.

Here, $x=N_{\rm HI}/\zeta$ and our best fit parameters are:

$\displaystyle\noindent\alpha=1.45\times 10^{-1}$ $\displaystyle\delta=1.5$
$\displaystyle\noindent\zeta=2.85\times 10^{23}\rm{cm^{-2}}$
$\displaystyle\chi=\left\\{\begin{tabular}[]{ll}$1$,&${\rm
N_{HI}<10^{22}~{}cm^{-2}}$\\\ $0.95$,&${\rm N_{HI}\geq 10^{22}~{}cm^{-2}}$\\\
\end{tabular}\right.$

Figure 1 shows the new HI-shielding factor fit, $f_{\rm fit}$ (Equation 2),
and the exact shielding factor from the full calculation, $f_{\rm exact}$, at
fixed $\log(n/{\rm{cm^{-3}}})=3$, $\log(N_{\rm{H_{2}}}/{\rm{cm^{-2}}})=16$ and
a range of temperatures. In the lower part of each panel is the ratio between
$f_{\rm fit}$ and $f_{\rm exact}$. The fitting formula for the shielding
function of ${\rm H_{2}}$ by HI is robust in the range of temperatures studied
here, $500-8000$K and is accurate to within a factor of two for column
densities $10^{20-24}~{}{\rm cm^{-2}}$.

##  3 Results

Table 1: Critical fluxes with and without HI–shielding in $J_{21}$ units. Virial masses and collapse redshifts indicated for $J<$$J_{\rm crit}$ runs with HI–shielding. Halo | $M/10^{7}M_{\odot}$ | $z_{\rm{coll}}$ | $T_{vir}/\rm{K}$ |  $J_{\rm crit}$$/10^{3}$ |  $J_{\rm crit}$${}_{\rm{,HI}}/10^{3}\ $
---|---|---|---|---|---
A | $2.8$ | $13.0$ | $8,296$ | $11$ | $22$
B | $8.2$ | $10.9$ | $14,418$ | $17$ | $32$
C | $2.4$ | $18.0$ | $10,231$ | $10$ | $16$

Figure 2 shows spherically–averaged radial profile of density, temperature,
electron fraction, and ${\rm H_{2}}$ fraction for Halo A at the collapse
redshift. Results for both sub–critical ($J_{21}=20,000$) and super–critical
($J_{21}=32,000$) LW fluxes are shown. Our halos follow the well–known
behavior of ACHs cooling in the presence of a photodissociating flux: with
$J<J_{\rm crit}$, the ${\rm H_{2}}$ fraction in the halo’s dense core reaches
$\sim 10^{-3}$, the standard “freeze–out” value (Oh & Haiman, 2002); ${\rm
H_{2}}$ cooling is efficient and the gas temperature falls below $10^{3}$K in
the dense core. Irradiation by a super–critical flux results in a suppressed
${\rm H_{2}}$ fraction, $\sim 10^{-7}$, and the temperature remains at $T_{\rm
vir}\sim 10^{4}\rm{K}$ throughout the halo.

To determine $J_{\rm crit}$, we run the zoom-in simulations with varied levels
of incident $J_{\rm LW}$ to find minimum flux that suppress ${\rm
H_{2}}$–cooling and prevents cooling below the virial temperature. The
resulting $J_{\rm crit}$ values for each of the three haloes are listed in
Table 1. We find that the critical flux with HI–shielding of ${\rm H_{2}}$ is
$\sim$60-100 percent larger than without HI–shielding. The $J_{\rm crit,21}$
values without HI–shielding for these halos are within the range
$(10-17)\times 10^{3}$, comparable to those found in previous studies, and
with HI–shielding $J_{\rm crit,21}$ are within the range $(16-32)\times
10^{3}$. Figure 3 shows phase plots for both sub–critical and super–critical
fluxes for Halo B at the collapse redshift.

The increase in $J_{\rm crit}$ with HI–shielding of ${\rm H_{2}}$indicates
that the neutral column densities are sufficient in these ACHs for Lyman
series absorption to be important, which occurs at ${\rm N_{HI}\mathrel{\hbox
to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10^{22}~{}cm^{-2}}$ (see Figure 1). In order to verify this,
we find the column densities at the critical density333above which,
collisional dissociation dominates the total ${\rm H_{2}}$ destruction rate by
summing along sightlines extending from the densest point of the halo out to a
radius of $100$ pc. We show a histogram of 25 such sightlines in Figure 4 at
the time when the gas in the core has reached at the critical density, just
before runaway cooling occurs. (Results are shown here for Halo B and those
from the other two halos are similar.We see that indeed the neutral columns
have reached the threshold at which HI–shielding becomes significant.

##  4 Conclusions

In this study, we examined the impact of HI–shielding of ${\rm H_{2}}$ on the
thermal evolution of protogalactic atomic cooling haloes exposed to
photodissociating UV radiation using three–dimensional hydrodynamic
simulations. We find that incorporation of HI–shielding raised the value of
the critical flux to suppress ${\rm H_{2}}$–radiative cooling, $J_{\rm crit}$,
by $\sim 60-100\%$ in the three ACHs we studied. This increase may have
important implications for the predicted number of candidate halos that could
seed massive black holes at $z\sim 10$ via direct collapse, which is sensitive
to the critical flux.

We used an updated fitting formula to model the suppression of ${\rm H_{2}}$
photodissociation by HI, which can be used in future simulations. The modified
fitting function is accurate to within $\sim 30$ per cent at $T=500-8000$K,
$n=10^{(0-5)}~{}{\rm cm^{-3}}$ and $N_{\rm HI}=10^{(20-24)}{\rm cm^{-2}}$.

## Acknowledgments

We thank Zoltán Haiman and S. Peng Oh for helpful discussions during the
course of this work. Meredith Neyer acknowledges funding from an Edison STEM
summer research program grant at University of California Santa Barbara. This
material is based upon work supported by the National Science Foundation under
Award No. 1903935. This work used the Extreme Science and Engineering
Discovery Environment (XSEDE; allocation TG-PHY200043), which is supported by
National Science Foundation grant number ACI-1548562.

##  5 Data availability

The data underlying this paper will be shared on reasonable request to the
corresponding author.

## References

  * Agarwal & Khochfar (2014) Agarwal B., Khochfar S., 2014, MNRAS, submitted, e-print ArXiv:1407.4115
  * Agarwal et al. (2012) Agarwal B., Khochfar S., Johnson J. L., Neistein E., Dalla Vecchia C., Livio M., 2012, MNRAS, 425, 2854
  * Ahn et al. (2009) Ahn K., Shapiro P. R., Iliev I. T., Mellema G., Pen U., 2009, ApJ, 695, 1430
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
  * Bromm & Yoshida (2011) Bromm V., Yoshida N., 2011, ARA&A, 49, 373
  * Bryan et al. (2014) Bryan G. L., Norman M. L., O’Shea B. W., Abel T., Wise J. H., Turk M. J., Reynolds D. R., Collins D. C., Wang P., Skillman S. W., 2014, ApJS, 211, 19
  * Chon et al. (2016) Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, ApJ, 832, 134
  * Dijkstra et al. (2014) Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS, 442, 2036
  * Dijkstra et al. (2008) Dijkstra M., Haiman Z., Mesinger A., Wyithe J. S. B., 2008, MNRAS, 391, 1961
  * Draine & Bertoldi (1996) Draine B. T., Bertoldi F., 1996, ApJ, 468, 269
  * Fan et al. (2001) Fan X., Narayanan V. K., Lupton R. H., Strauss M. A., Knapp G. R., Becker R. H., White R. L., Pentericci L., Leggett S. K., Haiman Z., Gunn J. E., Ivezić Ž., Schneider D. P., Anderson S. F., Brinkmann J., Bahcall N. A., Connolly A. J., Csabai I., Doi M., Fukugita M., Geballe T., Grebel E. K., Harbeck D., Hennessy G., Lamb D. Q., Miknaitis G., Munn J. A., Nichol R., Okamura S., Pier J. R., Prada F., Richards G. T., Szalay A., York D. G., 2001, AJ, 122, 2833
  * Fan et al. (2003) Fan X., Strauss M. A., Schneider D. P., Becker R. H., White R. L., Haiman Z., Gregg M., Pentericci L., Grebel E. K., Narayanan V. K., Loh Y.-S., Richards G. T., Gunn J. E., Lupton R. H., Knapp G. R., Ivezić Ž., Brandt W. N., Collinge M., Hao L., Harbeck D., Prada F., Schaye J., Strateva I., Zakamska N., Anderson S., Brinkmann J., Bahcall N. A., Lamb D. Q., Okamura S., Szalay A., York D. G., 2003, AJ, 125, 1649
  * Ferland et al. (2017) Ferland G. J., Chatzikos M., Guzmán F., Lykins M. L., van Hoof P. A. M., Williams R. J. R., Abel N. P., Badnell N. R., Keenan F. P., Porter R. L., Stancil P. C., 2017, Rev. Mex. Astron. Astrofis, 53, 385
  * Galli & Palla (1998) Galli D., Palla F., 1998, A&A, 335, 403
  * Haemmerlé et al. (2018) Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS, 474, 2757
  * Hahn & Abel (2011) Hahn O., Abel T., 2011, MNRAS, 415, 2101
  * Haiman et al. (2000) Haiman Z., Abel T., Rees M. J., 2000, ApJ, 534, 11
  * Haiman et al. (1997) Haiman Z., Rees M. J., Loeb A., 1997, ApJ, 476, 458
  * Hartwig et al. (2015) Hartwig T., Glover S. C. O., Klessen R. S., Latif M. A., Volonteri M., 2015, MNRAS, 452, 1233
  * Inayoshi et al. (2020) Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27
  * Kulkarni et al. (2020) Kulkarni M., Visbal E., Bryan G. L., 2020, arXiv e-prints, p. arXiv:2010.04169
  * Machacek et al. (2001) Machacek M. E., Bryan G. L., Abel T., 2001, ApJ, 548, 509
  * Mazzucchelli et al. (2017) Mazzucchelli C., Bañados E., Venemans B. P., Decarli R., Farina E. P., Walter F., Eilers A. C., Rix H. W., Simcoe R., Stern D., Fan X., Schlafly E., De Rosa G., Hennawi J., Chambers K. C., Greiner J., Burgett W., Draper P. W., Kaiser N., Kudritzki R. P., Magnier E., Metcalfe N., Waters C., Wainscoat R. J., 2017, ApJ, 849, 91
  * Mesinger et al. (2006) Mesinger A., Bryan G. L., Haiman Z., 2006, ApJ, 648, 835
  * Morganson et al. (2012) Morganson E., De Rosa G., Decarli R., Walter F., Chambers K., McGreer I., Fan X., Burgett W., Flewelling H., Greiner J., Hodapp K., Kaiser N., Magnier E., Price P., Rix H.-W., Sweeney B., Waters C., 2012, AJ, 143, 142
  * Oh & Haiman (2002) Oh S. P., Haiman Z., 2002, ApJ, 569, 558
  * Omukai (2001) Omukai K., 2001, ApJ, 546, 635
  * O’Shea & Norman (2008) O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14
  * Planck Collaboration et al. (2018) Planck Collaboration Aghanim N., Akrami Y., Ashdown M., Aumont J., Baccigalupi C., Ballardini M., Banday A. J., Barreiro R. B., Bartolo N., 2018, A&A, submitted, e-print arXiv:1807.06209
  * Regan et al. (2017) Regan J. A., Visbal E., Wise J. H., Haiman Z., Johansson P. H., Bryan G. L., 2017, Nature Astronomy, 1, 0075
  * Schauer et al. (2017) Schauer A. T. P., Agarwal B., Glover S. C. O., Klessen R. S., Latif M. A., Mas-Ribas L., Rydberg C.-E., Whalen D. J., Zackrisson E., 2017, MNRAS, 467, 2288
  * Schauer et al. (2015) Schauer A. T. P., Whalen D. J., Glover S. C. O., Klessen R. S., 2015, MNRAS, 454, 2441
  * Shang et al. (2010) Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249
  * Sugimura et al. (2014) Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS, 445, 544
  * Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman John H. I., Howell L. H., Greenough J. A., 1997, ApJL, 489, L179
  * Visbal et al. (2014) Visbal E., Haiman Z., Bryan G. L., 2014, MNRAS, 445, 1056
  * Wang et al. (2021) Wang F., Yang J., Fan X., Hennawi J. F., Barth A. J., Banados E., Bian F., Boutsia K., Connor T., Davies F. B., Decarli R., Eilers A.-C., Farina E. P., Green R., Jiang L., Li J.-T., Mazzucchelli C., Nanni R., Schindler J.-T., Venemans B., Walter F., Wu X.-B., Yue M., 2021, ApJL, 907, L1
  * Wang et al. (2019) Wang F., Yang J., Fan X., Wu X.-B., Yue M., Li J.-T., Bian F., Jiang L., Bañados E., Schindler J.-T., Findlay J. R., Davies F. B., Decarli R., Farina E. P., Green R., Hennawi J. F., Huang Y.-H., Mazzuccheli C., McGreer I. D., Venemans B., Walter F., Dye S., Lyke B. W., Myers A. D., Haze Nunez E., 2019, ApJ, 884, 30
  * Wise & Abel (2007) Wise J. H., Abel T., 2007, ApJ, 671, 1559
  * Wolcott-Green & Haiman (2011) Wolcott-Green J., Haiman Z., 2011, MNRAS, 412, 2603
  * Wolcott-Green & Haiman (2019) Wolcott-Green J., Haiman Z., 2019, MNRAS, 484, 2467
  * Wolcott-Green et al. (2011) Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
  * Wolcott-Green et al. (2017) Wolcott-Green J., Haiman Z., Bryan G. L., 2017, MNRAS, 469, 3329
  * Wolcott-Green et al. (2021) Wolcott-Green J., Haiman Z., Bryan G. L., 2021, MNRAS, 500, 138
  * Yang et al. (2019) Yang J., Wang F., Fan X., Yue M., Wu X.-B., Li J.-T., Bian F., Jiang L., Bañados E., Beletsky Y., 2019, AJ, 157, 236
  * Yoshida et al. (2003) Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, ApJ, 592, 645


# Efficient Penalized Generalized Linear Mixed Models for Variable Selection
and Genetic Risk Prediction in High-Dimensional Data

JULIEN ST-PIERRE
Department of Epidemiology, Biostatistics and Occupational Health,
McGill University, Montreal, Quebec, Canada
<EMAIL_ADDRESS>
KARIM OUALKACHA
Département de Mathématiques,
Université du Québec à Montréal, Montreal, Quebec, Canada
SAHIR RAI BHATNAGAR
Department of Epidemiology, Biostatistics and Occupational Health,
McGill University, Montreal, Quebec, Canada To whom correspondence should be
addressed.

###### Abstract

Sparse regularized regression methods are now widely used in genome-wide
association studies (GWAS) to address the multiple testing burden that limits
discovery of potentially important predictors. Linear mixed models (LMMs) have
become an attractive alternative to principal components (PC) adjustment to
account for population structure and relatedness in high-dimensional penalized
models. However, their use in binary trait GWAS rely on the invalid assumption
that the residual variance does not depend on the estimated regression
coefficients. Moreover, LMMs use a single spectral decomposition of the
covariance matrix of the responses, which is no longer possible in generalized
linear mixed models (GLMMs). We introduce a new method called pglmm, a
penalized GLMM that allows to simultaneously select genetic markers and
estimate their effects, accounting for between-individual correlations and
binary nature of the trait. We develop a computationally efficient algorithm
based on PQL estimation that allows to scale regularized mixed models on high-
dimensional binary trait GWAS ($\sim 300,000$ SNPs). We show through
simulations that penalized LMM and logistic regression with PC adjustment fail
to correctly select important predictors and/or that prediction accuracy
decreases for a binary response when the dimensionality of the relatedness
matrix is high compared to pglmm. Further, we demonstrate through the analysis
of two polygenic binary traits in the UK Biobank data that our method can
achieve higher predictive performance, while also selecting fewer predictors
than a sparse regularized logistic lasso with PC adjustment. Our method is
available as a Julia package PenalizedGLMM.jl.

##  1 Introduction

Genome-wide association studies (GWAS) have led to the identification of
hundreds of common genetic variants, or single nucleotide polymorphisms
(SNPs), associated with complex traits (Visscher et al.,, 2017) and are
typically conducted by testing association on each SNP independently. However,
these studies are plagued with the multiple testing burden that limits
discovery of potentially important predictors. Moreover, GWAS have brought to
light the problem of missing heritability, that is, identified variants only
explain a low fraction of the total observed variability for traits under
study (Manolio et al.,, 2009). Multivariable regression methods, on the other
hand, simultaneously fit many SNPs in a single model and are exempt from the
multiple testing burden. In both simulations and analysis of high-dimensional
data, sparse regularized logistic models have shown to achieve lower false-
positive rates and higher precision than methods based on univariable GWAS
summary statistics in case-control studies (Hoggart et al.,, 2008; Privé et
al.,, 2019). Contrary to univariable methods which implicitly assume that SNPs
are independent, a regularized model makes use of the linkage disequilibrium
(LD) structure between different loci, assigning weights to SNPs based on
their relative importance after accounting for all other SNPs already in the
model.

Confounding due to population structure or subject relatedness is another
major issue in genetic association studies. Modern large scale cohorts will
often include participants from different ethnic groups as well as admixed
individuals, that is, subjects with individual-specific proportions of
ancestries, or individuals with known or unknown familial relatedness, defined
as cryptic relatedness (Sul et al.,, 2018). Confounding comes from the fact
that allele frequencies can differ greatly between individuals who do not
share similar ancestry. When ignored, population structure and subject
relatedness can decrease power and lead to spurious associations (Astle and
Balding,, 2009; Price et al.,, 2010). Common practice is still to drop samples
by applying filters for relatedness or genetic ancestry, which can result in
decreasing the sample size by nearly 30% (Loh et al.,, 2018) in the full UK
Biobank data set (Bycroft et al.,, 2018).

Principal component analysis (PCA) can control for the confounding effect due
to population structure by including the top eigenvectors of a genetic
similarity matrix (GSM) as fixed effects in the regression model (Price et
al.,, 2006). With admixture and population structure being low dimensional
fixed-effects processes, they can correctly be accounted for by using a
relatively small number of PCs (e.g. 10) (Astle and Balding,, 2009; Novembre
and Stephens,, 2008). However, using too few PCs can result in residual bias
leading to false positives, while adding too many PCs as covariates can lead
to a loss of efficiency (Zhao et al.,, 2018). Alternatively, using mixed
models (MMs), one can model population structure and/or closer relatedness by
including a polygenic random effect with variance-covariance structure
proportional to the GSM (Yu et al.,, 2005). Indeed, kinship is a high-
dimensional process, such that it cannot be fully captured by a few PCs
(Hoffman,, 2013). Thus, it would require the inclusion of too many PCs as
covariates, relative to the dimension of the sample size. Hence, while both
PCA and MMs share the same underlying model, MMs are more robust in the sense
that they do not require distinguishing between the different types of
confounders (Price et al.,, 2010). Moreover, MMs alleviate the need to
evaluate the optimal number of PCs to retain in the model as fixed effects.

Several authors have proposed to combine penalized quasi-likelihood (PQL)
estimation with sparsity inducing regularization to perform selection of fixed
and/or random effects in generalized linear mixed model (GLMMs) (Groll and
Tutz,, 2014; Hui et al.,, 2017). However, none of these methods are currently
scalable for modern large-scale genome-wide data, nor can they directly
incorporate relatedness structure through the use of a kinship matrix. Indeed,
the computational efficiency of recent multivariable methods for high-
dimensional MMs rely on performing a spectral decomposition of the covariance
matrix to rotate the phenotype and design matrix such that the transformed
data become uncorrelated (Bhatnagar et al., 2020b, ; Rakitsch et al.,, 2012).
These methods are typically restricted to linear models since in GLMMs, it is
no longer possible to perform a single spectral decomposition to rotate the
phenotype and design matrix, as the covariance matrix depends on the sample
weights which in turn depend on the estimated regression coefficients that are
being iteratively updated. This limits the application of high-dimensional MMs
to analysis of binary traits in genetic association studies.

In this paper, we introduce a new method called pglmm that allows to
simultaneously select variables and estimate their effects, accounting for
between-individual correlations and binary nature of the trait. We develop a
scalable algorithm based on PQL estimation which makes it possible, for the
first time, to fit penalized GLMMs on high-dimensional GWAS data. To speedup
computation, we estimate the variance components and dispersion parameter of
the model under the null hypothesis of no genetic effect. Secondly, we use an
upper-bound for the inverse variance-covariance matrix in order to perform a
single spectral decomposition of the GSM and greatly reduce memory usage.
Finally, we implement an efficient block coordinate descent algorithm in order
to find the optimal estimates for the fixed and random effects parameters. Our
method is implemented in an open source Julia programming language (Bezanson
et al.,, 2017) package called PenalizedGLMM.jl and freely available at
https://github.com/julstpierre/PenalizedGLMM.

The rest of this paper is structured as follows. In Section 2 we present our
model and describe the block coordinate gradient descent algorithm that is
used to estimate the model parameters. We also discuss several approaches to
select the optimal tuning parameter in regularized models, and we detail how
predictions are obtained in GLMs with PC adjustment versus our proposed mixed
model. In Section 3, we show through simulations that both LMM and logistic
model with PC adjustment fail to correctly select important predictors and
estimate their effects when the dimensionality of the kinship matrix is high.
Further, we demonstrate through the analysis of two polygenic binary traits in
the UKBB data that our method achieves higher predictive performance, while
also selecting consistently fewer predictors than a logistic lasso with PC
adjustment. We finish with a discussion of our results, some limitations and
future directions in Section 4.

##  2 Methods

###  2.1 Model

We consider the following GLMM

$\displaystyle
g(\mu_{i})=\eta_{i}=\bm{X}i\bm{\alpha}+\bm{G}i\bm{\gamma}+b_{i},$ (1)

for $i=1,..,n$, where
$\mu_{i}=\mathbb{E}(y_{i}=1|\bm{X}_{i},\bm{G}_{i},b_{i})$, $\bm{X}_{i}$ is a
$1\times m$ row vector of covariates for subject $i$, $\bm{\alpha}$ is a
$m\times 1$ column vector of fixed covariate effects including the intercept,
$\bm{G}_{i}$ is a $1\times p$ row vector of genotypes for subject $i$ taking
values $\\{0,1,2\\}$ as the number of copies of the minor allele, and
$\bm{\gamma}$ is a $p\times 1$ column vector of fixed additive genotype
effects. We assume that
$\bm{b}=(b_{1},...,b_{n})^{\intercal}\sim\mathcal{N}(0,\sum_{s=1}^{S}\tau_{s}\bm{V}_{s})$
is an $n\times 1$ column vector of random effects,
$\bm{\tau}=(\tau_{1},...,\tau_{S})^{\intercal}$ are variance component
parameters, and $\bm{V}_{s}$ are known $n\times n$ relatedness matrices. The
phenotypes $y_{i}$ are assumed to be conditionally independent and identically
distributed given $(\bm{X}_{i},\bm{G}_{i},\bm{b})$ and follow any exponential
family distribution with canonical link function $g(\cdot)$, mean
$\mathbb{E}(y_{i}|\bm{b})=\mu_{i}$ and variance $\text{Var}(y_{i}|\bm{b})=\phi
a_{i}^{-1}\nu(\mu_{i}),$ where $\phi$ is a dispersion parameter, $a_{i}$ are
known weights and $\nu(\cdot)$ is the variance function. In order to estimate
the parameters of interest and perform variable selection, we need to use an
approximation method to obtain a closed analytical form for the marginal
likelihood of model (1). Following the derivation of Chen et al., (2016), we
propose to fit (1) using a penalized quasi-likelihood (PQL) method, from where
the log integrated quasi-likelihood function is equal to

$\displaystyle
ql(\bm{\alpha},\bm{\gamma},\phi,\bm{\tau})=-\frac{1}{2}\text{log}\left|\sum_{s=1}^{S}\tau_{s}\bm{V}_{s}\bm{W}+\bm{I}_{n}\right|+\sum_{i=1}^{n}ql_{i}(\bm{\alpha},\bm{\gamma}|\bm{\tilde{b}})-\frac{1}{2}\bm{\tilde{b}}^{\intercal}\left(\sum_{s=1}^{S}\tau_{s}\bm{V}_{s}\right)^{-1}\bm{\tilde{b}},$
(2)

where
$\bm{W}=\textrm{diag}\left\\{\frac{a_{i}}{\phi\nu(\mu_{i})[g^{\prime}(\mu_{i})^{2}]}\right\\}$
is a diagonal matrix containing weights for each observation,
$ql_{i}(\bm{\alpha,\gamma}|\bm{b})=\int_{y_{i}}^{\mu_{i}}\frac{a_{i}(y_{i}-\mu)}{\phi\nu(\mu)}d\mu$
is the quasi-likelihood for the $ith$ individual given the random effects
$\bm{b}$, and $\tilde{\bm{b}}$ is the solution which maximizes (2).

In typical genome-wide studies, the number of predictors is much greater than
the number of observations ($p>n$), and the parameter vector $\bm{\gamma}$
becomes underdetermined when modelling all SNPs jointly. Thus, we propose to
add a lasso regularization term (Tibshirani,, 1996) to the negative quasi-
likelihood function in (2) to seek a sparse subset of $\bm{\gamma}$ that gives
an adequate fit to the data. Because
$ql(\bm{\alpha},\bm{\gamma},\phi,\bm{\tau})$ is a non-convex loss function, we
propose a two-step estimation method to reduce the computational complexity.
First, we obtain the variance component estimates $\hat{\phi}$ and
$\bm{\hat{\tau}}$ under the null hypothesis of no genetic effect
($\bm{\gamma}=\bm{0}$) using the AI-REML algorithm (Gilmour et al.,, 1995;
Chen et al.,, 2016) detailed in Appendix A of the supplementary material.
Assuming that the weights in $\bm{W}$ vary slowly with the conditional mean,
we drop the first term in (2) (Breslow and Clayton,, 1993) and define the
following objective function which we seek to minimize with respect to
$(\bm{\alpha},\bm{\gamma},\tilde{\bm{b}})$:

$\displaystyle(\hat{\bm{\alpha}},\hat{\bm{\gamma}},\hat{\bm{b}})$
$\displaystyle=\underset{\bm{\alpha},\bm{\gamma},\tilde{\bm{b}}}{\text{argmin
}}Q_{\lambda}(\bm{\alpha},\bm{\gamma},\tilde{\bm{b}}),$ $\displaystyle
Q_{\lambda}(\bm{\alpha},\bm{\gamma},\tilde{\bm{b}})$
$\displaystyle=-\sum_{i=1}^{n}ql_{i}(\bm{\alpha},\bm{\gamma}|\bm{\tilde{b}})+\frac{1}{2}\bm{\tilde{b}}^{\intercal}\left(\sum_{s=1}^{S}\hat{\tau}_{s}\bm{V}_{s}\right)^{-1}\bm{\tilde{b}}+\lambda\sum_{j}v_{j}|\gamma_{j}|$
$\displaystyle:=-\ell_{PQL}(\bm{\alpha},\bm{\gamma},\hat{\phi},\hat{\bm{\tau}})+\lambda\sum_{j}v_{j}|\gamma_{j}|,$
(3)

where $\lambda$ is a nonnegative regularization parameter, and $v_{j}$ is a
penalty factor for the $j^{th}$ predictor.

In Appendix B, we detail our proposed general purpose block coordinate
gradient descent algorithm (CGD) to solve (2.1) and obtain regularized PQL
estimates for $\bm{\alpha},\bm{\gamma}$ and $\tilde{\bm{b}}$. Briefly, our
algorithm is equivalent to iteratively solve the two penalized generalized
least squares (GLS)

$\underset{\tilde{\bm{b}}}{\textrm{argmin}}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}{\bm{\beta}}-\tilde{\bm{b}}\right)^{\intercal}\bm{W}^{-1}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}{\bm{\beta}}-\tilde{\bm{b}}\right)+\tilde{\bm{b}}^{\intercal}\left(\sum_{s=1}^{S}\hat{\tau}_{s}\bm{V}_{s}\right)^{-1}\tilde{\bm{b}},$

and

$\underset{\bm{\beta}}{\textrm{argmin}}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}\right)^{\intercal}\bm{\Sigma}^{-1}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}\right)+\lambda\sum_{j}v_{j}|\beta_{j}|,$

where $\bm{\Sigma}=\bm{W}^{-1}+\sum_{s=1}^{S}\hat{\tau}_{s}\bm{V}_{s}$ is the
covariance matrix of the working response vector $\tilde{\bm{Y}}$,
$\tilde{\bm{X}}=\left[\bm{X};\ \bm{G}\right]$ and
$\bm{\beta}=(\bm{\alpha}^{\intercal},\bm{\gamma}^{\intercal})^{\intercal}$. We
use the spectral decomposition of $\bm{\Sigma}$ to rotate $\tilde{\bm{Y}}$,
$\tilde{\bm{X}}$ and $\tilde{\bm{b}}$ such that the transformed data become
uncorrelated. For binary data, because the covariance matrix $\bm{\Sigma}$
depends on the sample weights $\bm{W}$, we use an upper-bound on
$\bm{\Sigma}^{-1}$ to ensure a single spectral decomposition is performed
(Böhning and Lindsay,, 1988). By cycling through the coordinates and
minimizing the objective function with respect to one parameter at a time,
$\tilde{\bm{b}}$ can be estimated by fitting a generalized ridge-like model
with a diagonal penalty matrix equal to the inverse of the eigenvalues of
$\sum_{s=1}^{S}\hat{\tau}_{s}\bm{V}_{s}$. Then, conditional on
$\tilde{\bm{b}}$, ${\bm{\beta}}$ is estimated by solving a weighed least
squares (WLS) with a lasso regularization term. All calculations and
algorithmic steps are detailed in Appendix B.

###  2.2 Model selection

Approaches to selecting the optimal tuning parameter in regularized models are
of primary interest since in real data analysis, the underlying true model is
unknown. A popular strategy is to select the value of the tuning parameter
that minimizes out-of-sample prediction error, e.g., cross-validation (CV),
which is asymptotically equivalent to the Akaike information criterion (AIC)
(Akaike,, 1998; Yang,, 2005). While being conceptually attractive, CV becomes
computationally expensive for very high-dimensional data. Moreover, in studies
where the proportion of related subjects is important, either by known or
cryptic relatedness, the CV prediction error is no longer an unbiased
estimator of the generalization error (Rabinowicz and Rosset,, 2020). Through
simulation studies and real data analysis, Wang et al., (2020) found that LD
and minor allele frequencies (MAF) differences between ancestries could
explain between 70 and 80% of the loss of relative accuracy of European-based
prediction models in African ancestry for traits like body mass index and type
2 diabetes. Thus, there is no clear approach to how multiple admixed and/or
similar populations should be split when using CV to minimize out-of-sample
prediction error.

Alternatively, we can use the generalized information criterion (GIC) to
choose the optimal value of the tuning parameter $\lambda$, defined as

$\displaystyle\textrm{GIC}_{\lambda}=-2\ell_{PQL}+a_{n}\cdot\hat{df}_{\lambda},$
(4)

where $\ell_{PQL}$ is defined in (2.1), and $\hat{df}_{\lambda}=|\\{1\leq
k\leq p:\hat{\beta}_{k}\neq 0\\}|+\textrm{dim}(\hat{\bm{\tau}})$ is the number
of nonzero fixed-effects coefficients (Zou et al.,, 2007) plus the number of
variance components. Special cases of the GIC include AIC ($a_{n}=2$) and the
Bayesian information criterion (BIC) (Schwarz,, 1978) ($a_{n}=\text{log}(n)$).

###  2.3 Prediction

It is often of interest in genetic association studies to make predictions on
a new set of individuals, e.g., the genetic risk of developing a disease for a
binary response or the expected outcome in the case of a continuous response.
In what follows, we compare how predictions are obtained in pglmm versus a GLM
with PC adjustment.

####  2.3.1 pglmm

For the sake of comparison with the GLM with PC adjustment, we suppose a
sampling design where a single variance component is needed such that
$\bm{b}\sim\mathcal{N}(\bm{0},\tau_{1}\bm{V_{1}})$ where $\bm{V_{1}}$ is the
GSM between $n$ subjects that are used to fit the GLMM (1). We iteratively fit
on a training set the working linear mixed model

$\tilde{\bm{Y}}=\tilde{\bm{X}}\bm{\beta}+\bm{\bm{b}}+\bm{\epsilon},$

where
$\bm{\epsilon}=g^{\prime}(\bm{\mu})(\bm{y}-\bm{\mu})\sim\mathcal{N}(0,\bm{W}^{-1})$.
Let $\tilde{\bm{Y}}_{s}$ be the latent working vector in a set of individuals
with predictor set $\tilde{\bm{X}}_{s}$ that were not used in the model
training, $n_{s}$ denote the number of observations in the testing set and $n$
the number of observations in the training set. Similar to (Bhatnagar et al.,
2020b, ), we assume that the marginal joint distribution of
$\tilde{\bm{Y}}_{s}$ and $\tilde{\bm{Y}}$ is multivariate Normal :

$\displaystyle\begin{bmatrix}\tilde{\bm{Y}}_{s}\\\
\tilde{\bm{Y}}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\tilde{\bm{X}}_{s}\bm{\beta}\\\
\tilde{\bm{X}}\bm{\beta}\end{bmatrix},\begin{bmatrix}\bm{\Sigma}_{11}&\bm{\Sigma}_{12}\\\
\bm{\Sigma}_{21}&\bm{\Sigma}_{22}\end{bmatrix}\right),$

where $\bm{\Sigma}_{12}=\tau_{1}\bm{V}_{12}$ and $\bm{V}_{12}$ is the
$n_{s}\times n$ GSM between the testing and training individuals. It follows
from standard normal theory that

$\displaystyle\tilde{\bm{Y}}_{s}|\tilde{\bm{Y}},\phi,\tau_{1},\bm{\beta},\tilde{\bm{X}},\tilde{\bm{X}}_{s}\sim\mathcal{N}\left(\tilde{\bm{X}}_{s}\bm{\beta}+\bm{\Sigma}_{12}\bm{\Sigma}_{22}^{-1}(\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}),\bm{\Sigma}_{11}-\bm{\Sigma}_{12}\bm{\Sigma}_{22}^{-1}\bm{\Sigma}_{21}\right).$

The estimated mean response $\hat{{\bm{\mu}}}_{s}$ for the testing set is
given by

$\displaystyle
g^{-1}\left(\mathbb{E}[\tilde{\bm{Y}}_{s}|\tilde{\bm{Y}},\hat{\phi},\hat{\tau}_{1},\hat{\bm{\beta}},\tilde{\bm{X}},\tilde{\bm{X}}_{s}]\right)$
$\displaystyle=g^{-1}\left(\tilde{\bm{X}}_{s}\hat{\bm{\beta}}+\bm{\Sigma}_{12}\bm{\Sigma}_{22}^{-1}(\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}})\right)$
$\displaystyle=g^{-1}\left(\tilde{\bm{X}}_{s}\hat{\bm{\beta}}+\hat{\tau}_{1}\bm{V}_{12}\left(\bm{W}^{-1}+\hat{\tau}_{1}\bm{V}_{1}\right)^{-1}(\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}})\right)$
$\displaystyle=g^{-1}\left(\tilde{\bm{X}}_{s}\hat{\bm{\beta}}+\bm{V}_{12}\bm{U}\left(\frac{1}{\hat{\tau}_{1}}\bm{D}+\tilde{\bm{U}}^{\intercal}\bm{W}\tilde{\bm{U}}\right)^{-1}\tilde{\bm{U}}^{\intercal}\bm{W}(\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}})\right),$
(5)

where $g(\cdot)$ is a link function and $\tilde{\bm{U}}=\bm{UD}$ is the
$n\times n$ matrix of PCs obtained from the spectral decomposition of the GSM
for training subjects.

####  2.3.2 GLM with PC adjustment

Another approach to control for population structure and/or subjects
relatedness is to use the first $r$ columns of $\tilde{\bm{U}}$ as unpenalized
fixed effects covariates. This leads to the following GLM

$\displaystyle
g(\bm{\mu})=\tilde{\bm{X}}\bm{\beta}+\tilde{\bm{U}}_{r}\bm{\delta},$

where $\tilde{\bm{X}}=[\bm{X};\bm{G}]$ is the $n\times(m+p)$ design matrix for
non-genetic and genetic predictors, $\bm{\beta}\in\mathbb{R}^{p}$ is the
corresponding sparse vector of fixed effects, $\tilde{\bm{U}}_{r}$ is the
$n\times r$ design matrix for the first $r$ PCs and $\delta\in\mathbb{R}^{r}$
is the corresponding vector of fixed effects. Letting
$\tilde{\bm{Y}}=\tilde{\bm{X}}\bm{\beta}+\tilde{\bm{U}}_{r}\bm{\delta}+g^{\prime}(\bm{\mu})(\bm{y-\mu})$
be the working response vector, one can show that

$\displaystyle\hat{\bm{\delta}}=\left(\tilde{\bm{U}}_{r}^{\intercal}\bm{W}\tilde{\bm{U}}_{r}\right)^{-1}\tilde{\bm{U}}_{r}^{\intercal}\bm{W}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}}\right),$
(6)

where $\bm{W}$ is the diagonal matrix of GLM weights. Let $\bm{V}_{12}$ be the
$n_{s}\times n$ GSM between a testing set of $n_{s}$ individuals and $n$
training individuals such that the projected PCs on the testing subjects are
equal to $\bm{V}_{12}\bm{U}_{r}$. Then, the estimated mean response
$\hat{{\bm{\mu}}}_{s}$ for the testing set is given by

$\displaystyle\hat{\bm{\mu}}_{s}=g^{-1}\left(\tilde{\bm{X}_{s}}\hat{\bm{\beta}}+\bm{V}_{12}\bm{U}_{r}\hat{\bm{\delta}}\right)=g^{-1}\left(\tilde{\bm{X}_{s}}\hat{\bm{\beta}}+\bm{V}_{12}\bm{U}_{r}\left(\tilde{\bm{U}}_{r}^{\intercal}\bm{W}\tilde{\bm{U}}_{r}\right)^{-1}\tilde{\bm{U}}_{r}^{\intercal}\bm{W}\left(\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}}\right)\right).$
(7)

By comparing (2.3.1) and (7), we see that both GLM with PC adjustment and
pglmm use a projection of the training PCs on the testing set to predict new
responses, but with different coefficients for the projected PCs. For the
former, the estimated coefficients for the first $r$ projected PCs in (6) are
obtained by iteratively solving generalized least squares (GLS) on the partial
working residuals $\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}}$. For pglmm,
the estimated coefficients for all projected PCs are also obtained by
iteratively solving GLS on the partial working residuals
$\tilde{\bm{Y}}-\tilde{\bm{X}}\hat{\bm{\beta}}$, with an extra ridge penalty
for each coefficient that is equal to $\hat{\tau_{1}}^{-1}\Lambda_{i}$ with
$\Lambda_{i}$ the $i^{th}$ eigenvalue of $\bm{V}$ that is associated with the
$i^{th}$ PC.

From a Bayesian point of view, the fixed effect GLM assumes that each of the
$r$ selected PCs have equal prior probability, while the remaining $n-r$
components are given a prior with unit mass at zero (Astle and Balding,,
2009). In contrast, the MM puts a Gaussian prior on regression coefficients
with variances proportional to the corresponding eigenvalues. This implies
that the PCs with the largest eigenvalues have a higher prior probability of
explaining the phenotype. Hence, pglmm shrinks PCs coefficients in a smooth
way, while the fixed effect GLM uses a tresholding approach; the first $r$
predictors with larger eigenvalues are kept intact, and the others are
completely removed. This implies that the confounding effect from population
structure and/or relatedness on the phenotype is fully captured by the first
$r$ PCs. As we show in simulations results, departure from this assumption
leads to less accurate coefficients estimation, lower prediction accuracy and
higher variance in predictions.

###  2.4 Simulation design

We evaluated the performance of our proposed method against that of a lasso
LMM, using the R package ggmix (Bhatnagar et al., 2020a, ), and a logistic
lasso, using the Julia package GLMNet which wraps the Fortran code from the
original R package glmnet (Friedman et al.,, 2010). We compared glmnet when we
included or not the first 10 PCs in the model (glmnetPC). We performed a total
of 50 replications for each simulation scenario, drawing anew genotypes and
simulated traits. Values for all simulation parameters are presented in Table
1.

####  2.4.1 Simulated genotype from the admixture model

In the first scenario, we studied the performance of all methods for different
population structures by simulating random genotypes from the BN-PSD admixture
model for 10 or 20 subpopulations with 1D geography or independent
subpopulations using the bnpsd package in R (Ochoa and Storey, 2016a, ; Ochoa
and Storey, 2016b, ). Sample size was set to $n=2500$. We simulated $p$
candidate SNPs and randomly selected $p\times c$ to be causal. The kinship
matrix $\bm{V}$ and PCs were calculated using a set of $50,000$ additional
simulated SNPs. We simulated covariates for age and sex using Normal and
Binomial distributions, respectively.

####  2.4.2 Real genotypes from the UK Biobank data

In the second scenario, we compared the performance of all methods when a high
proportion of related individuals are present, using real genotype data from
the UK Biobank. We retained a total of 6731 subjects of White British ancestry
having estimated 1st, 2nd or 3rd degree relationships with at least one other
individual. We sampled $p$ candidate SNPS among all chromosomes and randomly
selected $p\times c$ to be causal. We used PCs as provided with the data set.
These were computed using a set of unrelated samples and high quality markers
pruned to minimise LD (Bycroft et al.,, 2018). Then, all subjects were
projected onto the principal components using the corresponding loadings.
Since the markers that were used to compute the PCs were potentially sampled
as candidate causal markers in our simulations, we included all candidate SNPs
in the set of markers used for calculating the kinship matrix $\bm{V}$. We
simulated age using a Normal distribution and used the sex covariate as
provided with the data.

####  2.4.3 Simulation model

The number of candidate SNPs and fraction of causal SNPs were set to $p=5000$
and $c=0.01$ respectively. Let $S$ be the set of candidate causal SNPs, with
$|S|=p\times c$, then the causal SNPs fixed effects $\beta_{j}$ were generated
from a Gaussian distribution $\mathcal{N}(0,h^{2}_{g}\sigma^{2}/|S|)$, where
$h^{2}_{g}$ is the fraction of variance on the logit scale that is due to
total additive genetic fixed effects. That is, we assumed the candidate causal
markers explained a fraction of the total polygenic heritability, and the rest
was explained by a random polygenic effect
$b\sim\mathcal{N}(0,h^{2}_{b}\sigma^{2}\bm{V})$. We simulated a SNR equal to 1
for the fixed genetic effects ($h^{2}_{g}=50\%$) under strong random polygenic
effects ($h^{2}_{b}=40\%$). Finally, we simulated a binary phenotype using a
logistic link function

$\displaystyle\text{logit}(\pi)=\text{logit}(\pi_{0})-\text{log}(1.3)\times
Sex+\text{log}(1.05)Age/10+\sum_{j\in S}\beta_{j}\cdot\widetilde{G}_{j}+b,$
(8)

where the parameter $\pi_{0}$ was chosen to specify the prevalence under the
null, and $\widetilde{G}_{j}$ is the $j^{th}$ column of the standardized
genotype matrix $\tilde{g}_{ij}=(g_{ij}-2p_{i})/\sqrt{2p_{i}(1-p_{i})}$ and
$p_{i}$ is the MAF.

####  2.4.4 Metric of comparison

For each replication, subjects were partitioned into training and test sets
using an 80/20 ratio, ensuring all related individuals were assigned into the
same set in the second scenario. Variable selection and coefficient estimation
were performed on training subjects for all methods. We compared each method
at a fixed number of predictors, ranging from 5 to 50 which corresponds to the
number of true causal SNPs. Comparisons were based on three criteria: the
ability to retrieve the causal predictors, measured by the true positive rate
${\textrm{TPR}=|\\{1\leq k\leq p:\hat{\beta}_{k}\neq 0\cap{\beta}_{k}\neq
0\\}|/|\\{1\leq k\leq p:\hat{\beta}_{k}\neq 0\\}|}$; the ability to accurately
estimate coefficients, measured by the root mean squared error
${\textrm{RMSE}=\sqrt{\frac{1}{p}\sum_{k=1}^{p}(\hat{\beta}_{k}-\beta_{k})^{2}}}$;
and the ability to predict outcomes in the test sets, measured by the area
under the roc curve (AUC).

In addition, we evaluated the performance of our proposed method when using
either AIC, BIC or cross-validation as model selection criteria, rather than
fixing the number of predictors in the model. For this, subjects from the real
UKBB data were randomly split into training (40%), validation (30%) and test
(30%) sets, again ensuring all related individuals were assigned into the same
set. For cross-validation, the full lasso solution path was fitted on the
training set, and the regularization parameter was obtained on the model which
maximized AUC on the validation set. We compared methods performance on the
basis of TPR, AUC on the test sets and RMSE. Additionally, we compared each
model selection approach on the total number of predictors selected and on the
model precision, which is defined as the proportion of selected predictors
that are true positives.

###  2.5 Real data application

We used the real UK Biobank data set presented in Section 2.4 to illustrate
the potential advantages of pglmm over logistic lasso with PC adjustment for
constructing a PRS on two highly heritable binary traits, asthma and high
cholesterol, in a set of related individuals. Asthma is a common respiratory
disease characterized by inflammation and partial obstruction of the bronchi
in the lungs that results in difficulty breathing (Anderson,, 2008). High
cholesterol can form plaques and fatty deposits on the walls of the arteries,
and thus prevent the blood to circulate to the heart and brain. It is one of
the main controllable risk factors for coronary artery disease, heart attack
and stroke (Kathiresan et al.,, 2008).

After filtering for SNPs with missing rate greater than $0.01$, MAF above
$0.05$ and a p-value for the Hardy–Weinberg exact test above $10^{-6}$, a
total of 320K genotyped SNPs were remaining. To better understand the
contribution of the PRS for predicting asthma and high cholesterol, we fitted
for each trait a null model with only age, sex, genotyping array and the first
10 PCs as main effects. For our proposed pglmm method, we did not include any
PC since kinship is accounted for by a random effect. Finally, we compared
with a logistic lasso in which the top 10 PCs were included as unpenalized
covariates in addition to age, sex and genotyping array. To find the optimal
regularization parameter for both methods, we split the subjects in training
(60%), validation (20%) and test (20%) sets for a total of 40 times. For each
replication, the full lasso solution path was fitted on the training set, and
the regularization parameter was obtained on the model which maximized AUC on
the validation set. We compared mean prediction accuracy on the test sets as
well as the median number of predictors included in all models. Finally, we
also compared our method’s performance when the best model was chosen using
BIC on the training fit.

##  3 Results

###  3.1 Simulation results from the admixture model

Results for selection of important predictors in the first simulation
scenario, as measured by the mean TPR in 50 replications, are presented in
Figure 1. For both 1D linear admixture and independent subpopulations, glmnet
without PC adjustment failed to retrieve causal markers compared to all other
methods. This is expected under population stratification; SNPs that differ in
frequency between subpopulations are identified as important predictors
because prevalence is not constant across each group. When the first 10 PCs
were added as unpenalized covariates, glmnetPC’s ability to select causal
predictors was lesser to that of pglmm and ggmix for the 20 independent
subpopulations. In this case, the rank of the GSM used to infer the PCs is
close to 20, and including only 10 PCs in the model does not correctly capture
the confounding structure. Because there is less overlap between
subpopulations in the admixture data compared to the independent populations
(Reisetter and Breheny,, 2021), a greater proportion of the simulated
polygenic random effect is explained by the GSM and including only 10 PCs is
enough to correct for confounding even when $K=20$ (bottom-left panel of
Figure 1). On the other hand, including a random effect with variance-
covariance structure proportional to the GSM correctly adjusts for population
structure in all scenarios while alleviating the burden of choosing the right
number of fixed predictors to include in the model. Even though ggmix assumes
a standard LMM for the binary trait, it was able to identify causal markers at
the same rate as pglmm.

Results for estimation of SNP effects as measured by the mean RMSE in 50
replications are presented in Figure 2. Results are consistent with TPR
results in that glmnet without PC adjustment performed poorly in all
scenarios, while pglmm outperformed all other methods for the 20 independent
subpopulations and performed comparably with glmnetPC for all other settings.
As expected, ggmix had higher RMSE compared to pglmm and glmnetPC. Thus, even
though ggmix was able to identify causal markers at the same rate as other
methods that accounted for the binary nature of the response, resulting
estimates for the SNP effects were not accurate.

For both 1D linear admixture and independent subpopulations, ggmix and glmnet
had poor predictive performance for $K=10$ and $K=20$, as reported in Figure
3. Also, the predictive performance of glmnetPC was greatly reduced when
$K=20$ for both admixture and independent populations. In the case of the
admixture data, the RMSE for estimation of SNP effects was comparable for
glmnetPC and pglmm. This means that the observed discrepancy in predictive
accuracy is due to the difference in how each method handle the confounding
effects. Using only 10 PCs as fixed effects when $K=20$ likely results in
overfitted coefficients, which may potentially decrease prediction accuracy
and increase variance of predictions in independent subjects. By using a
ridge-like estimator for the random effects, pglmm is less likely to overfit
the confounding effects compared to glmnetPC. This is supported by the results
of Table 2, where the relative decrease in AUC standard deviation for the
predictions obtained by pglmm could be as high as $16\%$ for $K=20$
subpopulations.

###  3.2 Simulation results from real genotype data

Results for selection of important predictors, estimation of SNPs effects and
prediction accuracy in the second simulation scenario are presented in Figure
4. We compared the ability of glmnetPC and pglmm to adjust for potential
confounding stemming from subjects relatedness. Both methods’ ability to
retrieve important predictors were comparable as measured by mean TPR, with
pglmm having a slight advantage. In terms of predictor effect estimation,
pglmm had lower reported mean RMSE. Furthermore, pglmm outperformed glmnetPC
when making predictions in independent test sets. Once again, this is
explained by the fact that pglmm uses a random effect parameterized by the
$n-$dimensional kinship matrix, which we have shown in Section 2.3 to be
equivalent to include all PCs as predictors in the model and shrink their
coefficients proportionally to their relative importance in a smooth way. On
the other hand, for glmnetPC, only the first $10$ PCs with larger eigenvalues
are kept intact, and the others are completely removed. As the confounding
effect from relatedness on the phenotype can not be fully captured by using
only the first $10$ PCs, prediction accuracy is greatly reduced.

###  3.3 Model selection

Boxplots of the model selection simulations results are presented in Figure 5.
As expected, BIC tended to choose sparser models with very high precision
values, compared to AIC and CV which tended to select larger models with
negligibly higher prediction performance. Thus, using BIC as a model selection
criteria resulted in trading a bit of prediction accuracy for a large boost in
the model precision. In many situations where it is of interest to identify a
smaller set containing the most important predictors, BIC should be preferred
over AIC and CV. Moreover, BIC alleviates the computational challenge of
performing out-of-sample predictions, which includes identifying pedigrees to
ensure independence between training, validation and testing sets.

###  3.4 PRS for the UK Biobank

Results for asthma and high cholesterol PRSs are summarized in Table 3. For
asthma, pglmm with either BIC or CV as model selection criteria performed
better than glmnetPC and the null model with covariates only when comparing
mean AUC on the test sets. The median number of predictors selected by
glmnetPC was four times higher than for glmnet when using CV for both methods.
Moreover, the variability in predictors selected was more important for
glmnetPC, as reported by an IQR value equal to 486, compared to 145 for our
method. pglmm with BIC selected 1 predictor (IQR: 1) compared to 16 (IQR: 145)
for pglmm with CV. This is consistent with our simulation results showing that
BIC results in sparser models with comparable predictive power. For high
cholesterol, very few genetic predictors were selected by all models, which
suggests that it may not be a highly polygenic trait. In fact, using only the
non-genetic covariates and first 10 PCs resulted in the best model for high
cholesterol based on mean test sets AUC.

##  4 Discussion

We have introduced a new method called pglmm based on regularized PQL
estimation, for selecting important predictors and estimating their effects in
high-dimensional GWAS data, accounting for population structure, close
relatedness and binary nature of the trait. Through a variety of simulations,
using both simulated and real genotype data, we showed that pglmm was markedly
better than a logistic lasso with PC adjustment when the number of
subpopulations was greater than the number of PCs included, or when a high
proportion of related subjects were present. We also showed that a lasso LMM
was unable to estimate predictor effects with accuracy for binary responses,
which greatly decreased its predictive performance. Performance assessment was
based on TPR of selected predictors, RMSE of estimated effects and AUC of
predictions. These results strongly advocate for using methods that explicitly
account for the binary nature of the trait while effectively controlling for
population structure and relatedness in genetic studies.

When the dimensionality of the confounding structure was low, we showed that
pglmm was equivalent to logistic lasso with PC adjustement. Hence, adjusting a
GLM with PCA is at best equivalent to pglmm, but with the additional burden of
selecting an appropriate number of PCs to retain in the model. Estimating the
dimensionality of real datasets, and thus the number of PCs to include as
fixed effects in a regression model can reveal to be challenging because
estimated eigenvalues have biased distributions (Yao and Ochoa,, 2022).
Another strategy involves selecting the appropriate number of PCs based on the
Tracy-Widom test (Tracy and Widom,, 1994). However, it is known that this test
tends to select a very large number of PCs (Lin and Zeng,, 2011), causing
convergence problems when fitting too many predictors. On the other hand,
modeling the population structure by using a random polygenic effect correctly
accounts for low and high-dimensional confounding structures, while only
fitting one extra variance component parameter.

We used real genotype data from the UK Biobank to simulate binary responses
and showed that BIC effectively selected sparser models with very high
precision and prediction accuracy, compared to AIC and CV. Using the same data
set, we illustrated the potential advantages of pglmm over a logistic lasso
with PC adjustment for constructing a PRS on two highly heritable binary
traits in a set of related individuals. Results showed that pglmm had higher
predictive performance for asthma, while also selecting consistently fewer
predictors as reported by median and IQR values.

In this study, we focused solely on the lasso as a regularization penalty for
the genetic markers effects. However, it is known that estimated effects by
lasso will have large biases because the resulting shrinkage is constant
irrespective of the magnitude of the effects. Alternative regularization like
the Smoothly Clipped Absolute Deviation (SCAD) (Fan and Li,, 2001) penalty
function should be explored. Although, we note that the number of nonzero
coefficients in the SCAD estimates is no longer an unbiased estimate of its
degrees of freedom. Other alternatives include implementation of the relaxed
lasso, which has shown to produce sparser models with equal or lower
prediction loss than the regular lasso estimator for high-dimensional data
(Meinshausen,, 2007). It would also be of interest to explore if tuning of the
generalized ridge penalty term on the random effects that arises from the PQL
loss could result in better predictive performance.

A limitation of pglmm compared to a logistic lasso with PC adjustment is the
computational cost of performing multiple matrix calculations that comes from
the estimation of variance components under the null. Indeed, at each
iteration, we perform a matrix inversion based on Cholesky decomposition with
complexity of $O(n^{3})$ and matrix multiplications with complexity of
$O(mn^{2}+S^{2}n^{2}+p^{2}n)$, where $n$ is the sample size, $m$ is the number
of non-genetic covariates, and $S$ is the number of variance components. Then,
we need to perform a spectral decomposition of the covariance matrix with a
computation time $O(n^{3})$. These computations become prohibitive for large
cohorts such as the full UK Biobank with a total of 500$K$ samples. A solution
to explore to increase computation speed and decrease memory usage would be
the use of conjugate gradient methods with a diagonal preconditioner matrix,
as proposed by Zhou et al., (2018).

Finally, we can take advantage of the fact that it is possible to allow for
multiple random effects to account for complex sampling designs and extend
pglmm to a wide variety of models. For example, building a PRS for a bivariate
binary trait, explicitly accounting for the shared causal pathways of many
diseases or complex traits. Moreover, pglmm could be used in models where
there is interest in selecting over fixed genetic and gene-environment
interaction (GEI) effects. Due to the hierarchical structure between the main
genetic and GEI effects, we will have to consider using a lasso for
hierarchical structures (Zemlianskaia et al.,, 2022).

## Software

Our Julia package called PenalizedGLMM is available on Github
https://github.com/julstpierre/PenalizedGLMM.

## Funding

This work was supported by the Fonds de recherche Québec-Santé [267074 to
K.O.]; and the Natural Sciences and Engineering Research Council of Canada
[RGPIN-2019-06727 to K.O., RGPIN-2020-05133 to S.B.].

## Acknowledgments

This research has been conducted using the UK Biobank Resource under
Application Number 20802. This study was enabled in part by support provided
by Calcul Québec (https://www.calculquebec.ca/) and Compute Canada
(https://www.computecanada.ca/). We thank the UK Biobank and all participants
for providing information.

## Data availability statement

Simulated data are available on Github
https://github.com/julstpierre/PenalizedGLMM/data. UK Biobank data are
available via application directly to UK Biobank
(https://www.ukbiobank.ac.uk/enable-your-research). The current study was
conducted under UK Biobank application number 20802.

Conflict of Interest: None declared.

## References

  * Akaike, (1998) Akaike, H. (1998).  Information Theory and an Extension of the Maximum Likelihood Principle, pages 199–213.  Springer New York, New York, NY.
  * Anderson, (2008) Anderson, G. P. (2008).  Endotyping asthma: new insights into key pathogenic mechanisms in a complex, heterogeneous disease.  The Lancet, 372(9643):1107–1119.
  * Astle and Balding, (2009) Astle, W. and Balding, D. J. (2009).  Population structure and cryptic relatedness in genetic association studies.  Statistical Science, 24(4):451–471.
  * Bezanson et al., (2017) Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017).  Julia: A fresh approach to numerical computing.  SIAM review, 59(1):65–98.
  * (5) Bhatnagar, S. R., Yang, Y., and Greenwood, C. M. T. (2020a).  ggmix: Variable selection in linear mixed models for snp data. R package version 0.0.1.
  * (6) Bhatnagar, S. R., Yang, Y., Lu, T., Schurr, E., Loredo-Osti, J., Forest, M., Oualkacha, K., and Greenwood, C. M. T. (2020b).  Simultaneous SNP selection and adjustment for population structure in high dimensional prediction models.  PLOS Genetics, 16(5):e1008766.
  * Böhning and Lindsay, (1988) Böhning, D. and Lindsay, B. G. (1988).  Monotonicity of quadratic-approximation algorithms.  Annals of the Institute of Statistical Mathematics, 40(4):641–663.
  * Breslow and Clayton, (1993) Breslow, N. E. and Clayton, D. G. (1993).  Approximate Inference in Generalized Linear Mixed Models.  Journal of the American Statistical Association, 88(421):9–25.
  * Bycroft et al., (2018) Bycroft, C., Freeman, C., Petkova, D., Band, G., Elliott, L. T., Sharp, K., Motyer, A., Vukcevic, D., Delaneau, O., O’Connell, J., Cortes, A., Welsh, S., Young, A., Effingham, M., McVean, G., Leslie, S., Allen, N., Donnelly, P., and Marchini, J. (2018).  The UK biobank resource with deep phenotyping and genomic data.  Nature, 562(7726):203–209.
  * Chen et al., (2016) Chen, H., Wang, C., Conomos, M. P., Stilp, A. M., Li, Z., Sofer, T., Szpiro, A. A., Chen, W., Brehm, J. M., Celedón, J. C., Redline, S., Papanicolaou, G. J., Thornton, T. A., Laurie, C. C., Rice, K., and Lin, X. (2016).  Control for population structure and relatedness for binary traits in genetic association studies via logistic mixed models.  The American Journal of Human Genetics, 98(4):653–666.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360.
  * Friedman et al., (2010) Friedman, J., Hastie, T., and Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1):1–22.
  * Gilmour et al., (1995) Gilmour, A. R., Thompson, R., and Cullis, B. R. (1995).  Average information REML: An efficient algorithm for variance parameter estimation in linear mixed models.  Biometrics, 51(4):1440.
  * Groll and Tutz, (2014) Groll, A. and Tutz, G. (2014).  Variable selection for generalized linear mixed models by L 1-penalized estimation.  Statistics and Computing, 24(2):137–154.
  * Hoffman, (2013) Hoffman, G. E. (2013).  Correcting for population structure and kinship using the linear mixed model: Theory and extensions.  PLoS ONE, 8(10):e75707.
  * Hoggart et al., (2008) Hoggart, C. J., Whittaker, J. C., Iorio, M. D., and Balding, D. J. (2008).  Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies.  PLoS Genetics, 4(7):e1000130.
  * Hui et al., (2017) Hui, F. K. C., Müller, S., and Welsh, A. H. (2017).  Joint Selection in Mixed Models using Regularized PQL.  Journal of the American Statistical Association, 112(519):1323–1333.
  * Kathiresan et al., (2008) Kathiresan, S., Melander, O., Anevski, D., Guiducci, C., Burtt, N. P., Roos, C., Hirschhorn, J. N., Berglund, G., Hedblad, B., Groop, L., Altshuler, D. M., Newton-Cheh, C., and Orho-Melander, M. (2008).  Polymorphisms Associated with Cholesterol and Risk of Cardiovascular Events.  New England Journal of Medicine, 358(12):1240–1249.
  * Lin and Zeng, (2011) Lin, D. Y. and Zeng, D. (2011).  Correcting for population stratification in genomewide association studies.  Journal of the American Statistical Association, 106(495):997–1008.
  * Loh et al., (2018) Loh, P.-R., Kichaev, G., Gazal, S., Schoech, A. P., and Price, A. L. (2018).  Mixed-model association for biobank-scale datasets.  Nature Genetics, 50(7):906–908.
  * Manolio et al., (2009) Manolio, T. A., Collins, F. S., Cox, N. J., Goldstein, D. B., Hindorff, L. A., Hunter, D. J., McCarthy, M. I., Ramos, E. M., Cardon, L. R., Chakravarti, A., Cho, J. H., Guttmacher, A. E., Kong, A., Kruglyak, L., Mardis, E., Rotimi, C. N., Slatkin, M., Valle, D., Whittemore, A. S., Boehnke, M., Clark, A. G., Eichler, E. E., Gibson, G., Haines, J. L., Mackay, T. F. C., McCarroll, S. A., and Visscher, P. M. (2009).  Finding the missing heritability of complex diseases.  Nature, 461(7265):747–753.
  * Meinshausen, (2007) Meinshausen, N. (2007).  Relaxed lasso.  Computational Statistics & Data Analysis, 52(1):374–393.
  * Novembre and Stephens, (2008) Novembre, J. and Stephens, M. (2008).  Interpreting principal component analyses of spatial population genetic variation.  Nature Genetics, 40(5):646–649.
  * (24) Ochoa, A. and Storey, J. D. (2016a).  FST and kinship for arbitrary population structures i: Generalized definitions.
  * (25) Ochoa, A. and Storey, J. D. (2016b).  FST and kinship for arbitrary population structures II: Method-of-moments estimators.
  * Price et al., (2006) Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006).  Principal components analysis corrects for stratification in genome-wide association studies.  Nature Genetics, 38(8):904–909.
  * Price et al., (2010) Price, A. L., Zaitlen, N. A., Reich, D., and Patterson, N. (2010).  New approaches to population stratification in genome-wide association studies.  Nature Reviews Genetics, 11(7):459–463.
  * Privé et al., (2019) Privé, F., Aschard, H., and Blum, M. G. B. (2019).  Efficient implementation of penalized regression for genetic risk prediction.  Genetics, 212(1):65–74.
  * Rabinowicz and Rosset, (2020) Rabinowicz, A. and Rosset, S. (2020).  Cross-validation for correlated data.  Journal of the American Statistical Association, pages 1–14.
  * Rakitsch et al., (2012) Rakitsch, B., Lippert, C., Stegle, O., and Borgwardt, K. (2012).  A lasso multi-marker mixed model for association mapping with population structure correction.  Bioinformatics, 29(2):206–214.
  * Reisetter and Breheny, (2021) Reisetter, A. C. and Breheny, P. (2021).  Penalized linear mixed models for structured genetic data.  Genetic Epidemiology.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2).
  * Sul et al., (2018) Sul, J. H., Martin, L. S., and Eskin, E. (2018).  Population structure in genetic studies: Confounding factors and mixed models.  PLOS Genetics, 14(12):e1007309.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288.
  * Tracy and Widom, (1994) Tracy, C. A. and Widom, H. (1994).  Level-spacing distributions and the airy kernel.  Communications in Mathematical Physics, 159(1):151–174.
  * Visscher et al., (2017) Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M. I., Brown, M. A., and Yang, J. (2017).  10 years of GWAS discovery: Biology, function, and translation.  The American Journal of Human Genetics, 101(1):5–22.
  * Wang et al., (2020) Wang, Y., Guo, J., Ni, G., Yang, J., Visscher, P. M., and Yengo, L. (2020).  Theoretical and empirical quantification of the accuracy of polygenic scores in ancestry divergent populations.  Nature Communications, 11(1).
  * Yang, (2005) Yang, Y. (2005).  Can the strengths of aic and bic be shared? a conflict between model indentification and regression estimation.  Biometrika, 92(4):937–950.
  * Yao and Ochoa, (2022) Yao, Y. and Ochoa, A. (2022).  Limitations of principal components in quantitative genetic association models for human studies.
  * Yu et al., (2005) Yu, J., Pressoir, G., Briggs, W. H., Bi, I. V., Yamasaki, M., Doebley, J. F., McMullen, M. D., Gaut, B. S., Nielsen, D. M., Holland, J. B., Kresovich, S., and Buckler, E. S. (2005).  A unified mixed-model method for association mapping that accounts for multiple levels of relatedness.  Nature Genetics, 38(2):203–208.
  * Zemlianskaia et al., (2022) Zemlianskaia, N., Gauderman, W. J., and Lewinger, J. P. (2022).  A scalable hierarchical lasso for gene-environment interactions.  Journal of Computational and Graphical Statistics, pages 1–36.
  * Zhao et al., (2018) Zhao, H., Mitra, N., Kanetsky, P. A., Nathanson, K. L., and Rebbeck, T. R. (2018).  A practical approach to adjusting for population stratification in genome-wide association studies: principal components and propensity scores (PCAPS).  Statistical Applications in Genetics and Molecular Biology, 17(6).
  * Zhou et al., (2018) Zhou, W., Nielsen, J. B., Fritsche, L. G., Dey, R., Gabrielsen, M. E., Wolford, B. N., LeFaive, J., VandeHaar, P., Gagliano, S. A., Gifford, A., Bastarache, L. A., Wei, W.-Q., Denny, J. C., Lin, M., Hveem, K., Kang, H. M., Abecasis, G. R., Willer, C. J., and Lee, S. (2018).  Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies.  Nature Genetics, 50(9):1335–1341.
  * Zou et al., (2007) Zou, H., Hastie, T., and Tibshirani, R. (2007).  On the “degrees of freedom” of the lasso.  The Annals of Statistics, 35(5).

## Tables

Table 1: Simulation parameters Parameter | Definition | Value
---|---|---
$M$ | Number of replications | 50
$h^{2}_{g}$ | Fraction of variance due to fixed | 0.5
| genetic effects (logit scale) |
$h^{2}_{b}$ | Fraction of variance due to random | 0.4
| genetic effects (logit scale) |
$\pi_{0}$ | Prevalence under the null | 0.1
$p$ | Number of snps | 5,000
$c$ | Fraction of causal SNPs | 0.01
Table 2: Mean and standard deviation of AUCs in test sets for 50 replications of the simulated genotype data. Model size represents the number of genetic predictors that are selected by each model. $K$ represents the number of intermediate subpopulations in the 1d linear admixture data, and the number of independent subpopulations in the independent data. $\%\Delta_{std}$ represents the relative decrease in AUC standard deviation for the predictions obtained by pglmm. |  | 1d linear admixture | independent
---|---|---|---
K | Model size | glmnetPC | pglmm | $\%\Delta_{std}$ | glmnetPC | pglmm | $\%\Delta_{std}$
10 | 5 | 0.765 (0.0456) | 0.769 (0.0443) | 3.0 | 0.801 (0.0496) | 0.802 (0.0494) | 0.4
| 10 | 0.790 (0.0350) | 0.794 (0.0344) | 1.5 | 0.817 (0.0445) | 0.817 (0.0454) | -2.1
| 15 | 0.804 (0.0313) | 0.808 (0.0305) | 2.6 | 0.826 (0.0417) | 0.827 (0.0425) | -2.0
| 20 | 0.814 (0.0272) | 0.817 (0.0275) | -1.1 | 0.831 (0.0404) | 0.832 (0.0418) | -3.7
| 25 | 0.820 (0.0253) | 0.821 (0.0262) | -3.4 | 0.834 (0.0395) | 0.835 (0.0409) | -3.7
| 30 | 0.823 (0.0247) | 0.824 (0.0248) | -0.5 | 0.836 (0.0390) | 0.837 (0.0401) | -2.8
| 35 | 0.825 (0.0243) | 0.827 (0.0245) | -0.7 | 0.838 (0.0386) | 0.839 (0.0400) | -3.6
| 40 | 0.827 (0.0241) | 0.828 (0.0242) | -0.3 | 0.840 (0.0382) | 0.840 (0.0395) | -3.3
| 45 | 0.829 (0.0239) | 0.830 (0.0238) | 0.2 | 0.841 (0.0381) | 0.842 (0.0394) | -3.5
| 50 | 0.830 (0.0238) | 0.831 (0.0238) | -0.1 | 0.842 (0.0380) | 0.843 (0.0390) | -2.7
20 | 5 | 0.751 (0.0431) | 0.764 (0.0419) | 2.7 | 0.771 (0.0430) | 0.807 (0.0387) | 9.9
| 10 | 0.775 (0.0383) | 0.788 (0.0358) | 6.5 | 0.789 (0.0387) | 0.822 (0.0355) | 8.3
| 15 | 0.789 (0.0356) | 0.802 (0.0313) | 12.1 | 0.801 (0.0375) | 0.830 (0.0333) | 11.0
| 20 | 0.798 (0.0336) | 0.811 (0.0301) | 10.4 | 0.808 (0.0368) | 0.835 (0.0316) | 14.0
| 25 | 0.803 (0.0327) | 0.816 (0.0299) | 8.5 | 0.815 (0.0367) | 0.838 (0.0308) | 16.1
| 30 | 0.807 (0.0321) | 0.819 (0.0295) | 8.0 | 0.819 (0.0361) | 0.840 (0.0305) | 15.6
| 35 | 0.810 (0.0315) | 0.821 (0.0297) | 5.6 | 0.822 (0.0354) | 0.842 (0.0303) | 14.4
| 40 | 0.812 (0.0310) | 0.823 (0.0293) | 5.5 | 0.826 (0.0346) | 0.843 (0.0301) | 13.0
| 45 | 0.814 (0.0309) | 0.824 (0.0293) | 5.1 | 0.829 (0.0341) | 0.844 (0.0298) | 12.4
| 50 | 0.816 (0.0302) | 0.825 (0.0290) | 3.9 | 0.831 (0.0336) | 0.845 (0.0297) | 11.9
Table 3: PRS results for the UK Biobank AUC values for asthma. We report mean of AUC and standard deviation for a total of 40 random splits. For model size, we report median and interquartile range for the number of genetic predictors selected. For pglmm, we compare performance when model is selected using BIC or CV. For BIC, the best model is chosen based on training fit. For CV, the best model is chosen based on maximum AUC on the validation set. Model | AUCval | AUCtest | Size
---|---|---|---
Asthma |  |  |
Covariates + 10PCs | 0.5232 (0.019) | 0.5254 (0.026) | -
glmnetPC (CV) | 0.5410 (0.018) | 0.5253 (0.027) | 67.5 (486)
pglmm (CV) | 0.5539 (0.023) | 0.5385 (0.026) | 16 (145)
pglmm (BIC) | - | 0.5452 (0.025) | 1 (1)
High cholesterol |  |  |
Covariates + 10PCs | 0.7183 (0.017) | 0.7215 (0.017) | -
glmnetPC (CV) | 0.7196 (0.018) | 0.7196 (0.018) | 0.5 (15.5)
pglmm (CV) | 0.7213 (0.018) | 0.7202 (0.019) | 3 (33.5)
pglmm (BIC) | - | 0.7212 (0.017) | 0 (0)

## Figures

Figure 1: Mean of 50 TPRs for the simulated genotype data. $K$ represents the
number of intermediate subpopulations in the 1d linear admixture data (left
panel), and the number of independent subpopulations in the independent data
(right panel).  Figure 2: Mean of 50 RMSEs for the simulated genotype data.
$K$ represents the number of intermediate subpopulations in the 1d linear
admixture data (left panel), and the number of independent subpopulations in
the independent data (right panel). Figure 3: Mean of 50 AUCs in test sets of
the simulated genotype data. $K$ represents the number of intermediate
subpopulations in the 1d linear admixture data (left panel), and the number of
independent subpopulations in the independent data (right panel).  Figure 4:
Mean of 50 AUCs, RMSEs and TPRs for the UK Biobank genotype data with related
subjects. Figure 5: Boxplots of the model selection simulations results for 50
replications of the UK Biobank genotype data with related subjects. For each
replication, the best model for pglmm was chosen using either AIC, BIC or CV.