code: def _fit_firth(
      beta_init: NDArray[(Any,), Float],
      X: NDArray[(Any, Any), Float],
      y: NDArray[(Any,), Float],
      offset: NDArray[(Any,), Float],
      convergence_limit: float = 1e-5,
      deviance_tolerance: float = 1e-6,
      max_iter: int = 250,
      max_step_size: int = 5,
      max_half_steps: int = 25,
  ) -> Optional[FirthFit]:
      """
      Firth’s bias-reduced penalized-likelihood logistic regression, based on the regenie implementation:
      https://www.nature.com/articles/s41588-021-00870-7
  
      Args:
          beta_init : Initial beta values
          X : Independent variable (covariate for null fit, genotype for SNP fit)
          y : Dependent variable (phenotype)
          offset : Phenotype offset only for null fit, phenotype + covariate effects for SNP fit
          convergence_limit : Convergence is reached if all entries of the penalized score have smaller magnitude
          deviance_tolerance : Non-inferiority margin when halving step size (default from regenie)
          max_iter : Maximum number of Firth iterations (default from regenie)
          max_step_size : Maximum step size during a Firth iteration (default from regenie)
          max_half_steps : Maximum number of half-steps during a Firth iteration (default from regenie)
  
      Returns:
          None if the fit failed. Otherwise, a FirthFit object containing the fit information.
      """
  
      n_iter = 0
      beta = beta_init.copy()
      model = Model(X, y, offset)
      log_likelihood = _calculate_log_likelihood(beta, model)
  
      while n_iter < max_iter:
          invI = np.linalg.pinv(log_likelihood.I)
  
          # build hat matrix
          rootG = np.sqrt(log_likelihood.pi * (1 - log_likelihood.pi))
          rootG_X = (
              rootG[:, None] * X
          )  # equivalent to sqrt(diagflat(pi * (1 - pi))) @ X
          h = np.sum((rootG_X @ invI) * rootG_X, axis=1)
  
          # modified score function
          U = X.T @ (
              y - log_likelihood.pi + h * (0.5 - log_likelihood.pi)
          )
  
          # f' / f''
          delta = invI @ U
  
          # force absolute step size to be less than max_step_size for each entry of beta
          mx = np.amax(np.abs(delta)) / max_step_size
          if mx > 1:
              delta /= mx
  
          # if the penalized log likelihood decreased, recompute with step-halving
          n_half_steps = 0
          while n_half_steps < max_half_steps:
              new_log_likelihood = _calculate_log_likelihood(
                  beta + delta, model
              )
              if (
                  new_log_likelihood.deviance
                  < log_likelihood.deviance + deviance_tolerance
              ):
                  break
              delta /= 2
              n_half_steps += 1
  
          beta += delta
          log_likelihood = new_log_likelihood
  
          if np.amax(np.abs(U)) < convergence_limit:
              break
  
          n_iter += 1
  
      if n_iter == max_iter:
          # Failed to converge
          return None
  
      return FirthFit(beta, log_likelihood)
  
file_path: glow/python/glow/gwas/approx_firth.py
imports: from typing import Any, Optional
  import numpy
  from dataclasses import dataclass
  import scipy
  from typeguard import typechecked
  from nptyping import Float, NDArray
  from scipy import stats
label_desc: Use image node_id 2 for calling a global function with example usage: _fit_firth(beta_init, X, y, offset, convergence_limit, deviance_tolerance, max_iter, max_step_size, max_half_steps) and returns: FirthFit, None