Small group size inference for logistic regression Please note this analysis is being done retrospectively (i.e., all data have been collected).
I have a binary response variable $y_i \in \{0, 1\}$ and two covariates $x_{1i}$ and $x_{2i}$. $x_{1i}$ is meant to be a control, and $x_{2i}$ is the covariate with which we'd like to perform inference against the $y_i$. Both $x_{1i}$ and $x_{2i}$ are binary explanatory variables.
With each combination of binary indicators $(x_{1i}, x_{2i}, y_i)$, counts range from 4 to about 15,000. In particular, when $y_i = 1$, the counts get extremely small (all less than 50). When $y_i = 0$, counts are all above 1,000.
My goal is to perform inference to get an idea of the association between $x_{2i}$ and $y_i$ when controlling for $x_{1i}$, and not prediction. How should one handle such a situation? Out-of-the-box logistic regression, I suspect, would not work for this case due to the small group sizes.
I would strongly prefer citations to literature that point to how to deal with such a situation, as this is for a paper.

As discussed in the comments, the main problem with inference from GLMs (such as logistic regression) with small sample sizes (at least in some categories) is the asymptotic nature of Wilks' theorem (which underlies p-values based on likelihood ratio tests and confidence intervals based on likelihood profiles). (This asymptotic constraint is much weaker than that used in Wald tests, which rely on the asymptotic normality of the sampling distribution of the parameters rather than the asymptotic $\chi^2$ distribution of differences in deviances.)
There is a literature on finite-size corrections under the rubric of Bartlett corrections, but the methods are often model-specific and are not widely implemented. Das et al. 2018 review a series of small-sample corrections; they do not provide software ("SAS code available upon request"). One of the options they explore is Firth correction, which is implemented in several R packages (logistf, brgml2) [I know Firth correction provides bias-corrected point estimates, I'm not sure how well it extends to inference.]
Another option is exact logistic regression; see Zamar et al. 2007 (their elrm package, which implements an MCMC computation, is archived on CRAN but can be installed via remotes::install_version("elrm", "1.2.4")). They refer to Hirji (2005) for theoretical details.

Das, Ujjwal, Subhra Sankar Dhar, and Vivek Pradhan. “Corrected Likelihood-Ratio Tests in Logistic Regression Using Small-Sample Data.” Communications in Statistics - Theory and Methods 47, no. 17 (September 2, 2018): 4272–85. https://doi.org/10.1080/03610926.2017.1373815.
Hirji, Karim F. Exact Analysis of Discrete Data. 1st edition. Boca Raton: Chapman and Hall/CRC, 2005.
Zamar, David, Brad McNeney, and Jinko Graham. “Elrm: Software Implementing Exact-Like Inference for Logistic Regression Models.” Journal of Statistical Software 21, no. 1 (October 8, 2007): 1–18. https://doi.org/10.18637/jss.v021.i03.
