docs: def lrt(data, model_full, model_reduced, df=None):
      """
      Compare two nested models using a likelihood ratio test
  
      Parameters
      ----------
      {0}
      model_full : obj
          A frozen scipy distribution object representing the full model
          (more complex model).
      model_reduced : scipy distribution object
          A frozen scipy distribution object representing the reduced model
          (simpler model).
      df : int
          The degrees of freedom for the lrt (optional). If none, df is
          calculated as the difference between the number of parameters in the
          full and reduced models.
  
      Returns
      -------
      tuple
          G^2 statistic, p-value
  
      Notes
      -----
      Parameters of distribution objects must be given as keyword arguments. Ex.
      ``norm = stats.norm(loc=0, scale=1)``
  
      A p-value < alpha suggests significant evidence for the full (more complex)
      model. In other words, the null hypothesis is that the reduced model is
      correct
  
      The LRT only applies to nested models. The G^2 statistic and G-test rely on
      the assumption that the test statistic is approximately chi-squared
      distributed. This assumption breaks down for small samples sizes.
  
      Examples
      --------
  
      >>> import macroeco.models as md
      >>> import macroeco.compare as comp
  
      >>> # Generate random data
      >>> rand_samp = md.nbinom_ztrunc.rvs(20, 0.5, size=100)
  
      >>> # Fit Zero-truncated NBD (Full model)
      >>> mle_nbd = md.nbinom_ztrunc.fit_mle(rand_samp)
  
      >>> # Fit a logseries (limiting case of Zero-truncated NBD, reduced model)
      >>> mle_logser = md.logser.fit_mle(rand_samp)
  
      >>> # Compare models with LRT
      >>> comp.lrt(rand_samp, md.nbinom_ztrunc(mu=mle_nbd[0], k_agg=mle_nbd[1]), md.logser(p=mle_logser[0]))
      (15.33429080890221, 9.0066719644695982e-05)
  
      >>> # Reject the null hypothesis that the logseries is a better model and
      >>> # choose the Zero-truncated NBD.
  
      """
  
      # Calculate G^2 statistic
      ll_full = nll(data, model_full) * -1
      ll_reduced = nll(data, model_reduced) * -1
      test_stat = 2 * (ll_full - ll_reduced)
  
      # Set df if necessary
      if not df:
          df =  ( len(model_full.args) + len(model_full.kwds)
                - len(model_reduced.args) - len(model_reduced.kwds) )
  
      return test_stat, stats.chisqprob(test_stat, df)