code: def lm_robust_subset(
      score, k_constraints, score_deriv_inv, cov_score
  ):
      """general formula for score/LM test
  
      generalized score or lagrange multiplier test for constraints on a subset
      of parameters
  
      `params_1 = value`, where params_1 is a subset of the unconstrained
      parameter vector.
  
      It is assumed that all arrays are evaluated at the constrained estimates.
  
      Parameters
      ----------
      score : ndarray, 1-D
          derivative of objective function at estimated parameters
          of constrained model
      k_constraint : int
          number of constraints
      score_deriv_inv : ndarray, symmetric, square
          inverse of second derivative of objective function
          TODO: could be OPG or any other estimator if information matrix
          equality holds
      cov_score B :  ndarray, symmetric, square
          covariance matrix of the score. This is the inner part of a sandwich
          estimator.
      not cov_params V :  ndarray, symmetric, square
          covariance of full parameter vector evaluated at constrained parameter
          estimate. This can be specified instead of cov_score B.
  
      Returns
      -------
      lm_stat : float
          score/lagrange multiplier statistic
      p-value : float
          p-value of the LM test based on chisquare distribution
  
      Notes
      -----
      The implementation is based on Boos 1992 section 4.1. The same derivation
      is also in other articles and in text books.
  
      """
  
      # Notation in Boos
      # score `S = sum (s_i)
      # score_obs `s_i`
      # score_deriv `I` is derivative of score (hessian)
      # `D` is covariance matrix of score, OPG product given independent observations
  
      # k_params = len(score)
  
      # Note: I reverse order between constraint and unconstrained compared to Boos
  
      # submatrices of score_deriv/hessian
      # these are I22 and I12 in Boos
      # h_uu = score_deriv[-k_constraints:, -k_constraints:]
      h_uu = score_deriv_inv[:-k_constraints, :-k_constraints]
      h_cu = score_deriv_inv[-k_constraints:, :-k_constraints]
  
      # TODO: pinv or solve ?
      tmp_proj = h_cu.dot(np.linalg.inv(h_uu))
      tmp = np.column_stack(
          (-tmp_proj, np.eye(k_constraints))
      )  # , tmp_proj))
  
      cov_score_constraints = tmp.dot(cov_score.dot(tmp.T))
  
      # lm_stat2 = wscore.dot(np.linalg.pinv(inner).dot(wscore))
      # Let's assume inner is invertible, TODO: check if usecase for pinv exists
      lm_stat = score.dot(np.linalg.solve(cov_score_constraints, score))
      pval = stats.chi2.sf(lm_stat, k_constraints)
  
      #     # check second calculation Boos referencing Kent 1982 and Engle 1984
      #     # we can use this when robust_cov_params of full model is available
      #     #h_inv = np.linalg.inv(score_deriv)
      #     hinv = score_deriv_inv
      #     v = h_inv.dot(cov_score.dot(h_inv)) # this is robust cov_params
      #     v_cc = v[:k_constraints, :k_constraints]
      #     h_cc = score_deriv[:k_constraints, :k_constraints]
      #     # brute force calculation:
      #     h_resid_cu = h_cc - h_cu.dot(np.linalg.solve(h_uu, h_cu))
      #     cov_s_c = h_resid_cu.dot(v_cc.dot(h_resid_cu))
      #     diff = np.max(np.abs(cov_s_c - cov_score_constraints))
      return lm_stat, pval
  
imports: import warnings
  import numpy
  from scipy import stats
  from statsmodels.tools.decorators import cache_readonly
  from statsmodels.regression.linear_model import OLS
label_desc: Use image node_id 6 for calling a global function with example usage: lm_robust_subset(score, k_constraints, score_deriv_inv, cov_score) and returns: lm_stat,  pval