def _covmat(self):
        """
        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.
        cov_robust_bc : array-like
           The "bias corrected" robust covariance of Mancl and
           DeRouen.
        cmat : array-like
           The center matrix of the sandwich expression, used in
           obtaining score test results.
        """

        endog = self.endog_li
        exog = self.exog_li

        varfunc = self.family.variance
        cached_means = self.cached_means

        # Calculate the naive (model-based) and robust (sandwich)
        # covariances.
        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.

            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))

        # Calculate the bias-corrected sandwich estimate of Mancl and
        # DeRouen (requires cov_naive so cannot be calculated
        # in the previous loop).
        bcm = 0
        for i in range(self.num_group):

            expval, lpr = cached_means[i]
            resid = expval - endog[i]
            dmat = self.mean_deriv(exog[i], lpr)
            sdev = np.sqrt(varfunc(expval))

            rslt = self.cov_struct.covariance_matrix_solve(expval,
                                                  i, sdev, (dmat,))
            if rslt is None:
                return None, None, None, None
            vinv_d = rslt[0]
            vinv_d /= scale

            hmat = np.dot(vinv_d, cov_naive)
            hmat = np.dot(hmat, dmat.T).T

            f = self.weights_li[i] if self.weights is not None else 1.

            aresid = np.linalg.solve(np.eye(len(resid)) - hmat, resid)
            rslt = self.cov_struct.covariance_matrix_solve(expval, i,
                                                     sdev, (aresid,))
            if rslt is None:
                return None, None, None, None
            srt = rslt[0]
            srt = f * np.dot(dmat.T, srt) / scale
            bcm += np.outer(srt, srt)

        cov_robust_bc = np.dot(cov_naive, np.dot(bcm, cov_naive))

        cov_naive *= self.scaling_factor
        cov_robust *= self.scaling_factor
        cov_robust_bc *= self.scaling_factor

        return (cov_robust, cov_naive, cov_robust_bc, cmat)
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/genmod/generalized_estimating_equations.py/GEE._covmat/swapped operands of `-`