'Use Lagrange Multiplier test to test whether restricted model is correct
Parameters
restricted : Result instance
The restricted model is assumed to be nested in the
current model. The result instance of the restricted model
is required to have two attributes, residual sum of
squares, `ssr`, residual degrees of freedom, `df_resid`.
demean : bool
Flag indicating whether the demean the scores based on the
residuals from the restricted model.  If True, the
covariance of the scores are used and the LM test is
identical to the large sample version of the LR test.
Returns
lm_value : float
test statistic, chi2 distributed
p_value : float
p-value of the test statistic
df_diff : int
degrees of freedom of the restriction, i.e. difference in df
between models
Notes
TODO: explain LM text'

def compare_lm_test(self, restricted, demean=True, use_lr=False):

    import statsmodels.stats.sandwich_covariance as sw
    from numpy.linalg import inv
    if (not self._is_nested(restricted)):
        raise ValueError('Restricted    model    is    not    nested    by    full    model.')
    wresid = restricted.wresid
    wexog = self.model.wexog
    scores = (wexog * wresid[:, None])
    n = self.nobs
    df_full = self.df_resid
    df_restr = restricted.df_resid
    df_diff = (df_restr - df_full)
    s = scores.mean(axis=0)
    if use_lr:
        scores = (wexog * self.wresid[:, None])
        demean = False
    if demean:
        scores = (scores - scores.mean(0)[None, :])
    cov_type = getattr(self, 'cov_type', 'nonrobust')
    if (cov_type == 'nonrobust'):
        sigma2 = np.mean((wresid ** 2))
        XpX = (np.dot(wexog.T, wexog) / n)
        Sinv = inv((sigma2 * XpX))
    elif (cov_type in ('HC0', 'HC1', 'HC2', 'HC3')):
        Sinv = inv((np.dot(scores.T, scores) / n))
    elif (cov_type == 'HAC'):
        maxlags = self.cov_kwds['maxlags']
        Sinv = inv((sw.S_hac_simple(scores, maxlags) / n))
    elif (cov_type == 'cluster'):
        groups = self.cov_kwds['groups']
        Sinv = inv(sw.S_crosssection(scores, groups))
    else:
        raise ValueError(('Only    nonrobust,    HC,    HAC    and    cluster    are    ' + 'currently    connected'))
    lm_value = (n * chain_dot(s, Sinv, s.T))
    p_value = stats.chi2.sf(lm_value, df_diff)
    return (lm_value, p_value, df_diff)
