'Returns the sampling covariance matrix of the regression
parameters and related quantities.
Returns
cov_robust : array-like
The robust, or sandwich estimate of the covariance, which
is meaningful even if the working covariance structure is
incorrectly specified.
cov_naive : array-like
The model-based estimate of the covariance, which is
meaningful if the covariance structure is correctly
specified.
cmat : array-like
The center matrix of the sandwich expression, used in
obtaining score test results.'

def _covmat(self):

    endog = self.endog_li
    exog = self.exog_li
    varfunc = self.family.variance
    cached_means = self.cached_means
    (bmat, cmat) = (0, 0)
    for i in range(self.num_group):
        (expval, lpr) = cached_means[i]
        resid = (endog[i] - expval)
        dmat = self.mean_deriv(exog[i], lpr)
        sdev = np.sqrt(varfunc(expval))
        rslt = self.cov_struct.covariance_matrix_solve(expval, i, sdev, (dmat, resid))
        if (rslt is None):
            return (None, None, None, None)
        (vinv_d, vinv_resid) = tuple(rslt)
        f = (self.weights_li[i] if (self.weights is not None) else 1.0)
        bmat += (f * np.dot(dmat.T, vinv_d))
        dvinv_resid = (f * np.dot(dmat.T, vinv_resid))
        cmat += np.outer(dvinv_resid, dvinv_resid)
    scale = self.estimate_scale()
    bmati = np.linalg.inv(bmat)
    cov_naive = (bmati * scale)
    cov_robust = np.dot(bmati, np.dot(cmat, bmati))
    cov_naive *= self.scaling_factor
    cov_robust *= self.scaling_factor
    return (cov_robust, cov_naive, cmat)
