To understand how galaxies reionized the Universe we must determine how the escape fraction of Lyman continuum LyC photons f SUBescSUB depends on galaxy properties. Using the z  0.3 Lowredshift Lyman Continuum Survey LzLCS we develop and analyze new multivariate predictors of f SUBescSUB. These predictions use the Cox proportional hazards model a survival analysis technique that incorporates both detections and upper limits. Our best model predicts the LzLCS f SUBescSUB detections with an rms scatter of 0.31 dex better than singlevariable correlations. According to ranking techniques the most important predictors of f SUBescSUB are the equivalent width EW of Lymanseries absorption lines and the UV dust attenuation which track lineofsight absorption due to H I and dust. The H I absorption EW is uniquely crucial for predicting f SUBescSUB for the strongest LyC emitters which show properties similar to weaker LyC emitters and whose high f SUBescSUB may therefore result from favorable orientation. In the absence of H I information star formation rate surface density SUBSFRSUB and O IIIO II ratio are the most predictive variables and highlight the connection between feedback and f SUBescSUB. We generate a model suitable for z gt 6 which uses only the UV slope SUBSFRSUB and O IIIO II. We find that SUBSFRSUB is more important in predicting f SUBescSUB at higher stellar masses whereas O IIIO II plays a greater role at lower masses. We also analyze predictions for other parameters such as the ionizingtononionizing flux ratio and Ly escape fraction. These multivariate models represent a promising tool for predicting f SUBescSUB at high redshift. SUPSUP Based on observations made with the NASAESA Hubble Space Telescope obtained at the Space Telescope Science Institute which is operated by the Association of Universities for Research in Astronomy Inc. under NASA contract NAS 526555. These observations are associated with programs GO15626 GO13744 GO14635 GO15341 and GO15639.
['Astrostatistics', 'Reionization', 'High-redshift galaxies', 'Starburst galaxies', 'Interstellar medium', 'Ultraviolet astronomy', 'Radiative transfer', '1882', '1383', '734', '1570', '847', '1736', '1335', 'Astrophysics - Astrophysics of Galaxies']
Multivariate Predictors of Lyman Continuum Escape. I. A Survival Analysis of the Lowredshift Lyman Continuum Survey
{'Multivariate Predictors of LyC Escape I: A Survival Analysis of the Low-redshift Lyman Continuum Survey ∗': "Anne E. Jaskot, 1 Anneliese C. Silveyra, 1, 2 Anna Plantinga, 3 Sophia R. Flury, 4 Matthew Hayes, 5 John Chisholm, 6 Timothy Heckman, 7 Laura Pentericci, 8 Daniel Schaerer, 9 Maxime Trebitsch, 10 Anne Verhamme, 9, 11 Cody Carr, 12, 13 Henry C. Ferguson, 14 Zhiyuan Ji, 15 Mauro Giavalisco, 4 Alaina Henry, 14 Rui Marques-Chaves, 9 Goran Ostlin, 5 Alberto Saldana-Lopez, 5 Claudia Scarlata, 16 G'abor Worseck, 17 and Xinfeng Xu 18  \n1 Department of Astronomy, Williams College, Williamstown, MA 01267, USA  \n2 Department of Physics, University of Nevada, Reno, NV 89557, USA  \n3 Department of Mathematics & Statistics, Williams College, Williamstown, MA 01267, USA  \n4 Department of Astronomy, University of Massachusetts Amherst, Amherst, MA 01002, USA  \n5 Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden  \n6 Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA  \n7 Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  \n8 INAF - Osservatorio Astronomico di Roma, via Frascati 33, 00078, Monteporzio Catone, Italy  \n9 Observatoire de Gen'eve, Universit'e de Gen'eve, Chemin Pegasi 51, 1290 Versoix, Switzerland  \n10 Astronomy, Kapteyn Astronomical Institute, Landleven 12, 9747 AD Groningen, The Netherlands  \n- Univ. Lyon, Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France  \n11  \n12 Center for Cosmology and Computational Astrophysics, Institute for Advanced Study in Physics, Zhejiang University, Hangzhou 310058, China  \n13 Institute of Astronomy, School of Physics, Zhejiang University, Hangzhou 310058, China  \n14 Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  \n15 Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  \n16 Minnesota Institute for Astrophysics, School of Physics and Astronomy, University of Minnesota, 316 Church St. SE, Minneapolis, MN 55455, USA  \n17 VDI/VDE Innovation+Technik, Berlin, Germany  \n18 Center for Interdisciplinary Exploration and Research in Astrophysics, Northwestern University, Evanston, IL 60201, USA", 'ABSTRACT': 'To understand how galaxies reionized the universe, we must determine how the escape fraction of Lyman Continuum (LyC) photons ( f esc ) depends on galaxy properties. Using the z ∼ 0 . 3 Low-redshift Lyman Continuum Survey (LzLCS), we develop and analyze new multivariate predictors of f esc . These predictions use the Cox proportional hazards model, a survival analysis technique that incorporates both detections and upper limits. Our best model predicts the LzLCS f esc detections with a rootmean-square (RMS) scatter of 0.31 dex, better than single-variable correlations. According to ranking techniques, the most important predictors of f esc are the equivalent width (EW) of Lyman-series absorption lines and the UV dust attenuation, which track line-of-sight absorption due to H i and dust. The H i absorption EW is uniquely crucial for predicting f esc for the strongest LyC emitters, which show properties similar to weaker LyC emitters and whose high f esc may therefore result from favorable orientation. In the absence of H i information, star formation rate surface density (Σ SFR ) and [O iii ]/[O ii ] ratio are the most predictive variables and highlight the connection between feedback and f esc . We generate a model suitable for z > 6, which uses only the UV slope, Σ SFR , and [O iii ]/[O ii ]. We find that Σ SFR is more important in predicting f esc at higher stellar masses, whereas [O iii ]/[O ii ] plays a greater role at lower masses. We also analyze predictions for other parameters, such as the  \nionizing-to-non ionizing flux ratio and Ly α escape fraction. These multivariate models represent a promising tool for predicting f esc at high redshift.', '1. INTRODUCTION': "Star-forming galaxies likely caused one of the most significant transformations in cosmic history: the reionization of hydrogen in the intergalactic medium (IGM) at z ≳ 6. Current constraints on the galaxy and quasar luminosity functions at z > 6 and on galaxies' ionizing photon production efficiencies favor stars as the dominant source of ionizing, Lyman continuum (LyC) photons (e.g., Bouwens et al. 2015, 2016; Finkelstein et al. 2015, 2019; Robertson et al. 2015; Ricci et al. 2017; Shen et al. 2020; Faucher-Gigu'ere 2020; De Barros et al. 2019; Endsley et al. 2021). Exactly which star-forming galaxies contribute these LyC photons remains unclear, however. Studies disagree as to whether low-, intermediate-, or high-luminosity galaxies dominate reionization (e.g., Razoumov & Sommer-Larsen 2010; Wise et al. 2014; Paardekooper et al. 2015; Finkelstein et al. 2019; Cain et al. 2021; Begley et al. 2022; Rosdahl et al. 2022; Saldana-Lopez et al. 2023; Wyithe & Loeb 2013; Naidu et al. 2020; Izotov et al. 2021; Ma et al. 2020; Matthee et al. 2022) or whether other properties such as concentrated star formation, nebular ionization, or starburst age demarcate the LyC-emitting galaxy population (e.g., Heckman et al. 2001; Clarke & Oey 2002; Alexandroff et al. 2015; Sharma et al. 2016; Marchi et al. 2018; Jaskot & Oey 2013; Nakajima & Ouchi 2014; Izotov et al. 2018b; Zastrow et al. 2013; Ma et al. 2015; Trebitsch et al. 2017; Naidu et al. 2022).  \nObservations with the James Webb Space Telescope (JWST) are revealing that proposed LyC-emitting galaxy populations exist at high redshift (e.g., Schaerer et al. 2022a; Williams et al. 2023; Endsley et al. 2023; Fujimoto et al. 2023; Tang et al. 2023; Mascia et al. 2024; Atek et al. 2024), although it remains to be seen whether they are sufficiently numerous at z > 6. However, identifying the galaxies responsible for reionization also requires knowing f esc , the fraction of LyC photons that escape into the IGM. The value of f esc for z > 6 galaxies is not known, and its physical connection with galaxy properties may be complex. LyC escape can depend on a galaxy's interstellar gas geometry, dust content, stellar feedback, and gravitational potential and may vary with time (e.g., Heckman et al. 2001, 2011; Wise et al. 2014; Ma et al. 2015; Sharma et al. 2016; Trebitsch et al. 2017; Chisholm et al. 2018; Barrow et al. 2020; Gazagnes et al. 2020; Mauerhofer et al. 2021; Saldana-Lopez et al. 2022).  \nLyC escape is an inherently multi-parameter problem.  \nBecause of attenuation in the IGM, detecting the LyC flux from galaxies becomes unlikely above z ∼ 4 (e.g.,  \nInoue et al. 2014). As a result, the astronomy community has relied on lower-redshift samples in order to investigate f esc and its dependence on galaxy properties. Studies of z ∼ 2 -4 galaxies find that f esc may be enhanced in galaxies with lower UV luminosities, lower dust attenuation, higher Ly α equivalent widths (EWs), strong [O iii ] λ 5007 emission, and/or compact sizes (e.g., Marchi et al. 2018; Steidel et al. 2018; Bassett et al. 2019; Fletcher et al. 2019; Nakajima et al. 2020; Begley et al. 2022; Saxena et al. 2022). More local samples, at z ∼ 0 . 3, likewise identify strong Ly α emission, elevated [O iii ]/[O ii ] ratios, high star formation rate surface densities (Σ SFR ), and low dust attenuation as characteristics of LyC emitters (LCEs) (e.g., Borthakur et al. 2014; Izotov et al. 2016b, 2018b; Verhamme et al. 2017; Chisholm et al. 2018). Nevertheless, f esc can have a range of values, even for galaxies with similar properties (e.g., Izotov et al. 2018b; Vanzella et al. 2016; Rutkowski et al. 2017; Fletcher et al. 2019; Bian & Fan 2020; Marques-Chaves et al. 2021).  \nThe recently completed Low-redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022a) offers a new opportunity to investigate all these physical properties simultaneously in a large statistical sample of LCEs and non-emitters. Through the LzLCS and archival programs (Izotov et al. 2016a,b, 2018a,b, 2021; Wang et al. 2019), 89 galaxies at z ∼ 0 . 3 have LyC observations from the Hubble Space Telescope (HST) Cosmic Origins Spectrograph (COS), and 50 of these galaxies are detected in the LyC. This combined sample, hereafter the LzLCS+, covers a wide range of luminosities, metallicities, H β and Ly α EWs, [O iii ]/[O ii ], and Σ SFR (Flury et al. 2022a), which enables it to more clearly reveal the trends and scatter between f esc and galaxy properties. By quantifying the relationship between f esc and observables, the LzLCS+ can generate predictions for f esc in z > 6 galaxies, where direct LyC detections are inaccessible.  \nThe results from the LzLCS+ confirm that f esc correlates with a variety of galaxy properties, from line-ofsight measurements such as H i covering fraction, dust attenuation, and Ly α escape fraction to global properties such as [O iii ]/[O ii ] and Σ SFR (e.g., Saldana-Lopez et al. 2022; Flury et al. 2022b; Chisholm et al. 2022; Xu et al. 2023). The latter properties hint at the possible role of mechanical or radiative feedback in creating low optical depth sightlines that allow ionizing photons to escape. Nevertheless, all correlations between f esc and observable properties show significant scatter  \n(Wang et al. 2021; Saldana-Lopez et al. 2022; Flury et al. 2022b; Chisholm et al. 2022; Xu et al. 2023).  \nA combination of properties might more accurately predict f esc than a single variable alone. For instance, several studies have generated multivariate predictions of Ly α properties for low- and high-redshift galaxies (e.g., Yang et al. 2017; Trainor et al. 2019; Runnholm et al. 2020). Runnholm et al. (2020) present such an analysis for Ly α luminosity using galaxies in the Lyman Alpha Reference Sample (LARS; Hayes et al. 2013; Ostlin et al. 2014). By applying multivariate regression using physical and observable galaxy properties, they predict the Ly α luminosities of the LARS galaxies with a root-mean-square residual scatter of 0.2-0.3 dex.  \nSimilarly, Maji et al. (2022) and Choustikov et al. (2024) perform multivariate linear fits for f esc using simulated galaxies from the SPHINX cosmological simulations (Rosdahl et al. 2018). Maji et al. (2022) find that a combination of four variables (escaping Ly α luminosity, gas mass, gas metallicity, and recent star formation rate) can explain 85% of the variance in escaping LyC luminosity, while three variables (escaping Ly α luminosity, recent star formation rate, and gas mass) explain 66% of the variance in f esc . However, some of the relevant variables (e.g., Ly α luminosity and gas mass) will be difficult or impossible to obtain for many galaxies at z > 6. In contrast, Choustikov et al. (2024) focus on multivariate f esc predictions using observable properties, such as [O iii ]/[O ii ], UV magnitude, and UV slope. While promising, the Maji et al. (2022) and Choustikov et al. (2024) simulation results require further testing and confirmation using observational samples.  \nOther recent studies have taken a more empirical approach by using the z ∼ 0 . 3 LyC observations to derive relationships between f esc and observable properties. Lin et al. (2024) fit for the probability of a galaxy having detectable vs. undetected f esc as a function of UVmagnitude, UV slope, and [O iii ]/[O ii ] and use these results to assess the likelihood of LyC escape from galaxies at z > 6. Mascia et al. (2023) develop a model that predicts f esc based on some of the strongest correlating variables in the LzLCS: β 1550 , UV half-light radius, and [O iii ]/[O ii ].  \nIn this paper, we investigate a variety of empirical multivariate predictions for f esc using the LzLCS+, the largest observational sample of LyC measurements at low redshift, which enables both reliable individual f esc measurements and comprehensive rest-frame UV and optical ancillary data. Whereas the Ly α data used by Runnholm et al. (2020) includes only Ly α -emitters and the simulated LyC data used by Maji et al. (2022) and Choustikov et al. (2024) include measurements down to  \nf esc = 0, the LzLCS+ sample includes 39 galaxies with non-detected LyC corresponding to f esc upper limits of 0.03-5.9%. Standard multivariate regression models do not account for upper limits, but a valid statistical analysis of the full LzLCS+ dataset should include information from both detections and non-detections (e.g., Isobe et al. 1986). Hence, we turn to the statistical approach of survival analysis, which can handle censored data like that of the LzLCS. We apply the Cox proportional hazards model (Cox 1972), a semi-parametric survival analysis method, to the full LzLCS+ dataset to derive multivariate predictions of f esc . By testing different sets of variables and assessing their predictive ability, we explore which physical properties combine to set a galaxy's f esc . We adopt a cosmology of H 0 = 70 km s -1 Mpc -1 , Ω m = 0 . 3, and Ω Λ = 0 . 7.", '2.1. Sample: The Low-redshift Lyman Continuum Survey': "To derive multivariate predictors of f esc , we use the Low-redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022a), the largest sample of LCEs at low redshift. Flury et al. (2022a) describe the survey sample, observations, data processing, and measurements in full, but we review some of the key details here. The LzLCS is a 134-orbit Cycle 26 HST program (GO15626; PI Jaskot) that obtained far-ultraviolet spectra, including the rest-frame LyC, with the COS G140L grating for 66 star-forming galaxies at z ∼ 0 . 3. The targeted galaxies each fulfill one or more hypothesized selection criteria for LCEs: high nebular ionization (O32 = [O iii ] λ 5007/[O ii ] λ 3727 ≥ 3), blue UV slopes (power law index β < -2), and/or concentrated star formation (Σ SFR > 0 . 1 M ⊙ yr -1 kpc -2 ). We combine the LzLCS with archival COS observations from Izotov et al. (2016a,b, 2018a,b, 2021) and Wang et al. (2019). We exclude one galaxy, J1333+6246 (Izotov et al. 2016b), whose [O iii ] λλ 5007,4959, H α , and H β line fluxes may be inaccurate; these lines appear truncated in the SDSS spectrum and the galaxy's H α /H β and H β /H γ line ratios are unphysical. Excluding this problematic galaxy, the combined LzLCS and archival galaxies (hereafter the LzLCS+ sample) consists of 88 low-redshift galaxies with LyC observations. The sample spans a range of 10 8 -10 10 M ⊙ in stellar mass, 12 + log 10 (O / H) = 7 . 5 -8 . 5, and -21.5 to -18.3 in observed (not corrected for internal reddening) 1500 ˚ A absolute magnitude.  \nWe process the raw COS spectra using the calcos pipeline (v3.3.9) and the FaintCOS software routines (Worseck et al. 2016; Makan et al. 2021) to reduce and calibrate the spectra and account for the backgrounds  \ndue to dark current and scattered geocoronal Ly α . We correct all spectra for Milky Way attenuation using the Green et al. (2018) dust maps and Fitzpatrick (1999) attenuation law. The LyC flux measurements represent the flux in a 20 ˚ A-wide region near rest-frame 900 ˚ A; we exclude wavelengths above an observed wavelength of 1180 ˚ A because of telluric contamination. For consistency, we re-process and re-measure the archival observations using this same methodology. Following the criteria in Flury et al. (2022a), we define LyC detections as observations that have a probability < 0 . 02275 of originating from background counts. By this definition, 49 of the 88 total galaxies have detected LyC. As in Flury et al. (2022a), we adopt the 84th percentile of the background count distribution as the upper limit for non-detections.  \nUsing these LyC fluxes, we then calculate the absolute f esc from different estimates of the intrinsic LyC, as described in Flury et al. (2022a); for this paper, we adopt the f esc values from Saldana-Lopez et al. (2022), derived from Starburst99 (Leitherer et al. 2011, 2014) spectral energy distribution (SED) fits to the UV continuum, following the methods outlined in Chisholm et al. (2019). As discussed in Flury et al. (2022a), the UV SED-fitting method is more reliable than H β -derived f esc estimates for the diverse LzLCS sample, because the presence of underlying, older ( > 10 Myr old) stellar populations that emit optical emission can bias the H β f esc results. The H β estimates of intrinsic LyC also do not account for LyC photons absorbed by dust and assume isotropic escape. With our preferred UV-based f esc estimate, the LzLCS+ LyC detections cover f esc = 0 . 4 -88 . 9% and the non-detections have upper limits ranging from 0.035.9%. We also consider an alternative, purely empirical measure of LyC escape, the F λ LyC /F λ 1100 flux ratio in Section § 4.1.  \nIn addition to LyC measurements, the LzLCS+ dataset contains a wealth of ancillary measurements, described in full in Flury et al. (2022a) and SaldanaLopez et al. (2022). Each galaxy in the survey has rest-frame optical photometry and spectroscopy through the Sloan Digital Sky Survey (SDSS; Blanton et al. 2017) and UV photometry via the Galaxy Evolution Explorer ( GALEX ; Martin et al. 2003). We derive stellar masses ( M ∗ ) using Prospector (Leja et al. 2017; Johnson et al. 2019) fits to the SDSS and GALEX photometry (Flury et al. 2022a, Ji et al. in prep.). We measure nebular line fluxes by fitting the SDSS spectra with multiple Gaussian profiles and then iteratively derive the electron temperature ( T e ), electron density, stellar absorption, and nebular E ( B -V ) from the [O iii ] λλ 5007,4959,4363, [S ii ] λλ 6716,6731, and Balmer  \nlines. Where the [S ii ] doublet or [O iii ] λ 4363 are undetected, we adopt n e = 100 cm -3 and the estimated [O iii ] λ 4363 flux from the Pilyugin et al. (2006) 'ffrelation', respectively. We then derive oxygen abundances via the direct method as implemented in pyneb (Luridiana et al. 2015). To estimate the star formation rate (SFR), we use the dust-corrected H β luminosities, Case B H α /H β ratio (Storey & Hummer 1995), and Kennicutt & Evans (2012) SFR calibration.  \nThe COS data provide estimates of additional physical parameters. We use the COS near-UV acquisition images to calculate the UV half-light radius. The sources all appear to be very compact within the central 1 '' diameter of the COS aperture, so that vignetting will not strongly affect our radius measurements. From the measured radius, we also derive Σ SFR as  \nΣ SFR = SFR 2 πr 2 50 , NUV . (1)  \nFrom the Starburst99 SED fits to the UV spectra, we fit for the dust excess E(B-V). We denote the E(BV) derived from the UV spectral fits as E(B-V) UV to distinguish it from E(B-V) neb derived from the Balmer line ratios. We also extrapolate the UV SED fits to infer the 'observed' (non-extinction corrected) absolute magnitude at 1500 ˚ A ( M 1500 ) and the power law index slope at 1550 ˚ A ( β 1550 ; Saldana-Lopez et al. 2022).  \nThe G140L spectra also cover Ly α . For our Ly α measurements, we use spectra extracted using a slightly wider spatial aperture (30 pixel vs. 25 pixel), because the scattered Ly α emission may be more spatially extended than the continuum light (Flury et al. 2022a). After masking the Si iii λ 1206 and N v λ 1240 regions, we linearly fit the continuum within 100 ˚ A of Ly α using an iterative sigma clipping algorithm and integrate the Ly α fluxes relative to this continuum level. This linear continuum fit typically agrees within < 5% with the continuum estimated from SED fits to the 25-pixel aperture extracted spectra (Saldana-Lopez et al. 2022). We do not correct for the small contribution of stellar Ly α absorption, but we conservatively adopt a 25% uncertainty on the Ly α continuum estimate (Flury et al. 2022a). From these fluxes, we derive Ly α luminosities, equivalent widths (EWs), and escape fractions ( f esc , Ly α ); we use the dust-corrected H β flux and the galaxies' measured electron temperatures and densities to infer the intrinsic Ly α flux. The reported Ly α measurements represent the sum of both absorption and emission along the line of sight. Even galaxies with net Ly α emission may have underlying Ly α absorption troughs, which are commonly detected in samples at z < 0 . 1 (e.g. McKinney et al. 2019; Hu et al. 2023). At the higher redshifts  \nof our sample, z ∼ 0 . 3, the COS aperture covers a larger physical area, and these troughs may experience more infilling due to scattered Ly α emission. Because the radii of our Ly α spectral apertures are ∼ 2.6 times larger than the UV continuum half-light radius (Flury et al. 2022a), we expect our Ly α measurements to capture most of the scattered emission (e.g. Hayes et al. 2013); future Ly α imaging of LzLCS+ galaxies will test this assumption (HST GO-17069, P.I. Hayes). Nine galaxies in our sample have significant Ly α absorption that overlaps with the Si iii λ 1206 absorption feature. To account for the uncertain strength of Si iii λ 1206 contamination in these galaxies, we re-measure the Ly α fluxes and EWs, excluding wavelengths within 500 km s -1 of Si iii λ 1206, and we increase our Ly α uncertainties to account for this possible flux difference. We note that this effect is minor, changing EWs by < 3 ˚ A and f esc , Ly α by ≤ 0 . 017, with a median change of only 0.001 in f esc , Ly α for the nine affected galaxies.  \nLastly, we use the G140L spectra to measure the EWs and residual fluxes near line center ( R l ) for a variety of low-ionization absorption lines, as described in SaldanaLopez et al. (2022). The UV absorption line measurements use the same (25 pixel aperture) spectral extraction as the LyC measurements, and the SED fits are used to estimate the continuum level. For optically thick gas, well-resolved absorption lines, and assuming a uniform dust screen the residual flux is related to the gas covering fraction ( C f ) by  \nR ( λ i ) = 1 -C f ( λ i ) (2)  \n(see e.g., Gazagnes et al. 2018 for a discussion of geometry effects on R ). This equation assumes that all lines are saturated; while this assumption appears to hold for most of the H i lines, the observational uncertainties are too high to draw conclusions about saturation for the metal lines (see Saldana-Lopez et al. 2022 for more information). In this work, we consider some of the strongest low-ionization metal absorption lines (Si ii λ 1260, C ii λ 1334, H i Ly β +O i ). To improve the S/N of the measurements, we also consider the inversevariance weighted average of the EWs and R l values for a given galaxy, calculated for the available H i Lymanseries lines (Ly β -Ly6) or for the available low-ionization metal (LIS) lines (Si ii λλ 1190 , 1193, Si ii λ 1260, [O i ] λ 1302, C ii λ 1334; see Saldana-Lopez et al. 2022 for details). We refer to these average absorption line measurements as EW(H i ,abs), R l (H i ,abs), EW(LIS), and R l (LIS). Because the measured residual intensity depends on resolution, we caution that the measured R l will not perfectly correspond to the true R l of the system; for the LzLCS+, this systematic error is of order  \n10-20% and is within the reported uncertainty (SaldanaLopez et al. 2022). While trends between R l and f esc will still be apparent (e.g., Saldana-Lopez et al. 2022), future investigations should use care when applying models derived from the LzLCS+ R l data to spectroscopic data with a different resolution.", '2.2. Survival Analysis': "As described above, one of the strengths of the LzLCS+ dataset is its large and diverse sample, spanning a wide range of measured f esc . The full sample of 88 galaxies includes 39 non-detections with f esc upper limits of only a few percent. These non-detections can convey essential information about which properties do and do not distinguish LCEs from non-LCEs, and omitting non-detections can bias fitted relations (c.f., Isobe et al. 1986). Standard multivariate linear regression does not account for censored data, that is, data with upper limits. Therefore, we instead apply survival analysis techniques to the LzLCS+ sample, as these techniques properly treat censored data.  \nAs implied in the name, survival analysis originated in the field of medicine (see reviews by Clark et al. 2003 and Bradburn et al. 2003, and see Feigelson & Nelson 1985 and Isobe et al. 1986 for applications to astronomy). In a medical context, the censored data often consist of known and unknown lifetimes for individuals participating in a medical study. People who are alive at the end of the study have an unknown lifetime; all we know is that they will live longer than their current age. Survival analysis techniques ultimately seek to describe the probability of a particular lifetime for a population or to compare how the survival probability changes between different populations, based on data that include some known lifetimes and some limits.  \nTo set up an analogous scenario for predicting the probability of a particular f esc , we do the following. Instead of a 'death' or known lifetime, we have a detection. Most survival analyses involve lower limits ('rightcensored' data). For simplicity in following common techniques, we therefore adopt the absorbed fraction of LyC ( f abs = 1 -f esc ) and its corresponding lower limits rather than f esc and its upper limits. With this setup, increasing f abs (decreasing f esc ) is equivalent to increasing the 'age' of the study participant. A non-detection, i.e., a lower limit to f abs (and upper limit on f esc ), is analogous to a lifetime greater than some threshold. Instead of starting at time = 0 and proceeding until we record a death or the study ends, we can imagine trying to observe a galaxy with f esc = 1 ( f abs = 0) and proceeding with deeper and deeper observations until we have a detection or cease our efforts. Just as the  \nthreshold between recording a measured lifetime vs. a limit will depend on factors such as the length of the study and age of the participant, the threshold between detected f esc and an upper limit will depend on aspects of our study such as the observation depth and galaxy brightness.", '2.2.1. The Cox Proportional Hazards Model': "The Cox proportional hazards regression model (Cox 1972) describes the probability of an event in an infinitesimally small window of time, given a combination of independent variables and assuming the event did not already occur. In the Cox model, this event probability (or 'hazard function') at time t for a set of variables x is modeled as  \nh ( t | x ) = h 0 ( t ) exp[ n ∑ i =1 b i ( x i -¯ x i )] , (3)  \nwhere b i are fitted coefficients for each variable x i and where h 0 ( t ) represents the baseline hazard, the probability of an event given average values of all variables (¯ x i ). An increase in variable x i relative to that variable's average value in the sample results in an exponential increase or decrease in the detection probability, depending on the sign of the coefficient b i . In our case, instead of modeling the probability of an event at time t , the hazard function represents the probability of a LyC detection in an infinitesimally small window of f esc for a set of independent variables, assuming the galaxy was not already detected at a higher value of f esc .  \nThe Cox model is semi-parametric. One of its strengths is that it does not assume a particular functional form for the baseline hazard function but rather estimates it non-parametrically. In other words, the event times or detection values do not have to obey a known statistical distribution (Bradburn et al. 2003). Because the Cox model predicts the detection probability at each possible f esc value, rather than predicting the value of f esc itself, the model results are identical for f esc and for log 10 ( f esc ). In other words, the non-parametric estimate of h 0 ( f abs ) for an input array of f esc values would be identical to h 0 (log 10 ( f abs )) for the equivalent input array of log 10 ( f esc ) values.  \nYet, the parametric nature of the hazard function equation (Equation 3) with its reported fitted coefficients makes it straightforward to analyze the relative importance of variables and generate predictions for future datasets. On the other hand, this fixed functional form implicitly assumes that the effect of the independent variables is multiplicative and does not depend on t , in our case, f abs . This functional form would not be appropriate, for example, if one variable tended to introduce a step function in f esc , while the others had more of  \na smooth correlation. However, the adopted functional form has some flexibility in that independent variables are allowed to take any functional form (e.g., x i can be replaced with log 10 ( x i ), x 2 i , etc). Another limitation of the Cox model is that, while it handles censored dependent variable values, it cannot handle censored data in the independent variables. Consequently, we have to limit our analysis to input variables, x i , that are available for all (or nearly all) of the sample galaxies, as we explain below. In summary, the Cox proportional hazards model serves as a multivariate regression model for describing and predicting dependent variables with upper limits. Many other common statistical methods either do not treat upper limits (e.g., multivariate linear regression, principal component analysis), are univariate analyses (e.g., Kaplan-Meier analysis), or are fully parametric, less-flexible models (e.g., the Weibull model).  \nIn fitting the Cox model to the LzLCS+, we experiment with different combinations of independent variables from the parent list given in Table 1. Where possible, we take the base 10 logarithm of the relevant variable. We cannot take the logarithm for variables that contain both positive and negative values, such as the UV LIS and Ly α lines; these lines range from absorption to emission within the sample. Putting variables on a logarithmic scale serves two purposes. First, it puts variables on a consistent scale of similar order of magnitude. In addition, the coefficients in Equation 3 have a simple interpretation; increasing a variable by an order of magnitude leads to a factor of e b i increase in the probability h ( f abs | x ). Exploring alternative functional forms for the input variables is outside the scope of this paper, but we test the effects of switching between logarithmic and linear scalings in the following subsection ( § 2.2.2). While some statistical methods require putting all variables on an identical scale, such as scaling from 0 to 1, this further rescaling is not necessary in the Cox model. Equation 3 effectively rescales variables; the sample mean is subtracted from each input variable measurement, and the b i parameter translates a linear change in variable x i to a change in the probability of measuring a given f esc . Any normalization would change the derived values of b i but not the statistical significance of input variables nor their ultimate effect on f esc .  \nIn selecting a subset of variables to use in our fits, we include multiple variables from each category in Table 1, but we avoid any variables that are closely related to each other and which may therefore be highly collinear. Collinearity can result in multiple best-fit solutions and can cause the Cox model to fail to converge. Consequently, we typically consider either Σ SFR alone or SFR and r 50 , NUV , but not all three variables at once. Sim-  \nlarly, we choose only one variable that traces nebular ionization, one variable to trace UV color, and limited subsets of the UV absorption line variables. We discuss the effect of selecting different variables in Section § 3.2. As mentioned above, the Cox proportional hazard model cannot account for missing or censored independent variable data. Hence, we exclude a galaxy from our fits if it is missing any of the required variables. Numbers in brackets in Table 1 indicate the number of galaxies missing a given measurement. Most of the models discussed in this paper exclude only one galaxy, J1046+5827, a non-leaker that lacks a reported Σ SFR and r 50 , NUV . We apply the Cox model to f esc as our primary dependent variable, but we also perform and discuss model fits to the F λ LyC /F λ 1100 ratio, ( § 4.1), LyC luminosity ( § 4.1), f esc , Ly α , and Ly α luminosity ( § 4.2).  \nWe apply the Cox proportional hazards model to our data using the python package lifelines (DavidsonPilon 2019). The lifelines CoxPHFitter routine returns the best-fit coefficients for our selected independent variable set, their p-values, and various goodnessof-fit statistics. As shown by Equation 3, the Cox proportional hazards model gives the probability of observing a particular f esc value given a set of physical or observable properties. Instead, our goal is to predict the expected value of f esc given a set of properties. As described below, to find the expected f esc , we find the median of the probability distribution, where the true f esc has an equal probability of lying above or below this adopted value (e.g., Bradburn et al. 2003; DavidsonPilon 2019); f esc is predicted to be above the median 50% of the time and below it the other 50% of the time.  \nThe median f esc = 1 f abs represents the f esc value where the survival function S ( f abs ), the probability that there is no detection at f abs , detect < f abs , reaches 0.5. This probability is equivalent to the probability that f esc , detect is detected at f esc , detect < f esc . The survival function is calculated as (e.g., Cox 1972; Bradburn et al. 2003; Davidson-Pilon 2019; McLernon et al. 2023)  \nS ( f abs ) = exp[ -HF 0 ( f abs ) · ph( x )] , (4)  \nwhere HF 0 is the baseline cumulative hazard function:  \nHF 0 ( f abs ) = ∫ f abs 0 h 0 ( f ) df (5)  \nand ph( x ) is the partial hazard function for a set of variables x:  \nph( x ) = exp[ n ∑ i =1 b i ( x i -¯ x i )] . (6)  \nOccasionally, S ( f abs ) for a set of parameters will never reach 0.5 and the predicted median f esc is indeterminate. This situation corresponds to an arbitrarily small predicted f esc ∼ 0.  \nThe best fit coefficients for the Cox proportional hazards model are found by maximizing the partial likelihood, which compares ph( x ) for each detection with the sum of ph( x ) for all galaxies with higher f abs (lower f esc ), including both detections and non-detections. The resulting coefficients b i are those that best sort f abs in order. Once the coefficients are determined, Breslow's estimator (Breslow 1972) determines the baseline cumulative hazard function, HF 0 ( f abs ), using the number of detections with values lower than each f abs (i.e., detections with higher values of f esc ) plus the ph( x ) values for all detections and limits higher than f abs (i.e., with lower f esc ). The model fits HF 0 ( f abs ) non-parametrically, reporting a value of HF 0 for each of the input f abs values of the LzLCS+ sample.  \nOne can also use the survival function to evaluate the values at which S reaches 0.159 and 0.841, the bounds corresponding to the Normal-theory 1 σ uncertainty of the predicted f esc . These probabilities account for the scatter in the correlations between f esc and the independent variables. The scatter may arise from observational uncertainty in the measurements as well as inherent variation among the population. In addition to this method, we have also used a Monte Carlo (MC) method to explore how the f esc predictions change if we vary each variable within its uncertainties. We randomly resample each independent variable measurement and each dependent variable detection using their uncertainties, re-run the Cox fit, and obtain new estimates of the median predicted f esc for each galaxy. Using the distribution of predicted f esc , we then calculate the 15.9 and 84.1 percentiles. For nearly all galaxies, the uncertainties estimated using the survival function are greater than or equal to the uncertainties determined from this MC method, which indicates that the scatter in the correlations is the dominant effect. We conclude that the survival function sufficiently represents the uncertainty in predicted f esc in most cases.", '2.2.2. Interpreting Cox Model f esc Predictions': 'As shown by Equation 4, the Cox model does not represent an equation that shows how f esc itself depends on particular variables. Rather, it describes how the probability of observing a given f esc changes for different sets of galaxy properties. The baseline cumulative hazard function HF 0 ( f abs ) describes the expected probability distribution for galaxies that have the average properties of the LzLCS+ dataset. In this case, x i = ¯ x i for all vari-  \nTable 1. Complete List of Independent Variables  \nNote -Numbers within brackets denote the number of galaxies within the LzLCS+ sample missing these measurements. Positive values of EW(H β ) and EW(Ly α ) denote net emission; positive values of the UV absorption lines denote net absorption.  \nables, and the reference probability that the observed f esc , detect is less than a particular value of f esc is simply S ref ( f esc ) = exp[ -HF 0 ( f esc )]. If we increase one variable by an increment of 1, x i -¯ x i = 1, ph( x ) = exp( b i ), and the new probability S ( f esc ) = S ref ( f esc ) exp( b i ) . Because the probability S ref ( f esc ) < 1, raising it to the power of exp( b i ) decreases the probability, which means that low values of f esc , detect < f esc are less likely. The probability distribution therefore shifts such that the probability corresponding to S =0.5 occurs at higher f esc , resulting in a larger predicted median f esc . Changing a second variable, x j , by an increment of 1 changes the new probability again by a power of exp( b j ), such that S ( f esc ) = [ S ref ( f esc ) exp( b i ) ] exp( b j ) ].  \nWe illustrate an example of the survival function probabilities, the associated median f esc , and the dependence of these quantities on input variables in Figure 1. In the left panel, we show a simple Cox model with three input variables: β 1550 , log 10 (Σ SFR ), and log 10 (O32). Dots represent the predicted median f esc where S ( f esc )= 0 . 5 for each value of O32. For log 10 (O32)= 1, S ( f esc )= 0 . 5 at f esc = 0 . 005. Changing only the log 10 (O32) value by 1, from O32 = 1, the blue dotted line, to O32 = 10, the magenta long-dashed line, shifts S ( f esc = 0 . 005) from 0.5 down to 0 . 5 exp(0 . 996) = 0 . 153, where 0.996 is the best-fit coefficient b i for the log 10 (O32) variable. Thus, according to this model, high f esc values are more common among galaxies with high O32. In this example, we changed only one parameter, but changing the other variables at same time could either further shift the probability distribution to higher f esc or counteract the change caused by increasing O32. We note that, unlike an equation for f esc as a function of these input variables, which could conceivably reach unphysical values of f esc > 1 for extreme parameters, the probability dis-  \nibution can shift toward higher values of f esc , but it never extends beyond f esc = 1.  \nFigure 1 also shows how the predicted probability responds to linear vs. logarithmic variables. As outlined above, the probabilities shift by a power of exp( b i ) for a step size of 1 in x i , which corresponds to an order of magnitude increase if x i represents the base 10 logarithm of a variable. In the left panel, the model uses log 10 (O32) as the input variable, and in the right panel, we show the effect of changing the form of the input variable from log 10 (O32) to linear O32. The differences are most apparent where a large linear spacing (e.g., changing O32 from 10 to 25) is not equivalent to a large logarithmic spacing and vice versa (e.g., O32 from 0.1 to 1). However, the change in the predicted median f esc values is largely minor. In changing to a linear scale, f esc changes from 0 to 0.006 for O32 = 0 . 1 and from 0.027 to 0.059 for O32 = 25. For the other plotted values, the change in f esc is < 0 . 002. The best-fit coefficients b i and HF 0 ( f esc ) values also differ between the two models and are optimized to match the observed distribution of f esc in the LzLCS+. Consequently, although the choice of variable form does affect predictions for galaxies at the extremes of the input parameter space, our main results, including the overall quality of our model predictions, are not highly sensitive to the functional form of the input variables.  \nWe do not thoroughly explore functional forms here, but we experiment with changing the ISM absorption measurements (EW(H i ,abs) and R l (H i ,abs)) and Ly α variables (EW(Ly α ), L (Ly α ), f esc , Ly α ) to logarithmic forms. For the logarithmic Ly α variables, we exclude all galaxies with net Ly α absorption from the dataset. We find the same variables are statistically significant in our models, although their exact coefficients and pvalues necessarily change. The goodness-of-fit metrics  \n<!-- image -->  \nFigure 1. Examples of the survival function probabilities predicted from the Cox model. The left panel uses β 1550 , log 10 (Σ SFR ), and log 10 (O32) as input variables, and the right panel uses β 1550 , log 10 (Σ SFR ), and O32 instead of log 10 (O32). Lines show the probability that the detected f esc , detect is less than a particular value of f esc or equivalently, the probability that f abs , detect is not less than a given f abs . The lines are plotted at f esc values corresponding to the measured f esc for the LzLCS+ dataset, from which the baseline cumulative hazard function is estimated non-parametrically. Data points show the median predicted f esc , where S ( f esc )= 0 . 5. The different lines show predictions for different input values of O32. As O32 increases, the probability distribution shifts to higher values of f esc , and the extent of the shift depends on whether the variable is linear or logarithmic. The black solid line in the left panel (log 10 (O32)= -1) shows an example where S never reaches 0.5, indicating an inferred median f esc ∼ 0.  \n<!-- image -->  \n( § 2.3) are also similar, typically changing by ≤ 0 . 02. The change in variable form also does not affect our main conclusions regarding the most important ranked variables ( § 2.4), although it can result in minor changes in ranked variable order.  \nObservationally, for most variables, log 10 ( f esc ) does seem to change significantly with the logarithm of a variable (Flury et al. 2022b). These trends could indicate that f esc has a power law dependence on many variables, although, as noted above, such relationships could only apply over the range where the predicted f esc ≤ 1. One exception to this logarithmic dependence is the UV absorption line measurements, which do show a dependence of log 10 ( f esc ) on the linear form of the variables (Saldana-Lopez et al. 2022). This functional form may indicate a more complex dependence of f esc on the gas geometry. We find that our predictions using the linear EW(H i ,abs) better match the observed f esc by every metric compared to models using log 10 (EW(H i ,abs)), and we therefore choose to keep EW(H i ,abs) in the linear form. Models using the linear or logarithmic form of R l (H i ,abs) perform comparably well, but with the linear form better matching the observed f esc for the strongest LCEs. Given their comparable or improved performance and for consistency with the LIS line measurements, which can reach negative values, we therefore  \nkeep all absorption line variables in their linear form. We also choose to keep the Ly α measurements in a linear form so as to not bias our dataset by excluding nonLy α emitters. We emphasize that our fit quality and main results are not sensitive to these choices regarding variable form.  \nIn summary, the Cox model predicts the probability of observing particular f esc values for different galaxy populations and can identify which variables most affect that probability. Statistically significant variables signify that the f esc probability is highly responsive to an incremental change in that variable: either a linear or logarithmic increase, depending on the variable form. While this work identifies these significant variables, future theoretical or observational programs could endeavor to derive the exact functional dependence of f esc on these parameters.', '2.3. Model Assessment': "We assess the goodness-of-fit for our models in a few different ways. Our primary method is the concordance index (Harrell et al. 1982), which applies to censored data. The concordance index considers all possible pairs of data points and how the observed rank order of f esc compares to the order predicted by the model. The model classifies each pair as concordant if the data point  \nwith a higher observed f esc also has a higher predicted f esc and discordant if it does not. Pairs can also be tied, if their predicted f esc values are identical. For some pairs with upper limits, their rank is ambiguous. The concordance index is  \nC = n c +0 . 5 n t n c + n d + n t , (7)  \nwhere n c is the number of concordant pairs, n t is the number of tied pairs, and n d is the number of discordant pairs. The concordance index ranges from 0 (perfect disagreement) to 0.5 (perfectly random) to 1.0 (perfect concordance; e.g., Davidson-Pilon 2019). 1 We also calculate the R 2 statistic for the galaxies with f esc detections as  \nR 2 = 1 -∑ i ( y i -f i ) 2 ∑ i ( y i -¯ y ) 2 , (8)  \nwhere y i are the observed values of log 10 ( f esc ), ¯ y is their mean value, and f i are the model-predicted log 10 ( f esc ) values. Following Maji et al. (2022), we also calculate the adjusted R 2 :  \nR 2 adj = 1 -(1 -R 2 ) n -1 n -p -1 , (9)  \nwhere n is the number of data points and p is the number of independent variables in the model. The R 2 adj parameter measures whether additional variables improve the model fit more than would be expected by random chance. Lastly, we also calculate the root-mean-square (RMS) dispersion about the model predictions:  \nRMS = √ ∑ i ( y i -f i ) 2 n . (10)  \nWe can only calculate these three quantities, R 2 , R 2 adj , and RMS, for the detections in our sample. Although they cannot provide a complete picture of model performance, they can tell us how well our predictions work for the galaxies with known f esc and can quantify our model performance at the high f esc end. Each of these metrics correlates strongly with the others and with the concordance index, and our main results are insensitive to the exact metric used.  \nWe choose to evaluate both R 2 and RMS using log 10 ( f esc ) rather than the linear f esc , because the scatter in the predicted f esc values is approximately constant in logarithmic space over the full observed f esc  \nrange. A single RMS value is therefore representative of both low and high f esc galaxies. In contrast, in linear space, the RMS changes systematically from low to high observed f esc . For instance, for our fiducial model ( § 3.1), the RMS is 0.02 for galaxies with f esc in the lowest third of the sample, 0.16 for the next third, and 0.22 for the final third, whereas the logarithmic RMS is ∼ 0.3-0.4 dex across the full sample range. Because R 2 and RMS use log 10 ( f esc ), these statistics only include galaxies with both detected f esc and non-zero predicted f esc . In contrast, the concordance index C incorporates the full dataset, including galaxies with observed upper limits and those with predicted f esc =0.  \nWe also assess the predictive ability of our models using cross validation, applying the model to data that was not used in generating the model itself. We perform a k -fold cross validation, in which we randomly split the dataset into k groups. We combine k -1 of the groups to become our training set, which we use to fit the model. We then test the model on the remaining group using the C , R 2 , and RMS metrics described above. ( R 2 adj is often undefined, due to runs where n for the small test set is equal to p ). We then repeat this process until each group has been used as the test group with the remaining groups used as the training set. Following Runnholm et al. (2020), we adopt k = 3, repeat the group selection 100 times, and average the final results. A value of k = 3 gives us sufficient galaxies in the test set to calculate statistics, and the repetition of the process ensures our results are not sensitive to the exact group selected (Runnholm et al. 2020). We provide the results of these tests in our evaluation of the models and distinguish the metrics derived from this cross-validation analysis by the subscript 'CV'.", '2.4. Variable Selection': 'We explore different subsets of independent variables and discuss the resulting model fits in Section § 3.2. With these different selections, we test the effect of switching the absorption line probed, of including or excluding Ly α , and of limiting the variable set to observables accessible to JWST (Section § 3.4). However, we also wish to analyze the relative importance of each independent variable and identify the variable subsets that most accurately predict f esc . To evaluate the importance of each variable, we turn to the tools of forward and backward selection (e.g., Runnholm et al. 2020; Maji et al. 2022).  \nIn forward selection, we run the Cox model with each independent variable individually in turn and determine which variable provides a fit with the highest concordance. We select this variable and then combine it with each remaining variable in turn to determine which  \ncombination of two variables provides the highest concordance. We proceed in this manner until we have a ranked order of the most significant variables. Conversely, with backward selection, we start with the full set of independent variables, and remove each one in turn. The combination of variables that gives the highest concordance identifies the least important excluded variable. We exclude it and proceed with the remaining variables, identifying and removing the next-least important variable each time to obtain a rank-ordered list. Importantly, each round of forward or backward selection only ranks variables with respect to the current best-performing model. Hence, some variables may be poorly ranked if they provide similar information as an already-selected variable. Variable rankings can reveal the most important predictors of f esc , but a poor ranking does not necessarily imply that a variable is uncorrelated with f esc .  \nWhen we perform our forward and backward selections, we do not use the full list of variables in Table 1 but limit it in the following ways. Firstly, we include only the EW(H i ,abs), EW(LIS), R l (H i ,abs), and R l (LIS) measurements, which are averages of multiple lines, rather than the measurements of individual absorption lines. This choice avoids overly limiting our sample sizes, since 22 galaxies are missing one or more individual absorption line measurements. Secondly, we avoid highly collinear variables and variables that measure very similar properties. We include only Σ SFR rather than SFR and r 50 , NUV , we use only one ionization-sensitive line ratio (O32 = [O iii ]/[O ii ]), and we use only one measure of UV dust attenuation ( β 1550 or E(B-V) UV ). When included separately, alternative measures of these properties end up with the same or nearly the same ranks, which suggests that they are indeed largely interchangeable. We explore three different sets of variables in our rankings: one set includes β 1550 and all remaining variables, one set corresponds to the variables in our best-performing Cox model, and one set is limited to variables accessible at high redshift. We also test our final variable rankings using an MC method. We re-run the forward and backward selection processes after randomly sampling the independent and dependent variables 100 times using their uncertainties. We present the results of the forward and backward selection methods in Section § 3.3.', '3.1. Fiducial Model': "For our fiducial model, we choose the following variables: log 10 ( M ∗ ), M 1500 , log 10 (EW(H β )), E(B-V) neb , 12+log 10 (O/H), log 10 (O32), log 10 (Σ SFR ), E(B-V) UV ,  \nf esc , Ly α , EW(LIS). We present the fitted coefficients and significance for these variables in Table 2 and compare the model-predicted median f esc with the observations in Figure 2. We list the goodness-of-fit metrics for this model in Table 3. Statistically significant variables, with p -values ≤ 0 . 05 are f esc , Ly α , E(B-V) UV , log 10 (Σ SFR ), E(B-V) neb , and log 10 (O32). Not surprisingly, these same variables individually correlate or anticorrelate with f esc , as discussed in Flury et al. (2022b) and Saldana-Lopez et al. (2022). We choose to include both O32 and EW(H β ) in the fiducial model, even though they correlate with each other (e.g., Flury et al. 2022b), because they may have subtly different relationships with f esc ; a young starburst may have both high EW(H β ) and high O32, but high global f esc will tend to increase O32, while decreasing EW( Hβ ) (e.g., Zackrisson et al. 2013; Nakajima & Ouchi 2014; Jaskot & Ravindranath 2016). However, we find that EW(H β ) is not significant and its coefficient is near zero, which indicates that EW(H β ) does not contribute meaningful information to the fit beyond what O32 already provides. For the fiducial model, we choose EW(LIS) instead of EW(H i ), because EW(LIS) is less affected by IGM and circumgalactic medium opacity and is therefore a viable indirect indicator across a wider redshift range. We consider the effect of substituting EW(H i ) or other measurements of absorption line strength in § 3.2.  \nAlthough most variables show the behavior expected from their individual relationship with f esc , one exception is the nebular attenuation, E(B-V) neb , which correlates with f esc in the Cox model instead of anticorrelating as might be expected. However, the model already contains the UV E(B-V) and f esc , Ly α , two parameters that show strong trends with E(B-V) neb . The nebular E(B-V) therefore essentially operates as a second-order effect, where the nebular attenuation would increase f esc at a fixed value of line-of-sight UV attenuation or f esc , Ly α . This extra nebular E(B-V) correlation could indicate a trend with some other physical property not included in the model, such as starburst age, global gas content, or gas clumpiness. In the case of f esc , Ly α , E(B-V) neb may help quantify how much of the Ly α escape is due to a low H i optical depth. At fixed f esc , Ly α , a higher E(B-V) neb implies a greater contribution of dust and weaker contribution of H i to the Ly α absorption. As expected, if we remove both Ly α and E(B-V) UV from the model, the correlation with E(BV) neb disappears.  \nAs discussed in § 2.2.2, the best-fit coefficients in Table 2 quantify how the probability of observing a given f esc responds to an incremental change in each variable.  \nThe probability S of observing f esc , detect <f esc changes from the original probability S ref by S = S exp( b i ∆ x i ) ref for a change in variable x i . The coefficients in Table 2 indicate that changing f esc , Ly α by 0.1 results in raising the probability to the power of exp(0 . 7) = 2 . 0, leading to a lower probability of low f esc values (see § 2.2.2. Conversely, increasing E(B-V) UV by 0.1 increases the probability of low f esc values by an even greater factor; the probability is raised to exp( -1 . 49) ∼ 1 / 5. For the other statistically significant variables, a 0.1 change in E(BV) neb also raises the probability to the power of ∼ 2, while 0.1 changes in either log 10 (Σ SFR ) or log 10 (O32) raise the probability to powers of 1.2 and 1.3, respectively.  \nAs illustrated in Figure 2 and Table 3, the fiducial model predicts the observed f esc values for detections with an RMS scatter of 0.36 dex for both the full sample and when the dataset is split into test and training sets (RMS CV = 0.36; see § 2.3). The model's high concordance, C = 0 . 89 ( C CV = 0 . 85), shows that it successfully ranks galaxies on the basis of f esc . The R 2 , R 2 adj , and R 2 CV metrics are 0.60, 0.49, and 0.55 respectively. This successful fit shows that including multiple variables results in a substantial improvement over predictions using a single variable. Even for some of the best single-variable correlations, f esc spans 2 dex at a given galaxy property (e.g., Flury et al. 2022b; Chisholm et al. 2022).  \nThe improved accuracy of the fiducial multivariate model is driven by only a handful of variables. Table 2 shows that only four variables are statistically significant in the fiducial model, and limiting the model to only these four variables achieves comparable accuracy, with the R 2 , RMS, and C metrics changing by 0.010.02. As can be seen in Table 2, highly insignificant variables tend to have best-fit coefficients near 0, which means the variable plays almost no role in the prediction. Excluding these variables therefore has little effect on the model, and because these extra parameters have negligible effects, overfitting is not a major concern.  \nDespite some overall improvement, the scatter in the fiducial Cox model is still significant, which demonstrates the difficulty in predicting the line-of-sight f esc even given a large amount of information. The residuals in the bottom panel of Figure 2 also show that the fiducial model systematically underpredicts the value of f esc for the strongest LCEs, even though it does typically identify them as having f esc > 0 . 1. Mascia et al. (2023) find a similar result and suggest that it may arise from the limited number of strong LCEs in the LzLCS+. We  \nsuggest that this under-prediction arises from a different effect, where the strongest LCEs represent a subset of galaxies with favorably oriented low-column density channels but whose global f esc may actually be substantially lower. We elaborate on this possibility in §§ 3.2 and 3.6. Because LyC may escape from narrow channels, the chance orientation of a galaxy can introduce scatter in relationships between global galaxy properties and the measured line-of-sight f esc (e.g., Cen & Kimm 2015; Flury et al. 2022b; Seive et al. 2022). Although the fiducial model explicitly includes variables that trace line-ofsight properties, such as EW(LIS) and E(B-V) UV , the LIS metals may imperfectly trace the LyC-absorbing H i gas.  \nOther effects may also contribute to the observed scatter. Simulations show that f esc may fluctuate in time, which may introduce scatter between f esc and observables that are sensitive to a different timescale (e.g., Trebitsch et al. 2017; Barrow et al. 2020). Other possible causes of the model scatter and underprediction of f esc are uncertainties in the observed f esc values or a disconnect between the global properties we measure and the local properties of the LyC-emitting region (e.g., Martin et al. 2015; Kim et al. 2023). Rivera-Thorsen et al. (2022) point out that a small number of massive stars can potentially dominate the escaping LyC emission; the properties of this LyC source region may differ from both the global galaxy properties and the properties inferred from the non-ionizing UV spectrum. Finally, because different properties may regulate LyC escape in different types of galaxies (e.g., Flury et al. 2022b), one relation may not suffice to predict f esc in all types of galaxies. We explore this possibility further in Section § 3.5.", '3.2. Modifications to the Fiducial Model': "Modifications to the fiducial model demonstrate that properties sensitive to dust and H i column density are the most useful predictors of f esc . We summarize the performance of the fiducial model and compare it with different modifications in Table 3. We do not test all possible combinations of the variables in Table 1, but we experiment with modifying the fiducial model by swapping alternative measures of some properties and by dropping each statistically significant variable.  \nWe first substitute different absorption line tracers (Table 1) for EW(LIS). Of these different tracers, only EW(Ly β ), EW(H i ,abs), and R l (H i ,abs) give statistically significant coefficients with p < 0 . 05. Because the sample size changes slightly for the other low-ionization metal absorption line selections (see Table 1), we cannot compare the metrics for the other models in detail. However, the fit quality appears comparable to the fidu-  \nTable 2. Fiducial Model Coefficients  \nNote -Variables are listed in order of their p -value significance. Positive values of EW(LIS) represent net absorption; positive values of EW(H β ) represent net emission.  \ncial model (RMS = 0 . 33 -0 . 39 and C = 0 . 88 -0 . 89). Only EW(H i ,abs) clearly improves the model fit, with RMS = 0 . 31 and C = 0 . 91.  \nThe fit with EW(H i ,abs) (Figure 3) is also our bestperforming model, and it scores better than the fiducial model on all metrics (Table 3). Moreover, this model more accurately predicts f esc for strong LCEs (Figure 3). Hence, the f esc underprediction in the fiducial model arises, at least in part, from the fact that EW(LIS) and the other variables imperfectly trace the line-of-sight H i gas. The Cox model is predicting the median f esc for a given set of parameters, and the strongest LCEs may not be noticeably distinct from more moderate LCEs in the fiducial model's set of parameters. In contrast, the low EW(H i ,abs) values of the strongest leakers do distinguish them from almost all of the moderate LCEs (e.g., Saldana-Lopez et al. 2022). The exact conversion from EW(LIS) to H i absorption and f esc depends on metallicity and gas geometry, since the trace metals will only show up in absorption for a sufficiently large H i column. While a range of physical conditions could cause weak LIS absorption, weak H i absorption more directly indicates LyC escape. Lastly, the observational uncertainty in EW(LIS) is also higher than in EW(H i ,abs), and these uncertainties will be the most problematic for the strongest LCEs with the weakest absorption lines. In conclusion, while LIS lines can identify LCEs, H i absorption lines appear necessary to accurately determine their line-of-sight f esc .  \nLike EW(H i ,abs), Ly α measurements can contribute key information about H i column density. Indeed, Maji et al. (2022) find that Ly α luminosity, L (Ly α ), is the  \nmost important variable in predicting the LyC f esc from simulations. In our fiducial model, f esc , Ly α was the most statistically significant coefficient, consistent with observational and theoretical studies that find a connection between LyC and Ly α (e.g., Dijkstra et al. 2016; Verhamme et al. 2017; Steidel et al. 2018; Izotov et al. 2020). In the fiducial model, f esc , Ly α serves as the best estimate of the H i optical depth. When EW(H i ,abs) is substituted in the model, f esc , Ly α remains statistically significant, but its p-value drops from 6E-9 to 3E-3 and its best-fit coefficient declines from 7 to 4. Because it is less sensitive to scattered emission, EW(H i ,abs) more directly traces the line-of-sight H i optical depth, which reduces the model's reliance on f esc , Ly α . Without EW(H i ,abs), however, f esc , Ly α provides crucial information in the fiducial model. If we exclude it, the fit quality worsens by all metrics; the RMS rises from 0.36 to 0.46, C drops from 0.89 to 0.84, and R 2 drops from 0.60 to 0.35. Substituting L (Ly α ) for f esc , Ly α also decreases the quality of the fit (RMS = 0.44, C = 0.85, R 2 = 0.41), whereas EW(Ly α ) performs almost as well as f esc , Ly α (RMS = 0.39, C = 0.88, R 2 = 0.54). Interestingly, when included together, both EW(Ly α ) and f esc , Ly α show up as significant coefficients. The two variables may provide slightly different information relevant to LyC escape. While f esc , Ly α may be more directly linked to the fraction of escaping LyC, the EW provides additional information about starburst properties such as the intrinsic Ly α production and underlying continuum. However, including both EW(Ly α ) and f esc , Ly α results in marginal, if any, improvement in the fit. All metrics change by only 0-0.01.  \nOf all the variables in the fiducial model, excluding E(B-V) UV or f esc , Ly α has the most detrimental effect  \nFigure 2. The f esc predictions from the fiducial Cox model compared with the observed f esc values for the LzLCS+ sample. Red circles indicate detections, and blue triangles represent the 1 σ upper limit on non-detections. The dashed line shows a one-to-one correspondence. The black cross in the upper left indicates the median size of the uncertainties in the observed and predicted f esc . The uncertainties in the predicted f esc represent the 15.9 and 84.1 percentiles of the f esc probability distribution. We plot galaxies with predicted f esc ∼ 0 at the bottom of the figure with their predicted f esc set to 10 -3 . The bottom panel shows a plot of the residuals on a linear scale. This model predicts f esc for the LyC detections with an RMS scatter of 0.36 dex.  \n<!-- image -->  \non the model fit. Without information on the line-ofsight dust attenuation from E(B-V) UV , R 2 drops to 0.30 and the RMS scatter rises by ∼ 0.1 dex to 0.48. The concordance also drops slightly to C = 0 . 86. Excluding f esc , Ly α has the strongest negative impact on the concordance ( C = 0 . 84) and the second largest impact on R 2 and RMS. The concordance metric includes nondetections, whereas the R 2 and RMS parameters do not. Hence, the larger effect of f esc , Ly α on concordance may illustrate that even without information on dust, a lack of Ly α emission is generally sufficient to identify nonleakers in the LzLCS+, even if the presence of Ly α does not necessarily imply LyC escape. Conversely, predicting accurate f esc for LCEs requires measurements of the line-of-sight dust attenuation.  \nSubstituting alternate measures of dust attenuation or ionization has little effect on the fiducial model. Using β 1550 instead of E(B-V) UV changes the goodness-offit metrics by only 0-0.03, consistent with the observed  \ntight correlation between β 1550 and E(B-V) UV in the LzLCS+ (Chisholm et al. 2022). However, the β 1200 parameter is not as successful at tracing dust content (e.g., Chisholm et al. 2022) and worsens the fit quality by all metrics. Likewise, alternate measures of ionization such as [Ne iii ]/[O ii ] or [O ii ]/H β work as well as O32, with no change in C and minor ( ≤ 0.03) changes in R 2 and RMS. Of the three measures of ionization, [O ii ]/H β performs best in all metrics, but this improvement is marginal.", '3.3. Most Important Variables': "As discussed in Section § 2.4, we perform forward and backward selection to determine which variables have the greatest effect on the model fit quality. In Table 4, we show a ranked ordering of variables based on forward and backward selection for three representative models. The 'Full Model' includes β 1550 instead of E(B-V) UV and includes some variables that measure similar but not highly collinear properties, in order to  \nFigure 3. The f esc predictions from the Cox model improve when we include EW(H i ,abs) instead of EW(LIS). Symbols are the same as in Figure 2.  \n<!-- image -->  \ntest which ones perform best. For instance, we include f esc , Ly α , EW(Ly α ), and L (Ly α ), which are related but may have different relationships with f esc . Our second set of variables is the more limited list from the Fiducial+HI model, our best-performing Cox model. This model also differs from the Full Model by using E(BV) UV . The third model, the 'JWST Model' excludes absorption line measurements, which are difficult or impossible to measure for most galaxies at z > 6, and Ly α measurements, which are heavily affected by the IGM at z > 6; we discuss this model in Section § 3.4. We also list the mean MC ranks for each variable obtained by sampling the observational uncertainties and rerunning the ranking process 100 times. We plot the distribution of these ranks in Figure 4.  \nTable 3. Goodness-of-Fit for Cox f esc Models  \nT able 4 . V ariables b y Rank ed Order  \n<!-- image -->  \na After this v ariable, add ing v ariables impro v es C b y < 0 . 01.  \nNote -V ariables in rank order from 1 (most imp ortan t) to 15 (least imp ortan t) b y forw ard or bac kw ard selection. Num b ers in brac k ets indicate the mean of the v ariable's ranks in the 100 MC runs. W e plot the full distribution of ranks in the MC runs in Figure 4 . The F ull Mo del includes most v ariables in T able 1 aside from those that are highly collinear. The Fiducial +HI Mo de l includes the v ariables in the b est-p erforming Co x Mo del. The JWST Mo del includes only v ariables acce ssible at z > 6.  \nWe find that the most important factors for predicting f esc are almost always the strength of H i absorption (EW(H i ,abs)) and dust attenuation (E(B-V) UV or β 1550 ). These factors are naturally connected to the lineof-sight f esc , since the two sources of LyC absorption are H i atoms and dust. Previous studies have used the combination of H i absorption lines and dust to derive predictions for f esc (e.g., Reddy et al. 2016; Chisholm et al. 2018; Saldana-Lopez et al. 2022), and the β 1550 parameter alone predicts much of the variation in f esc (Chisholm et al. 2022). Importantly, unlike other tracers of H i content and dust such as Ly α and the nebular E(BV), the H i absorption lines and UV measures of dust attenuation are more direct tracers of the H i and dust absorption along the line of sight to the UV-emitting stars and can thus substantially aid the prediction of the lineof-sight f esc . The combination of EW(H i ,abs) and one of the UV dust measures alone achieves a concordance C = 0 . 84.  \nThe MC method also finds that these two parameters, H i absorption and UV dust attenuation, have the top mean ranking. The MC mean ranks for the other variables are not necessarily consistent with their rank order in Table 4, which reflects the fact that the ranked order varies considerably among different MC runs (Figure 4). Since C barely changes after adding the first three to five variables, the order of most variables may be largely random, consistent with the broad ranking distribution in Figure 4. As also noted in § 3.1, only the most significant variables dominate the model predictions, and limiting the model to this smaller set of variables gives comparable results to a model generated using a larger variable set. Consequently, we recommend using only the most statistically significant or top-ranked variables when using Cox models to predict f esc .  \nTable 4 and Figure 4 show that Ly α measurements are often some of the more important variables, although generally not as important as EW(H i ,abs). For all rankings that involve Ly α , including the MC rankings, a Ly α measurement, either EW(Ly α ), L (Ly α ), or f esc , Ly α , achieves a rank between 1-5. In the MC runs, for the Fiducial+HI model, f esc , Ly α is the third most important variable after EW(H i ,abs) and E(B-V) UV . In Section 3.2, we found that excluding f esc , Ly α had one of the most detrimental effects on the fiducial model. Crucially, the fiducial model lacks EW(H i ,abs), so f esc , Ly α provided essential information on the H i optical depth in its place. With EW(H i ,abs) included, f esc , Ly α plays a lesser role but is still ranked among the top variables.  \nThe rankings for the Full Model in Table 4 show that EW(H i ,abs) is a better predictor of f esc in the LzLCS+ sample than the residual intensity R l (H i ,abs). This  \nranking may result from the low resolution of the FUV spectra. Low spectral resolution, which depends on the spatial size of the source, will artificially increase the observed residual intensity, R l (H i ,abs). In contrast, the EW(H i ,abs) is less sensitive to spectral resolution and may serve as a more accurate measure of the lineof-sight gas. Indeed, Saldana-Lopez et al. (2022) find that at a fixed R l , stronger LCEs have lower EWs. At the same time, the observational uncertainties are also higher in R l , such that f esc trends may appear more readily with EW. To investigate the effect of observational uncertainty on the rankings, we compare each variable's median uncertainty with standard deviation of that variable in the LzLCS+ sample. If the measurement uncertainty is comparable to the standard deviation, we may not be able to discern trends with that variable across the LzLCS+. For most variables, the median uncertainty is much lower, ≲ 20% of the standard deviation. These variables span the full range of possible ranks and MC ranks, which suggests that we can distinguish variables that do and do not affect the f esc predictions. However, a few variables ( R l (H i ,abs), log 10 ( M ∗ ), EW(LIS), and R l (LIS)) have higher ratios of median uncertainty to standard deviation (0.50-1.03). These variables all have high MC ranks in the Full Model (6.73-10.81) and any genuine trends with f esc may be hidden by the uncertainties in their measurements. As a test, we insert log 10 ( f esc ) as a dummy variable, as it should correlate perfectly with itself. Its MC rank begins to deviate from 1.00 when we give it an uncertainty ≥ 0 . 6 times the standard deviation. To further test the effect of uncertainty, we re-run the MC rankings after doubling the uncertainty in the EW(H i ,abs) measurements, such that the EW(H i ,abs) and R l (H i ,abs) variables have similar ratios of uncertainty to standard deviation. The EW(H i ,abs) MC ranks of 1.21-4.90 increase to 2.6-7.2 when the uncertainty is doubled. We conclude that the higher uncertainty in R l (H i ,abs), log 10 ( M ∗ ), EW(LIS), and R l (LIS) may prevent us from determining the importance of these variables. Higher resolution data are necessary to test whether EW(H i ,abs) or R l (H i ,abs) better predicts f esc .  \nIn their multivariate analysis of f esc from cosmological simulations, Maji et al. (2022) also use forward and backward selection to rank the importance of variables. They find that the three most important predictors of f esc are L (Ly α ), the SFR, and galaxy gas mass. Our rankings share some broad similarities with these simulation results. We also find that gas content (as measured by EW(H i ,abs)), SFR (in the form of the UV luminosity), and Ly α emerge as important parameters. However, our sample and the Maji et al. (2022) sim-  \nForward MC Rank Full Model  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank Fiducial+HI  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank JWST Model  \n<!-- image -->  \n<!-- image -->  \nBackward MC Rank Full Model  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->", 'Backward MC Rank Fiducial+HI': '<!-- image -->  \nBackward MC Rank JWST Model  \n<!-- image -->  \nFigure 4. The distribution of variable rankings from 100 forward selection (left) and backward selection (right) MC runs, after resampling each variable using its observational uncertainties. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.  \n<!-- image -->  \nulated galaxies have some crucial differences. First, we measure the line-of-sight f esc , whereas Maji et al. (2022) measure the total global f esc through all sightlines; consequently, the relevant parameter for our analysis is the line-of-sight H i rather than the global H i . Secondly, the galaxies in the Maji et al. (2022) sample are much less massive than the LzLCS+ sample, with a median stellar mass of log 10 ( M ∗ / M ⊙ )=6.41 vs. 8.8. Dust may play a greater role in determining f esc in the higher mass, more enriched galaxies in our sample, which would account for its higher importance in our predictions.', '3.4. Predictions at z > 6 with JWST': "Although the fiducial model and the models using the top-ranked variables predict f esc reasonably well, we cannot apply these models to JWST observations of galaxies in the epoch of reionization. The GunnPeterson trough and Lyman-series IGM absorption at lower redshifts will prevent H i absorption line measurements, and the partially neutral IGM at z > 6 can suppress Ly α emission (e.g., Stark et al. 2011; Schenker et al. 2014). Measuring LIS absorption lines instead requires high signal-to-noise observations of the restUV continuum, which will be difficult for faint galaxies. Consequently, we explore alternative models that use parameters that can be derived from z > 6 JWST observations.  \nWe create a JWST model by modifying the fiducial model to exclude EW(LIS) and f esc , Ly α . We also choose β 1550 instead of E(B-V) UV as it is more easily inferred from observations without any required stellar population modeling (e.g., Chisholm et al. 2022) or assumptions about the dust-attenuation law. The JWST model thus includes the following variables: log 10 ( M ∗ ), M 1500 , log 10 (EW(H β )), E(B-V) neb , 12+log 10 (O/H), log 10 (O32), log 10 (Σ SFR ), β 1550 . We show the resulting Cox model fit in Figure 5a and list the goodness-of-fit metrics in Table 3. The scatter is noticeably higher by 0.1 dex than in the fiducial model, with RMS=0.47 dex. Like the fiducial model, the JWST model also tends to under-predict f esc in several of the strongest LCEs with f esc > 0 . 2. This reduced fit quality shows that information about the line-of-sight H i content from absorption lines and Ly α is essential to precisely predict f esc . Nevertheless, other observable properties can still provide rough f esc estimates and distinguish strong LCEs from non-leakers. Since the JWST model relies more on global properties rather than line-of-sight H i to predict f esc , the underprediction for the most extreme LCEs suggests that their global properties are not distinct from more moderate LCEs. The model-predicted f esc may indicate the typical f esc value for this combination  \nof parameters. Moreover, if the most extreme f esc occurs only along favorable, nearly transparent sightlines, the f esc predicted from global parameters could be a better estimate of these galaxies' global average f esc . The strongest LCEs in the LzLCS+ sample do tend to have high nebular EWs (EW(H β ≳ 200 ˚ A; Flury et al. 2022b), which shows that they cannot be devoid of absorbing gas in all directions.  \nOnly three variables are statistically significant coefficients in the JWST model: β 1550 , log 10 (Σ SFR ), and log 10 (O32). These same three parameters are ranked as the most important in forward and backward selection and have the lowest mean ranks from forward and backward selection with MC sampling (Table 4, Figure 4). When we limit the model to just these three variables, the fit quality is similar to or better than the JWST model with the full set of variables (Table 3). We also find that the C CV for the cross-validation analysis is closer to the C derived from the full sample; with fewer variables in the fit, even a smaller training sample can generate a reliable model. We display the predicted f esc from this model in Figure 5b. Using [Ne iii ]/[O ii ] instead of O32 results in a similar quality fit (Table 3) and may be more suitable for observations covering a limited spectral range, observations at z ≳ 11 where [O iii ] is redshifted out of the NIRSpec wavelength range, or observations with uncertain nebular dust attenuation (Levesque & Richardson 2014). These fits demonstrate that the multivariate Cox model can predict f esc for high-redshift galaxies using JWST observables. We apply Cox models to high-redshift galaxy samples in a forthcoming paper (Jaskot et al., in prep.), where we also provide all parameters needed to apply these models to future samples.  \nOur JWST model includes similar variables as the multivariate f esc predictions derived by Choustikov et al. (2024) from the SPHINX simulations, but we find a different dependence on some of these variables. Choustikov et al. (2024) use the following observables: UV slope β , E(B-V) neb , H β luminosity, EW(H β ), UV magnitude, R23=([O iii ] λλ 5007,4959+[O ii ] λ 3727)/H β , O32, and half-light radius. Like the Choustikov et al. (2024) model, we find that β 1550 is statistically significant and that high f esc is associated with blue UV slopes. We likewise find that O32 is statistically significant. However, in the LzLCS+ sample, O32 correlates with f esc , whereas the SPHINX simulations find an anti-correlation. Strong LCEs with high O32 do not appear within the SPHINX galaxy population (Choustikov et al. 2024). This disagreement may result from the different properties of the observed vs. simulated galaxy populations, such as the stellar mass range. Alterna-  \n<!-- image -->  \nFigure 5. The f esc predictions from JWST models with the full set of variables (a) and top-ranked variables only (b). Symbols are the same as in Figure 2. The scatter in these models is higher than in the fiducial model, but they still reproduce the observed f esc with an RMS scatter of 0.46-0.47 dex.  \n<!-- image -->  \nit may indicate the need for more efficient radiative feedback and/or resolving smaller-scale turbulent gas structure, to allow LyC escape from younger clusters (e.g., Kimm et al. 2019; Kakiichi & Gronke 2021; Choustikov et al. 2024).  \nOur model also disagrees with Choustikov et al. (2024) regarding the role of Σ SFR . Whereas we find that Σ SFR is one of the most important predictors of f esc in the JWST model and is statistically significant, Choustikov et al. (2024) find that compactness is insignificant and shows little relationship with f esc . One of main differences between the LzLCS+ and SPHINX galaxy datasets is the mass range probed; the LzLCS+ sample has a median mass of log 10 ( M ∗ / M ⊙ )= 8 . 8 and extends up to log 10 ( M ∗ / M ⊙ )= 10 . 8, whereas the majority of the SPHINX galaxies have masses log 10 ( M ∗ / M ⊙ ) < 9. At even lower masses, log 10 ( M ∗ / M ⊙ ) < 7, simulations suggest that LyC escape may occur in a more extended mode, driven by star clusters in low-column density regions in the outskirts of galaxies (Kostyuk et al. 2023). The LzLCS+ likewise hints that the link between Σ SFR and f esc may be mass dependent. Several of the LCEs with high Σ SFR in the LzLCS+ sample have high masses (log 10 ( M ∗ / M ⊙ )= 9 . 7 -10 . 5), and as we discuss in the following section, Σ SFR may be a particularly significant predictor of f esc in high-mass galaxies.", '3.5. Variations with Galaxy Properties': "Flury et al. (2022b) discuss the possibility that the LyC escape process may vary among different types of galaxies. For example, the dominant feedback mecha-  \nnism may shift from radiative feedback to mechanical supernova feedback as a starburst ages or in galaxies with different masses or metallicities (e.g., Jaskot et al. 2019; Kimm et al. 2019; Jecmen & Oey 2023). In addition, the strongest LCEs could have a nearly density-bounded gas geometry, whereas LyC photons may escape along narrow channels in weaker LCEs (e.g., Gazagnes et al. 2020; Flury et al. 2022b). We therefore explore whether the Cox model predictions change if we limit our analysis to different subpopulations of the larger LzLCS+ sample. We first divide our galaxy sample into two bins, above and below the median value of stellar mass, log 10 ( M ∗ / M ⊙ )= 8 . 8.  \nDifferent variables appear important for the low- vs. high-mass sample fits. For the fiducial model, Σ SFR is statistically significant in the high-mass fit only, with a coefficient value ∼ 5 × higher than that found in the low-mass fit. Conversely, M 1500 and O32 are only statistically significant in the low-mass fit, with coefficients ∼ 8 and ∼ 5 × stronger than in the high-mass fit. The same pattern occurs with the JWST models, with Σ SFR significant in the high-mass fits and M 1500 and O32 significant in the low-mass fits.  \nThe variable ranks by forward and backward selection (Tables 5-6, Figures 6-7) likewise suggest that different variables are important for the two subsamples. As with the full sample (Table 4), variables tracing dust attenuation, H i absorption, and Ly α emission are some of the most important variables for both the high- and low-mass subsamples. However, for the high-mass subsample (Table 5), log 10 (Σ SFR ) consistently ranks as the  \nmost important variable in all models and in its mean MC rank values. In the low-mass subsample (Table 6), log 10 (O32) actually appears at the bottom of the ranked list for the full set of variables (rank 15), but by mean MC rank, log 10 (O32) is one of the four most important variables in all models and is the most important variable in the JWST model runs. This seemingly discrepant rank may suggest that O32 provides similar information as some of the other top-ranked variables. Figure 7 shows that O32 replaces EW(H i ,abs) as the top-ranked variable ∼ 10-25% of the time, and high O32 may generally imply weak H i absorption. In the JWST model, when EW(H i ,abs) is not included, O32 takes over as the most important variable in a majority of the MC runs. To test whether the smaller sample size of the high-mass and low-mass subsets affects these results, we re-run the MC variable ranking using bootstrap resampling to change the selected sample. The resulting ranking distributions are similar to those in Figures 6-7 with only minor changes to the top four ranked variables. Our main conclusions remain unchanged. Σ SFR consistently holds the top MC rank for every model for the high-mass subset. Dust, H i absorption, Ly α and O32 are important for the low-mass subset, and O32 is the top-ranked variable for the low-mass JWST model.  \nThe separate high- and low-mass models do a slightly better job at reproducing the observed f esc values in LCEs than the full fiducial model with all galaxies included. The separate models show a slightly higher R 2 = 0 . 64 -0 . 70 and lower RMS= 0 . 28 -0 . 32. However, the R 2 adj values of 0.31-0.43 are lower than the full model R 2 adj = 0 . 49, which suggests the improvement mainly comes from fitting a smaller sample. One exception is the JWST Model, where the high-mass subsample predictions are substantially improved compared to the full model; R 2 rises to 0.54 from 0.29, R 2 adj is 0.23 instead of 0.14, the RMS drops from 0.47 to 0.34, and C increases from 0.83 to 0.90. The Σ SFR variable takes on greater importance in the JWST model, since it lacks H i and Ly α information. The high-mass subsample f esc depends strongly on Σ SFR , and the best-fit coefficient for Σ SFR increases by a factor of four compared to the full sample JWST model. Evidently Σ SFR is key for accurately predicting f esc in the higher mass galaxies within the LzLCS+.  \nInterestingly, the importance of Σ SFR vs. O32 for different subsamples only emerges when we split the sample by stellar mass, not by UV luminosity. Almost all the highly ionized galaxies, with O32 ≳ 5 are in the lower stellar mass bin. Some of these galaxies can be quite luminous for their mass, so they are not as distinguishable by luminosity. The galaxies with the highest O32 ratios  \nalso have high nebular emission line EWs, suggestive of young ages. Lower-mass galaxies may have more bursty star formation histories in general (e.g., Lee et al. 2007). Among the low-mass galaxies in the LzLCS+, the high O32, high EW galaxies may represent currently bursting galaxies with extremely young average ages for their UVemitting stellar populations (e.g., Izotov et al. 2011).  \nThe noticeably different roles of the Σ SFR and O32 variables in the high- vs. low-mass samples suggest that the causes of LyC escape may indeed differ between these galaxy populations. In the low-mass sample, radiative feedback in young starbursts may drive high LyC escape over wide opening angles (e.g., Flury et al. 2022b). Galaxies with high luminosities and highly ionized gas would have higher f esc . At the same time, if the youngest stellar population dominates the UV light output, the measured line-of-sight dust attenuation and H i absorption would also be representative of the conditions near the LyC-emitting stars and would be sufficient to constrain f esc . In contrast, escape in higher-mass galaxies, with potentially longer-duration starburst episodes may rely on the cumulative effect of supernova feedback punching holes in the interstellar medium (ISM; e.g., Flury et al. 2022b). The ability to drive these required outflows may be linked to the galaxy's Σ SFR (e.g., Heckman et al. 2001, 2011; Kim et al. 2020). If the size of these holes is sufficiently small, the measured dust attenuation and H i EW across the full starburst may not reflect the conditions at the LyC escape site, especially if the youngest LyC-emitting populations do not dominate the total UV light. This disconnect could account for the reduced importance of β 1550 and EW(H i ,abs) and enhanced importance of Σ SFR in the high-mass models.  \nT able 5 . Most Imp ortan t V ariables for High-Mass Sample (log 10 ( M ∗ ) ≥ 8 . 8)  \n7]  \n[1.1  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[1.10]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n]  \n[2.98  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[1.79]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[4.35]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n24]  \n[3.  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[4.61]  \n(O32)  \n10  \nlog  \n[5.06]  \n)  \n∗  \nM  \n(  \n10  \nlog  \n[4.37]  \nα  \nLy  \n,  \nesc  \nf  \n[3.68]  \nUV  \nE(B-V)  \n[7.17]  \n)  \nα  \n(Ly  \nL  \n[5.38]  \n,abs)  \ni  \nEW(H", 'JWST Mo del': "del  \nMo  \nFiducial+HI  \ndel  \nMo  \null  \nF  \nRank  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nBac  \nard  \norw  \nF  \nT op-Rank ed V ariables  \nOrder  \nRank  \nMC  \ny  \nb  \nariables  \nV  \nop  \nT  \n[2.93]  \n1550  \nβ  \n[2.64]  \n1550  \nβ  \n[4.33]  \n,abs)  \ni  \nEW(H  \n[3.53]  \nα  \nLy  \n,  \nesc  \nf  \n5]  \n[6.7  \n(O32)  \n10  \nlog  \n[5.34]  \n1550  \nβ  \n[5.14]  \nneb  \nE(B-V)  \n[5.36]  \n))  \nβ  \n(EW(H  \n10  \nlog  \n[4.76]  \nUV  \n)  \nE(B-V  \n3.99]  \n[  \n,abs)  \ni  \nEW(H  \n[7.41]  \n1550  \nβ  \n[6.18]  \nα  \nLy  \n,  \nesc  \nf  \nNote -Num b ers in brac k ets indicate the mean of the v ariable's ranks in the 100 MC runs.  \nForward MC Rank Full Model High Mass  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank Fiducial+HI High Mass  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank JWST Model High Mass  \n<!-- image -->  \n<!-- image -->  \nBackward MC Rank Full Model High Mass  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->", 'Backward MC Rank Fiducial+HI High Mass': '<!-- image -->  \n<!-- image -->', 'Backward MC Rank JWST Model High Mass': 'Figure 6. The distribution of MC variable rankings from forward selection (left) and backward selection (right) for galaxies with log 10 ( M ∗ / M ⊙ ) ≥ 8 . 8. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.  \n<!-- image -->', 'T able 6 . Most Imp ortan t V ariables for Lo w-Mass Sample (log 10 ( M ∗ ) < 8 . 8)': "bles  \naria  \nV  \ned  \nop-Rank  \nT  \nOrder  \nRank  \nMC  \ny  \nb  \nariables  \nV  \nop  \nT  \nNote -Num b ers in brac k ets indicate the mean of the v ariable's ranks in the 100 MC runs.  \nTo explore variations related to feedback, we also split the sample based on the median values of log 10 (Σ SFR )=0.77 and O32=3.4. While splitting the LzLCS+ sample on mass resulted in slightly different variable selections, the high O32 and high Σ SFR subsample fiducial fits more closely resemble the full model, likely because they contain nearly the same population of strong LCEs; the high Σ SFR and high O32 subsets each include all but two of the strong LCEs with f esc > 0 . 1. The two most important predictors of f esc remain H i absorption (either R l (H i ,abs) or EW(H i ,abs)) and dust attenuation. Interestingly, Σ SFR is no longer a statistically significant variable in the fiducial and JWST models for the high Σ SFR and high O32 subsamples. Σ SFR may be more useful for identifying non-leakers than in actually quantifying f esc . Indeed, LCEs in the LzLCS+ sample nearly all lie above a threshold value of Σ SFR = 10 M ⊙ yr -1 kpc -2 , but above this value, no discernible relationship exists between f esc and Σ SFR (Flury et al. 2022b).  \nFor the low Σ SFR subsample, only f esc , Ly α and E(BV) UV are statistically significant in the fiducial model, while these same variables plus E(B-V) neb are significant for the low O32 subset. As usual, in the variable rankings, EW(H i ,abs) and dust measurements are typically important variables. However, Ly α measurements also appear to be relevant predictors. For the Full Model variable set, L (Ly α ) is the top-ranked variable for the low O32 subsample. Similarly, for the low Σ SFR subsample, f esc , Ly α is the top-ranked variable for the Fiducial+HI model and the top-ranked variable in forward modeling for the Full Model. Because these subsamples contain few strong LCEs, this result emphasizes that Ly α emission may be useful in distinguishing weak LCEs from non-leakers.", '3.6. Outliers in f esc Trends': "Many studies have identified relationships between f esc and physical or observable parameters (e.g., Verhamme et al. 2017; Izotov et al. 2018b; Flury et al. 2022b; Saldana-Lopez et al. 2022; Chisholm et al. 2022), yet even the best trends show significant scatter, with individual galaxies showing f esc significantly above or below the expected value. Schaerer et al. (2022b) highlight one outlier galaxy in particular, J1248+4259. J1248+4259's high nebular ionization and strong Ly α emission resemble those of other strong low-redshift LCEs, yet its f esc is ≤ 1.3%. Here, we investigate whether the Cox model improves f esc predictions for J1248+4259 and other outlier galaxies.  \nWe select outliers based on the variables that show the strongest trends with f esc in the Flury et al.  \n(2022b) analysis of the LzLCS+ sample: r 50 , NUV , O32, Σ SFR / M ∗ , EW(Ly α ), Σ SFR , and f esc , Ly α . From the strong LCEs ( f esc ≥ 0 . 05) in the LzLCS+ sample, we choose the three with the highest r 50 , NUV , the three lowest in O32, the three lowest in Σ SFR / M ∗ , the three lowest in EW(Ly α ), the three lowest in Σ SFR , and the three lowest in f esc , Ly α . Conversely, we choose three weak or non-LCEs ( f esc < 0 . 05) with the lowest r 50 , NUV , highest O32, highest Σ SFR / M ∗ , highest EW(Ly α ), highest Σ SFR , or highest f esc , Ly α . This selection gives us a list of 10 strong LCEs and 13 weak or non-LCEs that are outliers by one or more selection. In other words, applying a simple, single-variable trend as in Flury et al. (2022b) would incorrectly identify these galaxies as strong LCEs or as non-leakers.  \nFigure 8a shows that the fiducial model modified to include EW(H i ,abs) predicts f esc for these 'outliers' quite well. In fact, the outlier predictions are generally no more or less accurate than the predictions for other galaxies in the LzLCS+ sample. For example, the outlier galaxy J1248+4259, which 'should' be a strong LCE based on its high O32 ratio and high EW(Ly α ), has a predicted f esc of 0.006 according to the model, consistent with its observed limit of f esc ≤ 0.013.  \nCrucially, the fiducial+EW(H i ,abs) model includes information about the optical depth along the line of sight, which is necessary to predict the f esc observed along the same line of sight. As discussed in Sections 3.2 and 3.3, EW(H i ,abs), f esc , Ly α , and E(B-V) UV , all of which are sightline dependent, have some of the largest effects on the accuracy of the model predictions. If we remove these and other sightline-dependent parameters (E(B-V) neb , M 1500 ) from the fit (Figure 8b), we can no longer reproduce the line-of-sight f esc of many galaxies. Figure 8b indicates that we would expect a different f esc for these galaxies based on their global properties. Hence, if the only cause of the scatter in Figure 8b is chance galaxy orientation, the predicted f esc should represent the global average across all sightlines. By this model, J1248+4259 should indeed be an LCE, although not an exceptionally strong one, with f esc = 0 . 051.  \nAt high f esc , the model without line-of-sight information seems to systematically under-predict the observed f esc , which may suggest that the global properties of the strongest LCEs are not distinct from those of weaker LCEs. The observed properties of the strongest LCEs support this interpretation. For example, the strongest LCEs all tend to have high O32 (Flury et al. 2022b), yet the median observed f esc for the 11 galaxies with O32 ≥ 10 is only 0.05. For these same 11 galaxies, the model without line-of-sight information finds the median of their predicted f esc values to be 0.04, compa-  \nForward MC Rank Full Model Low Mass  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank Fiducial+HI Low Mass  \n<!-- image -->  \n<!-- image -->  \nForward MC Rank JWST Model Low Mass  \n<!-- image -->  \n<!-- image -->  \nBackward MC Rank Full Model Low Mass  \n<!-- image -->  \n<!-- image -->  \nBackward MC Rank Fiducial+HI Low Mass  \n<!-- image -->  \n<!-- image -->  \nBackward MC Rank JWST Model Low Mass  \n<!-- image -->  \nFigure 7. The distribution of MC variable rankings from forward selection (left) and backward selection (right) for galaxies with log 10 ( M ∗ / M ⊙ ) < 8 . 8. Variable combinations are the same as in Table 4. Thick solid lines indicate the four variables with the lowest (best) mean MC ranks.  \n<!-- image -->  \n<!-- image -->  \nFigure 8. The f esc predictions from the fiducial model (a) modified to use EW(H i ,abs) and (b) modified to exclude information about the line of sight from Ly α , absorption lines, dust attenuation, or the observed UV luminosity. Black crosses show strong LCEs that are outliers by one or more selection criteria; gray circles and triangles show weak and non-LCEs that are outliers by one or more selection criteria. Small red circles and blue triangles represent the other detections and non-detections in the LzLCS+ sample. The fiducial+EW(H i ,abs) model reproduces the f esc of the outliers, whereas the model without line-of-sight information does not.  \n<!-- image -->  \nrable to the observations. The only model that successfully reproduces the f esc of the strongest LCEs is the fiducial+EW(H i ,abs) model. Hence, of the variables we consider, only EW(H i ,abs) distinguishes the strongest LCEs from other galaxies. Either the strongest LCEs differ in some other unknown global property or their high f esc results from a favorable orientation compared to weaker LCEs. The high O32 galaxies may have f esc ∼ 0 . 04 -0 . 05 along most sightlines but may also contain several extremely optically thin channels, which lead to much higher observed f esc , ≳ 0 . 5 in some cases.  \nAlthough orientation bias likely causes some of the scatter in f esc predictions, the Cox modeling results also suggest that some of the 'outlier' galaxies are not outliers at all, at least not when multiple global variables are considered. For example, J0919+4906 and J02320426 have high O32 ratios of 12.7 and 9.9 but f esc of 0.049 and < 0 . 047, respectively. However, both galaxies also have low stellar masses (log 10 ( M ∗ )= 7 . 5 and 7.4), such that their O32 values are not extreme compared to other galaxies of similar mass (see Figure 23 in Flury et al. 2022b). Furthermore, other criteria, such as their moderate or low Σ SFR values (log 10 (Σ SFR )= 0 . 66 and -0.39) would argue against a high f esc . The combination of multiple properties leads to a moderate predicted f esc for J0919+4906 ( f esc =0.048) and little escaping LyC ( f esc =0.005) for J0232-0426, even when the model includes only global, and not sightline-dependent, properties. Similarly, several galaxies have concentrated star formation (J0910+6105, J1209+3053, J1248+1234, J1349+5631, J1503+3644, J1517+3705) or strong Ly α (J0122+0520, J0811+4141, J0901+2119, J1648+4957) but f esc < 0 . 05; these galaxies' other properties yield consistent predictions of f esc < 0 . 05, even in models that neglect line-of-sight information.  \nThe Cox model's success at reproducing f esc even for these outliers demonstrates that accurately predicting f esc requires multiple input variables. Moreover, because the observed f esc depends on our line of sight, the best models also require sightline-specific information about optical depth from absorption lines, UV dust attenuation, and Ly α . Although Ly α and H i absorption lines are inaccessible for most z > 6 galaxies, highredshift f esc predictions can include information about UV dust attenuation as in the JWST models discussed in Section 3.4. Information about dust attenuation is critical to predict f esc , and large sample sizes can probe multiple random orientations to estimate the average f esc for a galaxy population (e.g., Cen & Kimm 2015; Saldana-Lopez et al. 2023).", '4.1. Alternative Measures of LyC: Predicting F λ LyC /F λ 1100 and L (LyC)': "As shown above, multivariate Cox models can successfully predict f esc ; we now consider whether they can predict alternative measurements of the escaping LyC. To model how galaxies reionized the universe, we need to know the input rate of LyC photons into the IGM from different galaxy populations. Deriving this quantity from f esc requires knowing the production rate of LyC photons in each galaxy, which is typically derived from SED fitting or Balmer line observations. To avoid this process and its associated systematics, we could instead predict a more directly observable quantity, such as the F λ LyC /F λ 1100 ratio or the total escaping LyC luminosity log 10 ( L (LyC)). The F λ LyC /F λ 1100 flux ratio measures the ratio of observed LyC flux to the observed flux at 1100 ˚ A. As such, it depends on observed quantities and does not require an assumed stellar population model or dust attenuation law. To predict F λ LyC /F λ 1100 using the Cox proportional hazards model, we follow the same procedure as for f esc and convert the F λ LyC /F λ 1100 measurements to right-censored data (i.e., detections and lower limits) by predicting 1F λ LyC /F λ 1100 . Unlike f esc and F λ LyC /F λ 1100 , which cannot be higher than 1, L (LyC) has no obvious maximum allowed value. To convert it to a right-censored format, we therefore adopt a maximum of log 10 ( L (LyC))= 41, greater than the maximum luminosity in the sample of log 10 ( L (LyC))=40.6, and predict 41-log 10 ( L (LyC). We present the goodnessof-fit metrics for the fiducial and JWST model predictions of F λ LyC /F λ 1100 , L (LyC), and Ly α measurements in Table 7.  \nThe F λ LyC /F λ 1100 models perform worse than the f esc models for every variable combination we consider (see Figure 9) and by all three metrics ( R 2 , RMS, C ). Although F λ LyC /F λ 1100 is simpler to measure, f esc appears simpler to predict. It may be more directly linked to galaxy physical properties such as optical depth, porosity, and feedback. Furthermore, while measuring F λ LyC /F λ 1100 does not depend on assumptions about the galaxy SED, predicting it does. The galaxy SED shape changes with age and metallicity, which affects the ratio of F λ LyC /F λ 1100 even for f esc = 1 (see e.g., Chisholm et al. 2019).  \nAlthough we derive both E(B-V) UV and f esc from the COS UV spectra, the greater accuracy of the f esc Cox models does not appear to result from this interdependency. First, we see the same improvement in accuracy for f esc vs. F λ LyC /F λ 1100 when we use β 1550 instead of E(B-V) UV ; β 1550 is directly related to the observed spectral slope and is not sensitive to the choice of stellar population model or star formation history (Chisholm et al.  \n2022). Secondly, the f esc models have higher R 2 and C than the F λ LyC /F λ 1100 predictions, even for the models that contain no UV spectral information at all. We conclude that although F λ LyC /F λ 1100 is a useful quantity, accurately predicting it is difficult given its dependence on the intrinsic SED shape. We list and discuss the topranked variables for the F λ LyC /F λ 1100 predictions in the Appendix.  \nThe F λ LyC /F λ 1100 predictions rely on much the same variables as the f esc predictions, except for a weaker dependence on the UV dust attenuation, as expected. For the fiducial model, the same variables are statistically significant except for E(B-V) UV , and excluding E(B-V) UV from the fiducial model changes R 2 , RMS, and C by only ≤ 0 . 03. However, when we substitute EW(H i ,abs) in the fiducial model, E(B-V) UV does appear statistically significant, indicating that its influence is not always negligible. Table 8 in the Appendix lists the top-ranked variables and likewise indicates a reduced role for dust attenuation in predicting F λ LyC /F λ 1100 compared to its role in predicting f esc .  \nLike F λ LyC /F λ 1100 , the predictions for L (LyC) (Figure 10, Table 7) are also less accurate than for the f esc predictions. The concordance drops for every model, and R 2 is negative in all cases, meaning that simply adopting the mean L (LyC) for the detections would result in a better fit than using the model predictions. The negative R 2 arises from the fact that the R 2 calculation (Equation 8) uses only the detections and their mean, whereas the Cox model fits to the full set of galaxies; R 2 becomes positive if we use the median L (LyC) from the full sample of limits and detections, instead of the mean L (LyC) of the detections alone. Not surprisingly, other measures of luminosity show up as important variables for the L (LyC) models. M 1500 is statistically significant in the fiducial model and is one of top three ranked parameters by forward and backward selection, while L (Ly α ) holds the top place in forward selection for the full variable list (see Table 9 in the Appendix for a list of top-ranked variables). However, the L (LyC) models' poor performance indicates that key information is missing. Specifically, a lack of detailed information about the SED likely prevents an estimate of the intrinsic L (LyC) from M 1500 . The escaping L (Ly α ) also imperfectly traces the escaping L (LyC) because of scattering. In addition, the production of the Ly α recombination line requires LyC absorption, which could cause L (Ly α ) to drop for the few extreme LCEs in the LzLCS+ sample (e.g., Nakajima & Ouchi 2014). In conclusion, neither the L (LyC) nor the F λ LyC /F λ 1100 mod-  \nide satisfactory alternatives to the f esc predictions.", '4.2. Predicting f esc , Ly α': 'Although closely related to the H i optical depth, the escape of Ly α may depend on different parameters than the LyC escape fraction does. Here, we apply the Cox model to the f esc , Ly α data for the LzLCS+ sample and assess the differences between predictions for f esc , Ly α vs. f esc . We use the same input variables, with one exception: we exclude any variables containing Ly α measurements ( f esc , Ly α , EW(Ly α ), L (Ly α )) from the models. Also, unlike our f esc data, the f esc , Ly α measurements include no upper limits, since Ly α can appear in absorption or emission. For the negative values of f esc , Ly α , representing net absorption, we cannot calculate log 10 ( f esc , Ly α ) and do not include these galaxies in the reported R 2 or RMS; they are included in the concordance.  \nThe fiducial model is not very successful at predicting f esc , Ly α (Figure 11a, Table 7), with R 2 = 0 . 33, R 2 adj = 0 . 24, RMS=0.32, C = 0 . 75. The lack of upper limits and need to precisely rank galaxies with even very low f esc , Ly α may explain the reduced C values for the f esc , Ly α model. Including EW(H i ,abs) improves the f esc , Ly α predictions ( R 2 = 0 . 44, RMS=0.30, C = 0 . 78; Figure 11b). As with the f esc predictions, the JWST f esc , Ly α model performs worse than the fiducial model (Figure 11c), which demonstrates that the observed f esc , Ly α does depend on line-of-sight H i properties such as the neutral gas covering fraction or H i column density. Unlike the f esc models, for f esc , Ly α , the fiducial model without E(B-V) UV performs comparably to the fiducial model ( R 2 = 0 . 37, RMS=0.29, C = 0 . 75), because Ly α is more sensitive to the nebular attenuation rather than the stellar attenuation.  \nWe also investigate L (Ly α ), using a similar method as for L (LyC) and predicting 43.5-log 10 ( L (Ly α )). We exclude all galaxies with net Ly α absorption and predict only the luminosities of Ly α Emitters. Whereas we found that L (LyC) was harder to predict than f esc , for L (Ly α ), the situation is reversed: the L (Ly α ) models are actually slightly better than the f esc , Ly α models (Table 7 and Figure 12). The L (LyC) predictions suffered from the large uncertainty in estimating the intrinsic LyC production from the FUV continuum. The available Balmer line measurements make estimating the intrinsic Ly α luminosity less of a concern. Another explanation for the better performance of the L (Ly α ) predictions may be mundane; since it depends on predicting the intrinsic Ly α , f esc , Ly α is a more modeldependent measurement than L (Ly α ) and consequently  \n<!-- image -->  \n(a)  \nFigure 9. Examples of F λ LyC /F λ 1100 model predictions using the variable sets for the fiducial model (a) and the JWST model (b). Symbols are the same as in Figure 2. Predictions of F λ LyC /F λ 1100 are not as accurate as the model predictions for f esc .  \n<!-- image -->  \nTable 7. Goodness-of-Fit Metrics for Different Dependent Variables  \na The f esc , Ly α variable is not included in the fiducial model for the f esc , Ly α and L (Ly α ) predictions.  \nhas a greater uncertainty (26% on average for f esc , Ly α vs 21% for the Ly α flux alone). Again the L (Ly α ) JWST model (Table 7 and Figure 12b) has slightly worse R 2 and RMS metrics than the fiducial model because of the lack of information about line-of-sight gas.  \nTables 10 and 11 in the Appendix list the topranked variables for predicting f esc , Ly α and L (Ly α ). For f esc , Ly α , we find that H i absorption strength and nebular dust attenuation are two of the most important  \npredictive variables. In contrast, the UV dust attenuation is less important, always ranking after E(B-V) neb in forward and backward selection. As a nebular emission line, Ly α escape is more closely linked to the nebular attenuation rather than the attenuation experienced by the stellar light.  \nFor L (Ly α ), the most important variables include Σ SFR , M 1500 , E(B-V) neb , EW(H β ), and EW(H i ,abs). These five variables are also statistically signifi-  \n<!-- image -->  \n(a)  \n<!-- image -->  \n(b)  \nFigure 10. Examples of L (LyC) model predictions using the variable sets for the fiducial model (a) and the JWST model (b). We plot galaxies with predicted L (LyC) ∼ 0 at the bottom of each panel. Symbols are the same as in Figure 2. The models fail to accurately predict L (LyC).  \n<!-- image -->  \n<!-- image -->  \nFigure 11. Examples of f esc , Ly α model predictions using the variable sets for the fiducial model (a), the fiducial model with EW(H i ,abs) (b) and the JWST model (c). Symbols are the same as in Figure 2. We plot galaxies with negative f esc , Ly α at f esc , Ly α =2E-3. Of the three models, the fiducial model with EW(H i ,abs) predicts f esc , Ly α with the lowest RMS scatter.  \n<!-- image -->  \ncant in either the fiducial model (Σ SFR ), the fiducial+EW(H i ,abs) model (EW(H β ), EW(H i ,abs)), or both (E(B-V) neb and M 1500 ). The dependence on UV luminosity reflects the fact that higher SFRs enhance Ly α production, while the dependence on E(B-V) neb and EW(H i ,abs) arises because lower dust attenuation and lower H i column densities enhance Ly α escape. As a recombination line, the H β EW is also directly linked to nebular Ly α production. While the dependence on Σ SFR may reflect the role of feedback in Ly α escape (e.g. Kim et al. 2021), it may also serve to quantify aperture losses, as more extended galaxies may have more of their Ly α emission scatter outside the COS aperture.  \nRunnholm et al. (2020) also assess the most important variables in predicting L (Ly α ) using the Lymanα  \nReference Sample (LARS). They find that the top five ranked variables by both forward and backward selection are SFR, E(B-V) neb , M ∗ , r 50 , and the gas covering fraction. If using observables, the top variables are the FUV luminosity and the UV size. Despite the different methodologies between the LARS Ly α measurements and the LzLCS+ Ly α measurements, our results generally agree with their findings, as Σ SFR , E(B-V) neb , and the galaxy FUV luminosity also appear among our most important predictors.', '5. DISCUSSION': "Both theoretical and observational studies have long recognized that LyC escape is likely anisotropic due to the inherently inhomogeneous structure of the ISM (e.g., Gnedin et al. 2008; Wise & Cen 2009; Kim et al. 2013;  \n<!-- image -->  \n(a)  \nFigure 12. Examples of L (Ly α ) model predictions using the variable sets for the fiducial model (a) and the JWST model (b). Symbols are the same as in Figure 2. We exclude galaxies with negative L (Ly α ). The models do a better job at predicting L (Ly α ) compared to f esc , Ly α .  \n<!-- image -->  \nZastrow et al. 2013; Paardekooper et al. 2015; Cen & Kimm 2015; Trebitsch et al. 2017; Keenan et al. 2017; Rivera-Thorsen et al. 2017, 2019; Kim et al. 2023). LyC may escape in a 'picket-fence' geometry, emerging only through holes in the ISM (e.g., Heckman et al. 2001), and simulations find that supernova feedback may carve these low-density paths in the ISM (e.g., Wise & Cen 2009; Cen & Kimm 2015; Paardekooper et al. 2015; Trebitsch et al. 2017). Observationally, the detection of residual flux from saturated low-ionization absorption lines (e.g., Heckman et al. 2001, 2011; Gazagnes et al. 2018) and the co-existence of narrow Ly α emission with underlying broad absorption (e.g., McKinney et al. 2019) provide evidence for optical depth variations within galaxies. Highly ionized channels observed in some low-redshift galaxies (e.g., Zastrow et al. 2013; Bik et al. 2018) and spatially resolved LyC observations in the lensed Sunburst Arc at z = 2 . 37 (Rivera-Thorsen et al. 2019; Kim et al. 2023) likewise imply that LyC may escape through narrow opening angles, perhaps driven by outflows (e.g., Amor'ın et al. 2024). Even when f esc is large, with escape occurring in most or all directions, H i column density variations will still lead to variations in the observed f esc along different lines of sight (e.g., Gazagnes et al. 2020).  \nBecause of the anisotropy of LyC escape, estimating global f esc values from line-of-sight LyC detections is difficult, particularly since the observed f esc depends both on galaxy properties and on orientation. The multivariate Cox models present new evidence for the anisotropy of LyC escape by showing that the global properties of  \nthe strongest LCEs ( f esc ≳ 0 . 4) do not distinguish them from weaker LCEs ( f esc ∼ 0 . 04; see § 3.6). At z ∼ 3, Nakajima et al. (2020) reach a similar conclusion, finding that high O32 galaxies can have wide variations in f esc , despite sharing comparable spectral properties. We find that the only variable that successfully predicts the high observed f esc of the strongest LCEs is EW(H i ,abs), which traces the line-of-sight H i absorption. This result supports the picture described by Gazagnes et al. (2020), where galaxies with high O32 may have sufficient ionizing radiation to ionize most lines of sight, but where the highest f esc occurs in only the most diffuse channels.  \nMultivariate survival analysis models such as those presented here may offer a tool for quantifying the line-of-sight variation in f esc among galaxy populations. Given a large reference sample of galaxies with measurements of relevant global properties, the Cox models generate the probability distribution of f esc for a given set of variables. If the variation in the observed f esc for similar galaxies results only from orientation effects, the derived probability distributions should reflect the distribution of f esc across all sightlines.  \nAlthough the outsized importance of EW(H i ,abs) in accurately predicting high observed f esc suggests that favorable orientation accounts for the strongest LCEs in the LzLCS+, orientation alone may not fully explain the scatter between predicted and observed f esc across the entire LzLCS+ sample. LyC escape will depend on local properties, which global variables may not accurately capture (e.g., Rivera-Thorsen et al. 2019). Observational studies at higher spatial resolution and sim-  \nulations that track the evolution of f esc with high spatial resolution will help reveal the origin of f esc variations in otherwise similar galaxies.  \nThe variable rankings in § 3.3 show that the two most essential variables in predicting the line-of-sight f esc are variables sensitive to the line-of-sight H i gas (EW(H i ,abs)) and the line-of-sight dust (E(B-V) UV or β 1550 ). Although z > 6 observations will not be able to measure the former, β 1550 is easily accessible. JWST observations of the epoch of reionization are rapidly expanding the available measurements of β 1550 at high redshift (e.g., Schaerer et al. 2022a; Mascia et al. 2023; Saxena et al. 2024; Morishita et al. 2024), and this variable alone can provide useful constraints on f esc (Chisholm et al. 2022). By examining the spread of β 1550 values for a fixed intrinsic M 1500 or set of galaxy properties, high-redshift studies could gain insight into the possible variation in the line-of-sight f esc among a population, at least the variation due to changing levels of dust attenuation. To understand the H i component of the LyC absorption and its variation, high-redshift studies will need to rely on simulations and lower-redshift observations.  \nOur models provide some clues as to which physical properties most affect a galaxy's H i distribution and resulting f esc . In the absence of EW(H i ,abs) and Ly α measurements, which more directly trace H i , the models find that Σ SFR and O32 are the most important variables in predicting f esc , given constraints on the dust attenuation from β 1550 ( § 3.4). These variables suggest that feedback is key in determining a galaxy's average H i column density. Interestingly, we find in § 3.5 that O32 may be more important in lower mass galaxies (log 10 ( M ∗ / M ⊙ ) < 8 . 8), whereas Σ SFR takes on greater importance in the higher mass LzLCS+ galaxies. Previous studies have suggested that the mechanism responsible for LyC escape may be different in different galaxy populations (e.g., Flury et al. 2022b; Katz et al. 2023). Here, we find that differences in the multivariate f esc models are specifically associated with galaxy mass, rather than luminosity. We hypothesize that the typical star formation histories of the low- and high-mass LzLCS+ galaxies may differ, with the lower-mass starbursts exhibiting burstier star formation (e.g., Lee et al. 2007). Radiative feedback in extremely young, strong starbursts may drive LyC escape over wide opening angles in these low-mass galaxies, whereas the cumulative effect of supernovae in longer-lived concentrated starbursts could regulate LyC escape at higher masses (e.g., Bremer & Dayal 2023).  \nWhile the same physical processes, such as bursty star formation with radiative feedback or an extended star  \nformation episode with supernova feedback, may occur at high redshift, they may not occur in exactly the same type of galaxies as the LzLCS+. In particular, at a given stellar mass, a high-redshift galaxy could have a different gas mass, total mass, or star formation history compared to the same stellar mass galaxy at low redshift. Nevertheless, the distinct models for the lowvs. high-mass subsets in the LzLCS+ demonstrate that different parameters and physical mechanisms may regulate LyC escape in different types of galaxies. These differences could explain why some studies disagree regarding the properties of LCEs (e.g., Alexandroff et al. 2015; Rutkowski et al. 2017; Izotov et al. 2018b; Naidu et al. 2018). Depending on the sample probed, high O32 or high Σ SFR may or may not be a relevant predictor of f esc . Understanding the origin of the relationships between O32, Σ SFR , and f esc in low- vs. high-mass LzLCS+ galaxies will require more observational and theoretical studies examining the interplay between star formation history, geometry, and feedback in these types of galaxies.", '6. SUMMARY': "LyC escape depends on a range of galaxy properties, from gas and dust distributions in the ISM to galaxy star formation histories and the resulting stellar feedback. Because of this complexity, f esc correlates observationally with a variety of parameters, and correlations with individual variables show high scatter. Using the results of the LzLCS and other z ∼ 0 . 3 LyC observations from the literature, we have generated new multivariate models to predict f esc . We use the statistical technique of survival analysis, specifically the Cox proportional hazards model (Cox 1972), which accounts for data from both LyC detections and f esc upper limits from nondetections. With these models, we explore which combinations of variables best predict f esc and discuss the possible physical basis for these trends.  \nWe summarize our main findings below.  \n- 1. The f esc predictions from our fiducial model match the observed f esc of z ∼ 0 . 3 LCEs with an RMS of 0.36 dex, and our best-performing model, which includes H i absorption lines, reaches an RMS of 0.31 dex. These multivariate models offer a substantial improvement over correlations with single variables, which can span ≳ 2 dex in f esc for a given galaxy property (e.g., Wang et al. 2021; Flury et al. 2022b; Saldana-Lopez et al. 2022; Chisholm et al. 2022). ( §§ 3.1-3.2).\n- 2. Our best-performing model accurately predicts f esc even for galaxies that are outliers in a singlevariable selection. This model, which includes H i  \n- absorption line measurements, is the only model that reproduces the high observed f esc values of the strongest LCEs in the LzLCS+. Among the variables we consider, weak H i absorption is the only distinctive property of these strong LCEs, which suggests that these high f esc values may occur only along fortuitous sightlines. ( § 3.6).\n- 3. The most important variables in the models, as determined by forward and backward selection, are the EW of the H i Lyman series absorption lines and either E(B-V) UV or the UV slope β 1550 . These properties trace the line-of-sight H i absorption and UV dust attenuation, which determine the line-of-sight f esc that we measure (e.g., SaldanaLopez et al. 2022). In our fiducial model, E(BV) UV is a statistically significant variable, along with f esc , Ly α , Σ SFR , E(B-V) neb , and O32. Excluding E(B-V) UV or f esc , Ly α has the greatest adverse effects on model performance, and including f esc , Ly α appears to be particularly important for identifying non-LyC emitters. These analyses show which variables are the most important to include for accurate predictions; models using only the statistically significant variables or the 3-5 topranked variables achieve comparable accuracy to models with a longer list of inputs ( §§ 3.2-3.3).\n- 4. We generate a 'JWST model' for predicting f esc , which can apply to z > 6 observations. The model excludes Ly α and absorption line measurements, which will be inaccessible for most z > 6 galaxies. When compared with the z ∼ 0 . 3 LCEs, this model has a higher scatter (RMS= 0 . 47 dex) than the fiducial model, but it can distinguish LCEs from non-emitters. Three variables in this model are statistically significant: β 1550 , Σ SFR , and O32. These three variables alone can predict f esc with an RMS of 0.46 dex in the z ∼ 0 . 3 sample. ( § 3.4).\n- 5. We find that f esc may depend on different variables for high- vs. low-mass galaxies. Σ SFR is statistically significant for the high-mass subsample but not the low-mass subsample and is the top-ranked variable for the high-mass subsample by forward and backward selection. In contrast, M 1500 and O32 are statistically significant for the low-mass subsample but not for the high-mass subsample. These differences suggest that the physical causes of LyC escape and the feedback mechanisms involved could differ across the galaxy population, with ionization and radiative feedback potentially playing a greater role in low-mass galaxies. ( § 3.5).\n- 6. The multivariate Cox models are less successful at accurately predicting alternative measures of LyC escape, such as F λ LyC /F λ 1100 and the LyC luminosity. We hypothesize that accurately predicting these properties requires information about the unattenuated SED shape, which is not included in our set of variables. ( § 4.1).\n- 7. We also create Cox models to predict Ly α , specifically the f esc , Ly α and Ly α luminosity. Unlike the LyC f esc predictions, H i absorption EW and E(BV) neb , rather than E(B-V) UV , are two of the most important variables for predicting f esc , Ly α . For L (Ly α ), top-ranked and statistically significant variables are the H i absorption EW, E(B-V) neb , M 1500 , EW(H β ), and Σ SFR , which are similar to results found for the LARS sample (Runnholm et al. 2020). ( § 4.2).  \nOur models highlight the complexity of LyC escape and the dependence of f esc on multiple variables. Multivariate Cox models represent a new tool to investigate the relationships between f esc and galaxy properties and to predict f esc in the epoch of reionization. These models offer a substantial improvement over f esc predictions based on a single variable, and the models presented here can be refined and expanded to incorporate future observations. Our results show which variables have the greatest effect on f esc predictions for the LzLCS+ sample and for galaxies in different mass ranges. By exploring these same variables, future simulations could try to reproducing these scalings and determine their physical origin. These flexible multivariate models also offer an empirical approach to predicting f esc at high redshift. In a follow-up paper (Paper II, Jaskot et al., in prep.), we apply Cox models derived from the LzLCS+ to z ∼ 3 and z > 6 samples and provide the best-fit models to the community for application to future observations. 2  \n- 2 We have also created a python script, available at https://github. com/sflury/LyCsurv, that generates Cox model f esc predictions using any desired combination of variables from the LzLCS+. The LzLCS+ dataset used in this paper is available at the same page, and version 0.1.0 of the code is archived in Zenodo (Flury et al. 2024).  \nWe thank Kate Follette for interesting discussions that led to this project, and we thank the anonymous referee for comments that improved the clarity of this work. AEJ and SRF acknowledge support from NASA through grant HST-GO-15626. STScI is operated by AURA under NASA contract NAS-5-26555. ASL acknowledges support from Knut and Alice Wallenberg Foundation. MT acknowledges support from the NWO grant 0.16.VIDI.189.162 ('ODIN').", 'APPENDIX': "Tables 8-11 list the four top-ranked variables by forward selection, backward selection, and an MC rank ordering for the F λ LyC /F λ 1100 , L (LyC), f esc , Ly α , and L (Ly α ) predictions. We describe these ranking methods in § 2.4 and discuss the results in §§ 4.1 and 4.2. The 'Full Model' uses β 1550 as a measure of UV dust attenuation and includes multiple variables that trace Ly α or UV absorption line strength. The 'Fiducial+HI Model' uses E(B-V) UV to quantify dust attenuation, f esc , Ly α as a measure of Ly α , and EW(H i ,abs) as the only absorption line variable. The 'JWST Model' excludes absorption line and Ly α variables.", 'T able 8 . Most Imp ortan t V ariables for F λ LyC /F λ 1100 Predictions': "ariables  \nV  \ned  \nop-Rank  \nT  \nOrder  \nRank  \nMC  \ny  \nb  \nariables  \nV  \nop  \nT  \n]  \n[1.09  \n)  \nSFR  \n(Σ  \n10  \nlog  \n3]  \n[1.1  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[1.94]  \n,abs)  \ni  \n(H  \nEW  \n.00]  \n[1  \n,abs)  \ni  \nEW(H  \n[2.87]  \n,abs)  \ni  \nW(H  \nE  \n]  \n[1.01  \n,abs)  \ni  \nEW(H  \n2.14]  \n[  \n(O32)  \n10  \nlog  \n[3.22]  \n(O32)  \n10  \nlog  \n[3.64]  \nα  \nLy  \n,  \nesc  \nf  \n[3.57]  \nα  \nLy  \n,  \nesc  \nf  \n[5.83]  \n,abs)  \ni  \n(H  \nl  \nR  \n.86]  \n[2  \n,abs)  \ni  \n(H  \nl  \nR  \n[3.11]  \n))  \nβ  \n(EW(H  \n10  \nlog  \n[5.45]  \n))  \nβ  \n(EW(H  \n10  \nlog  \n[4.58]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[3.79]  \n1500  \nM  \n95]  \n[5.  \n))  \nβ  \n(EW(H  \n10  \nlog  \n[7.20]  \nα  \nLy  \n,  \nesc  \nf  \n[4.73]  \n1550  \nβ  \n[3.97]  \n1550  \nβ  \n.85]  \n[4  \n(O32)  \n10  \nlog  \n[4.21]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[6.68]  \n)  \nα  \nEW(Ly  \n[7.50]  \n)  \nSFR  \n(Σ  \n10  \nlog  \nor  \n)  \nα  \nEW(Ly  \nNote -Num b ers in brac k ets indicate the mean of the v ariable's ranks in the 100 MC runs.  \n.71]  \n[1  \n(O32)  \n10  \nlog  \n[1.17]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[2.95]  \n1500  \nM  \n[2.61]  \n1500  \nM  \n[3.05]  \n1500  \nM  \n[1.00]  \n)  \nα  \n(Ly  \nL  \n[2.87]  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[3.53]  \n1500  \nM  \n[4.27]  \nα  \nLy  \n,  \nesc  \nf  \n]  \n[3.30  \n)  \nSFR  \n(Σ  \n10  \nlog  \n6.08]  \n[  \n,abs)  \ni  \nEW(H  \n[6.21]  \nneb  \nE(B-V)  \n[4.58]  \n))  \nβ  \n(EW(H  \n10  \nlog  \n.38]  \n[4  \n(O32)  \n10  \nlog  \n]  \n[4.54  \n(O32)  \n10  \nlog  \n[3.32]  \nα  \nLy  \n,  \nesc  \nf  \n52]  \n[6.  \n(O32)  \n10  \nog  \nl  \n[6.96]  \n)  \nα  \nEW(Ly", 'T able 9 . Most Imp ortan t V ariables for L (LyC) Predi ctions': "l  \nde  \nMo  \nST  \nJW  \ndel  \nMo  \ncial+HI  \nFidu  \ndel  \nMo  \null  \nF  \nRank  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nBac  \nard  \norw  \nF  \nariables  \nV  \ned  \nop-Rank  \nT  \nOrder  \nRank  \nMC  \ny  \nb  \nariables  \nV  \nop  \nT  \n[2.55]  \n1500  \nM  \n[3.40]  \n1550  \nβ  \n9]  \n[3.3  \n,abs)  \ni  \nEW(H  \n[3.28]  \n,abs)  \ni  \nEW(H  \n[4.78]  \n)  \nα  \nEW(Ly  \n[4.43]  \n1500  \nM  \nNote -Num b ers in brac k ets indicate the mean of the v ariable's ranks in the 100 MC runs.", 'T able 10 . Most Imp ortan t V ariables for f esc , Ly α Predictions a': 'del  \nMo  \nJWST', 'Fiducial+HI Mo del': "del  \nMo  \null  \nF  \nRank  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nBac  \nard  \norw  \nF  \nard  \nkw  \nc  \na  \nB  \nard  \norw  \nF  \nT op-Rank ed V ariable s  \n))  \nβ  \n(EW(H  \n10  \nlog  \n)  \nSFR  \n(Σ  \n10  \nlog  \n))  \nβ  \n(EW(H  \n10  \nlog  \n)  \nSFR  \n(Σ  \n10  \nlog  \n1500  \nM  \n)  \nSFR  \n(Σ  \n10  \nlog  \n1  \n1500  \nM  \nneb  \nE(B-V)  \n1500  \nM  \nneb  \nE(B-V)  \n))  \nβ  \n(EW(H  \n10  \nlog  \nneb  \nE(B-V)  \n2  \n(O/H)  \n10  \n12+log  \n1500  \nM  \n,abs)  \ni  \nEW(H  \n1500  \nM  \n,abs)  \ni  \nEW(H  \n1500  \nM  \n3  \n)  \nSFR  \n(Σ  \n10  \nlog  \n(O/H)  \n10  \n12+log  \nUV  \nE(B-V)  \n(O/H)  \n10  \n12+log  \nneb  \nE(B-V)  \n(O/H)  \n10  \n12+log  \n4  \nOrder  \nRank  \nMC  \ny  \nb  \nariables  \nV  \nop  \nT  \n[1.94]  \n1500  \nM  \n0]  \n[1.0  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[1.98]  \n1500  \nM  \n0]  \n[2.1  \n)  \nSFR  \n(Σ  \n10  \nlog  \n[1.93]  \n1500  \nM  \n73]  \n[1.  \n)  \nSFR  \n(Σ  \n10  \nlog  \nruns.  \nMC  \n100  \nthe  \nin  \nranks  \nariable's  \nv  \nthe  \nof  \nmean  \nthe  \nindicate  \nets  \nk  \nbrac  \nin  \ners  \nb  \n-Num  \nNote  \na V ariable lists exclude an y measures of Ly α .", 'T able 11 . Most Imp ortan t V ariables for L (Ly α ) Predic tions a': 'del  \nMo  \nJWST', 'REFERENCES': 'Worseck, G., Prochaska, J. X., Hennawi, J. F., & McQuinn,  \nM. 2016, ApJ, 825, 144,  \ndoi: 10.3847/0004-637X/825/2/144  \nWyithe, J. S. B., & Loeb, A. 2013, MNRAS, 428, 2741,  \ndoi: 10.1093/mnras/sts242  \nXu, X., Henry, A., Heckman, T., et al. 2023, ApJ, 943, 94,  \ndoi: 10.3847/1538-4357/aca89a Yang, H., Malhotra, S., Gronke, M., et al. 2017, ApJ, 844, 171, doi: 10.3847/1538-4357/aa7d4d Zackrisson, E., Inoue, A. K., & Jensen, H. 2013, ApJ, 777,  \n39, doi: 10.1088/0004-637X/777/1/39  \nZastrow, J., Oey, M. S., Veilleux, S., & McDonald, M. 2013, ApJ, 779, 76, doi: 10.1088/0004-637X/779/1/76'}