class_docstr: Linear Mixed Effects Model
  
  Parameters
  ----------
  endog : 1d array_like
      The dependent variable
  exog : 2d array_like
      A matrix of covariates used to determine the
      mean structure (the "fixed effects" covariates).
  groups : 1d array_like
      A vector of labels determining the groups -- data from
      different groups are independent
  exog_re : 2d array_like
      A matrix of covariates used to determine the variance and
      covariance structure (the "random effects" covariates).  If
      None, defaults to a random intercept for each group.
  exog_vc : VCSpec instance or dict-like (deprecated)
      A VCSPec instance defines the structure of the variance
      components in the model.  Alternatively, see notes below
      for a dictionary-based format.  The dictionary format is
      deprecated and may be removed at some point in the future.
  use_sqrt : bool
      If True, optimization is carried out using the lower
      triangle of the square root of the random effects
      covariance matrix, otherwise it is carried out using the
      lower triangle of the random effects covariance matrix.
  missing : str
      The approach to missing data handling
  
  Notes
  -----
  If `exog_vc` is not a `VCSpec` instance, then it must be a
  dictionary of dictionaries.  Specifically, `exog_vc[a][g]` is a
  matrix whose columns are linearly combined using independent
  random coefficients.  This random term then contributes to the
  variance structure of the data for group `g`.  The random
  coefficients all have mean zero, and have the same variance.  The
  matrix must be `m x k`, where `m` is the number of observations in
  group `g`.  The number of columns may differ among the top-level
  groups.
  
  The covariates in `exog`, `exog_re` and `exog_vc` may (but need
  not) partially or wholly overlap.
  
  `use_sqrt` should almost always be set to True.  The main use case
  for use_sqrt=False is when complicated patterns of fixed values in
  the covariance structure are set (using the `free` argument to
  `fit`) that cannot be expressed in terms of the Cholesky factor L.
  
  Examples
  --------
  A basic mixed model with fixed effects for the columns of
  ``exog`` and a random intercept for each distinct value of
  ``group``:
  
  >>> model = sm.MixedLM(endog, exog, groups)
  >>> result = model.fit()
  
  A mixed model with fixed effects for the columns of ``exog`` and
  correlated random coefficients for the columns of ``exog_re``:
  
  >>> model = sm.MixedLM(endog, exog, groups, exog_re=exog_re)
  >>> result = model.fit()
  
  A mixed model with fixed effects for the columns of ``exog`` and
  independent random coefficients for the columns of ``exog_re``:
  
  >>> free = MixedLMParams.from_components(
                   fe_params=np.ones(exog.shape[1]),
                   cov_re=np.eye(exog_re.shape[1]))
  >>> model = sm.MixedLM(endog, exog, groups, exog_re=exog_re)
  >>> result = model.fit(free=free)
  
  A different way to specify independent random coefficients for the
  columns of ``exog_re``.  In this example ``groups`` must be a
  Pandas Series with compatible indexing with ``exog_re``, and
  ``exog_re`` has two columns.
  
  >>> g = pd.groupby(groups, by=groups).groups
  >>> vc = {}
  >>> vc['1'] = {k : exog_re.loc[g[k], 0] for k in g}
  >>> vc['2'] = {k : exog_re.loc[g[k], 1] for k in g}
  >>> model = sm.MixedLM(endog, exog, groups, vcomp=vc)
  >>> result = model.fit()
class_docstr_tok: ["Linear","Mixed","Effects","Model","Parameters","--","--","--","--","--","endog",":","1d","array_like","The","dependent","variable","exog",":","2d","array_like","A","matrix","of","covariates","used","to","determine","the","mean","structure","(","the","``","fixed","effects","''","covariates",")",".","groups",":","1d","array_like","A","vector","of","labels","determining","the","groups","--","data","from","different","groups","are","independent","exog_re",":","2d","array_like","A","matrix","of","covariates","used","to","determine","the","variance","and","covariance","structure","(","the","``","random","effects","''","covariates",")",".","If","None",",","defaults","to","a","random","intercept","for","each","group",".","exog_vc",":","VCSpec","instance","or","dict-like","(","deprecated",")","A","VCSPec","instance","defines","the","structure","of","the","variance","components","in","the","model",".","Alternatively",",","see","notes","below","for","a","dictionary-based","format",".","The","dictionary","format","is","deprecated","and","may","be","removed","at","some","point","in","the","future",".","use_sqrt",":","bool","If","True",",","optimization","is","carried","out","using","the","lower","triangle","of","the","square","root","of","the","random","effects","covariance","matrix",",","otherwise","it","is","carried","out","using","the","lower","triangle","of","the","random","effects","covariance","matrix",".","missing",":","str","The","approach","to","missing","data","handling","Notes","--","--","-","If","`","exog_vc","`","is","not","a","`","VCSpec","`","instance",",","then","it","must","be","a","dictionary","of","dictionaries",".","Specifically",",","`","exog_vc","[","a","]","[","g","]","`","is","a","matrix","whose","columns","are","linearly","combined","using","independent","random","coefficients",".","This","random","term","then","contributes","to","the","variance","structure","of","the","data","for","group","`","g","`",".","The","random","coefficients","all","have","mean","zero",",","and","have","the","same","variance",".","The","matrix","must","be","`","m","x","k","`",",","where","`","m","`","is","the","number","of","observations","in","group","`","g","`",".","The","number","of","columns","may","differ","among","the","top-level","groups",".","The","covariates","in","`","exog","`",",","`","exog_re","`","and","`","exog_vc","`","may","(","but","need","not",")","partially","or","wholly","overlap",".","`","use_sqrt","`","should","almost","always","be","set","to","True",".","The","main","use","case","for","use_sqrt=False","is","when","complicated","patterns","of","fixed","values","in","the","covariance","structure","are","set","(","using","the","`","free","`","argument","to","`","fit","`",")","that","can","not","be","expressed","in","terms","of","the","Cholesky","factor","L.","Examples","--","--","--","--","A","basic","mixed","model","with","fixed","effects","for","the","columns","of","``","exog","``","and","a","random","intercept","for","each","distinct","value","of","``","group","``",":",">",">",">","model","=","sm.MixedLM","(","endog",",","exog",",","groups",")",">",">",">","result","=","model.fit","(",")","A","mixed","model","with","fixed","effects","for","the","columns","of","``","exog","``","and","correlated","random","coefficients","for","the","columns","of","``","exog_re","``",":",">",">",">","model","=","sm.MixedLM","(","endog",",","exog",",","groups",",","exog_re=exog_re",")",">",">",">","result","=","model.fit","(",")","A","mixed","model","with","fixed","effects","for","the","columns","of","``","exog","``","and","independent","random","coefficients","for","the","columns","of","``","exog_re","``",":",">",">",">","free","=","MixedLMParams.from_components","(","fe_params=np.ones","(","exog.shape","[","1","]",")",",","cov_re=np.eye","(","exog_re.shape","[","1","]",")",")",">",">",">","model","=","sm.MixedLM","(","endog",",","exog",",","groups",",","exog_re=exog_re",")",">",">",">","result","=","model.fit","(","free=free",")","A","different","way","to","specify","independent","random","coefficients","for","the","columns","of","``","exog_re","``",".","In","this","example","``","groups","``","must","be","a","Pandas","Series","with","compatible","indexing","with","``","exog_re","``",",","and","``","exog_re","``","has","two","columns",".",">",">",">","g","=","pd.groupby","(","groups",",","by=groups",")",".groups",">",">",">","vc","=","{","}",">",">",">","vc","[","'","1","'","]","=","{","k",":","exog_re.loc","[","g","[","k","]",",","0","]","for","k","in","g","}",">",">",">","vc","[","'","2","'","]","=","{","k",":","exog_re.loc","[","g","[","k","]",",","1","]","for","k","in","g","}",">",">",">","model","=","sm.MixedLM","(","endog",",","exog",",","groups",",","vcomp=vc",")",">",">",">","result","=","model.fit","(",")"]
docstr: Returns the model's Hessian matrix.
  
  Calculates the Hessian matrix for the linear mixed effects
  model with respect to the parameterization in which the
  covariance matrix is represented directly (without square-root
  transformation).
  
  Parameters
  ----------
  params : MixedLMParams or array_like
      The model parameters at which the Hessian is calculated.
      If array-like, must contain the packed parameters in a
      form that is compatible with this model instance.
  
  Returns
  -------
  hess : 2d ndarray
      The Hessian matrix, evaluated at `params`.
  sing : boolean
      If True, the covariance matrix is singular and a
      pseudo-inverse is returned.
docstr_tok: ["Returns","the","model","'s","Hessian","matrix",".","Calculates","the","Hessian","matrix","for","the","linear","mixed","effects","model","with","respect","to","the","parameterization","in","which","the","covariance","matrix","is","represented","directly","(","without","square-root","transformation",")",".","Parameters","--","--","--","--","--","params",":","MixedLMParams","or","array_like","The","model","parameters","at","which","the","Hessian","is","calculated",".","If","array-like",",","must","contain","the","packed","parameters","in","a","form","that","is","compatible","with","this","model","instance",".","Returns","--","--","--","-","hess",":","2d","ndarray","The","Hessian","matrix",",","evaluated","at","`","params","`",".","sing",":","boolean","If","True",",","the","covariance","matrix","is","singular","and","a","pseudo-inverse","is","returned","."]
code: def hessian(self, params):
      """
      Returns the model's Hessian matrix.
  
      Calculates the Hessian matrix for the linear mixed effects
      model with respect to the parameterization in which the
      covariance matrix is represented directly (without square-root
      transformation).
  
      Parameters
      ----------
      params : MixedLMParams or array_like
          The model parameters at which the Hessian is calculated.
          If array-like, must contain the packed parameters in a
          form that is compatible with this model instance.
  
      Returns
      -------
      hess : 2d ndarray
          The Hessian matrix, evaluated at `params`.
      sing : boolean
          If True, the covariance matrix is singular and a
          pseudo-inverse is returned.
      """
  
      if type(params) is not MixedLMParams:
          params = MixedLMParams.from_packed(
              params,
              self.k_fe,
              self.k_re,
              use_sqrt=self.use_sqrt,
              has_fe=True,
          )
  
      fe_params = params.fe_params
      vcomp = params.vcomp
      cov_re = params.cov_re
      sing = False
  
      if self.k_re > 0:
          try:
              cov_re_inv = np.linalg.inv(cov_re)
          except np.linalg.LinAlgError:
              cov_re_inv = np.linalg.pinv(cov_re)
              sing = True
      else:
          cov_re_inv = np.empty((0, 0))
  
      # Blocks for the fixed and random effects parameters.
      hess_fe = 0.0
      hess_re = np.zeros(
          (self.k_re2 + self.k_vc, self.k_re2 + self.k_vc)
      )
      hess_fere = np.zeros((self.k_re2 + self.k_vc, self.k_fe))
  
      fac = self.n_totobs
      if self.reml:
          fac -= self.exog.shape[1]
  
      rvir = 0.0
      xtvix = 0.0
      xtax = [
          0.0,
      ] * (self.k_re2 + self.k_vc)
      m = self.k_re2 + self.k_vc
      B = np.zeros(m)
      D = np.zeros((m, m))
      F = [[0.0] * m for k in range(m)]
      for group_ix, group in enumerate(self.group_labels):
          vc_var = self._expand_vcomp(vcomp, group_ix)
          vc_vari = np.zeros_like(vc_var)
          ii = np.flatnonzero(vc_var >= 1e-10)
          if len(ii) > 0:
              vc_vari[ii] = 1 / vc_var[ii]
          if len(ii) < len(vc_var):
              sing = True
  
          exog = self.exog_li[group_ix]
          ex_r, ex2_r = self._aex_r[group_ix], self._aex_r2[group_ix]
          solver = _smw_solver(1.0, ex_r, ex2_r, cov_re_inv, vc_vari)
  
          # The residuals
          resid = self.endog_li[group_ix]
          if self.k_fe > 0:
              expval = np.dot(exog, fe_params)
              resid = resid - expval
  
          viexog = solver(exog)
          xtvix += np.dot(exog.T, viexog)
          vir = solver(resid)
          rvir += np.dot(resid, vir)
  
          for jj1, matl1, matr1, vsl1, vsr1, sym1 in self._gen_dV_dPar(
              ex_r, solver, group_ix
          ):
              ul = _dot(viexog.T, matl1)
              ur = _dot(matr1.T, vir)
              hess_fere[jj1, :] += np.dot(ul, ur)
              if not sym1:
                  ul = _dot(viexog.T, matr1)
                  ur = _dot(matl1.T, vir)
                  hess_fere[jj1, :] += np.dot(ul, ur)
  
              if self.reml:
                  ul = _dot(viexog.T, matl1)
                  ur = ul if sym1 else np.dot(viexog.T, matr1)
                  ulr = _dot(ul, ur.T)
                  xtax[jj1] += ulr
                  if not sym1:
                      xtax[jj1] += ulr.T
  
              ul = _dot(vir, matl1)
              ur = ul if sym1 else _dot(vir, matr1)
              B[jj1] += np.dot(ul, ur) * (1 if sym1 else 2)
  
              # V^{-1} * dV/d_theta
              E = [(vsl1, matr1)]
              if not sym1:
                  E.append((vsr1, matl1))
  
              for (
                  jj2,
                  matl2,
                  matr2,
                  vsl2,
                  vsr2,
                  sym2,
              ) in self._gen_dV_dPar(ex_r, solver, group_ix, jj1):
                  re = sum(
                      [
                          _multi_dot_three(matr2.T, x[0], x[1].T)
                          for x in E
                      ]
                  )
                  vt = 2 * _dot(
                      _multi_dot_three(vir[None, :], matl2, re),
                      vir[:, None],
                  )
  
                  if not sym2:
                      le = sum(
                          [
                              _multi_dot_three(matl2.T, x[0], x[1].T)
                              for x in E
                          ]
                      )
                      vt += 2 * _dot(
                          _multi_dot_three(vir[None, :], matr2, le),
                          vir[:, None],
                      )
  
                  D[jj1, jj2] += np.squeeze(vt)
                  if jj1 != jj2:
                      D[jj2, jj1] += np.squeeze(vt)
  
                  rt = _dotsum(vsl2, re.T) / 2
                  if not sym2:
                      rt += _dotsum(vsr2, le.T) / 2
  
                  hess_re[jj1, jj2] += rt
                  if jj1 != jj2:
                      hess_re[jj2, jj1] += rt
  
                  if self.reml:
                      ev = sum(
                          [_dot(x[0], _dot(x[1].T, viexog)) for x in E]
                      )
                      u1 = _dot(viexog.T, matl2)
                      u2 = _dot(matr2.T, ev)
                      um = np.dot(u1, u2)
                      F[jj1][jj2] += um + um.T
                      if not sym2:
                          u1 = np.dot(viexog.T, matr2)
                          u2 = np.dot(matl2.T, ev)
                          um = np.dot(u1, u2)
                          F[jj1][jj2] += um + um.T
  
      hess_fe -= fac * xtvix / rvir
      hess_re = hess_re - 0.5 * fac * (
          D / rvir - np.outer(B, B) / rvir**2
      )
      hess_fere = -fac * hess_fere / rvir
  
      if self.reml:
          QL = [np.linalg.solve(xtvix, x) for x in xtax]
          for j1 in range(self.k_re2 + self.k_vc):
              for j2 in range(j1 + 1):
                  a = _dotsum(QL[j1].T, QL[j2])
                  a -= np.trace(np.linalg.solve(xtvix, F[j1][j2]))
                  a *= 0.5
                  hess_re[j1, j2] += a
                  if j1 > j2:
                      hess_re[j2, j1] += a
  
      # Put the blocks together to get the Hessian.
      m = self.k_fe + self.k_re2 + self.k_vc
      hess = np.zeros((m, m))
      hess[0 : self.k_fe, 0 : self.k_fe] = hess_fe
      hess[0 : self.k_fe, self.k_fe :] = hess_fere.T
      hess[self.k_fe :, 0 : self.k_fe] = hess_fere
      hess[self.k_fe :, self.k_fe :] = hess_re
  
      return hess, sing
  
imports: import warnings
  import numpy
  import pandas
  import patsy
  from scipy import sparse
  from scipy.stats.distributions import norm
  from statsmodels.base._penalties import Penalty
  import statsmodels.base.model
  from statsmodels.tools import data
  from statsmodels.tools.decorators import cache_readonly
  from statsmodels.tools.sm_exceptions import ConvergenceWarning
label_desc: Use image node_id 16 for calling the MixedLM obj's underlying member method code with example usage: obj.hessian(params) and returns: hess,  sing