Repository: statsmodels/statsmodels

Issue:
ENH: confint_poisson only has two-sided confidence intervals
https://www.statsmodels.org/dev/generated/statsmodels.stats.rates.confint_poisson.html

confint_poisson does not have an one-sided confidence interval option.
As mentioned in #9249, add a new option 'alternative' to confint_poisson.

Maybe same for multinomial_proportions_confint, later.

Discussion:
mcmrc: I apologize for wasting computational resources. I checked the code with flake8 before committing, but no max-line-length errors were reported.
mcmrc: @josef-pkt 
Hi, could you review this PR?
josef-pkt: looks good, some corrections to docstring

then squash and merge

Thanks for the PR

Patch:
diff --git a/statsmodels/genmod/generalized_linear_model.py b/statsmodels/genmod/generalized_linear_model.py
index 6cd29fd76..4617ad4bc 100644
--- a/statsmodels/genmod/generalized_linear_model.py
+++ b/statsmodels/genmod/generalized_linear_model.py
@@ -41,7 +41,6 @@ from statsmodels.tools.decorators import (
     cached_data,
     cached_value,
 )
-from statsmodels.tools.data import _as_array_with_name
 from statsmodels.tools.docstring import Docstring
 from statsmodels.tools.sm_exceptions import (
     DomainWarning,
@@ -308,21 +307,15 @@ class GLM(base.LikelihoodModel):
                            f"{type(family).__name__} family."),
                           DomainWarning)
 
-        self._exposure_name = None
-        self._offset_name = None
-        self._freq_weights_name = None
-        self._var_weights_name = None
-
         if exposure is not None:
-            exposure_array, self._exposure_name = _as_array_with_name(exposure, "exposure")
-            exposure = np.log(exposure_array)
+            exposure = np.log(exposure)
         if offset is not None:  # this should probably be done upstream
-            offset, self._offset_name = _as_array_with_name(offset, "offset")
+            offset = np.asarray(offset)
 
         if freq_weights is not None:
-            freq_weights, self._freq_weights_name = _as_array_with_name(freq_weights, "freq_weights")
+            freq_weights = np.asarray(freq_weights)
         if var_weights is not None:
-            var_weights, self._var_weights_name = _as_array_with_name(var_weights, "var_weights")
+            var_weights = np.asarray(var_weights)
 
         self.freq_weights = freq_weights
         self.var_weights = var_weights
@@ -1565,39 +1558,6 @@ class GLM(base.LikelihoodModel):
         res._results.results_constrained = res_constr
         return res
 
-    @property
-    def offset_name(self):
-        """
-        Name of the offset variable if available. If offset is not a pd.Series,
-        defaults to 'offset'.
-        """
-        return self._offset_name
-
-    @property
-    def exposure_name(self):
-        """
-        Name of the exposure variable if available. If exposure is not a pd.Series,
-        defaults to 'exposure'.
-        """
-        return self._exposure_name
-
-    @property
-    def freq_weights_name(self):
-        """
-        Name of the freq weights variable if available. If freq_weights is not a
-        pd.Series, defaults to 'freq_weights'.
-        """
-        return self._freq_weights_name
-
-    @property
-    def var_weights_name(self):
-        """
-        Name of var weights variable if available. If var_weights is not a pd.Series,
-        defaults to 'var_weights'.
-
-        """
-        return self._var_weights_name
-
 
 get_prediction_doc = Docstring(pred.get_prediction_glm.__doc__)
 get_prediction_doc.remove_parameters("pred_kwds")
diff --git a/statsmodels/graphics/boxplots.py b/statsmodels/graphics/boxplots.py
index 526bfd911..2e142a2eb 100644
--- a/statsmodels/graphics/boxplots.py
+++ b/statsmodels/graphics/boxplots.py
@@ -144,13 +144,7 @@ def violinplot(data, ax=None, labels=None, positions=None, side='both',
         _single_violin(ax, pos, pos_data, width, side, plot_opts)
 
     if show_boxplot:
-        try:
-            ax.boxplot(
-                data, notch=1, positions=positions, orientation='vertical'
-            )
-        except TypeError:
-            # Remove after Matplotlib 3.10 is the minimum
-            ax.boxplot(data, notch=1, positions=positions, vert=1)
+        ax.boxplot(data, notch=1, positions=positions, vert=1)
 
     # Set ticks and tick labels of horizontal axis.
     _set_ticks_labels(ax, data, labels, positions, plot_opts)
diff --git a/statsmodels/robust/covariance.py b/statsmodels/robust/covariance.py
index f8e215082..b6897de3d 100644
--- a/statsmodels/robust/covariance.py
+++ b/statsmodels/robust/covariance.py
@@ -28,15 +28,9 @@ import numpy as np
 from scipy import stats, linalg
 from scipy.linalg.lapack import dtrtri
 from .scale import mad
-import statsmodels.robust.norms as rnorms
-import statsmodels.robust.scale as rscale
-import statsmodels.robust.tools as rtools
 from statsmodels.tools.testing import Holder
 
-from statsmodels.stats.covariance import corr_rank, corr_normal_scores
-
 mad0 = lambda x: mad(x, center=0)  # noqa: E731
-median = lambda x: np.median(x, axis=0)  # noqa: E731
 
 
 # from scikit-learn
@@ -317,7 +311,6 @@ def mahalanobis(data, cov=None, cov_inv=None, sqrt=False):
     """Mahalanobis distance squared
 
     Note: this is without taking the square root.
-    assumes data is already centered.
 
     Parameters
     ----------
@@ -412,15 +405,13 @@ def cov_gk(data, scale_func=mad):
 def cov_ogk(data, maxiter=2, scale_func=mad, cov_func=cov_gk,
             loc_func=lambda x: np.median(x, axis=0), reweight=0.9,
             rescale=True, rescale_raw=True, ddof=1):
-    """Orthogonalized Gnanadesikan and Kettenring covariance estimator.
+    """orthogonalized Gnanadesikan and Kettenring covariance estimator
 
     Based on Maronna and Zamar 2002
 
     Parameters
     ----------
-    data : array-like
-        Multivariate data set with observation in rows and variables in
-        columns.
+    data : array_like, 2-D
     maxiter : int
         Number of iteration steps. According to Maronna and Zamar the
         estimate doesn't improve much after the second iteration and the
@@ -599,10 +590,12 @@ def cov_tyler(data, start_cov=None, normalize=False, maxiter=100, eps=1e-13):
     ----------
     .. [1] Tyler, David E. “A Distribution-Free M-Estimator of Multivariate
        Scatter.” The Annals of Statistics 15, no. 1 (March 1, 1987): 234–51.
+
     .. [2] Soloveychik, I., and A. Wiesel. 2014. Tyler's Covariance Matrix
        Estimator in Elliptical Models With Convex Structure.
        IEEE Transactions on Signal Processing 62 (20): 5251-59.
        doi:10.1109/TSP.2014.2348951.
+
     .. [3] Ollila, Esa, Daniel P. Palomar, and Frederic Pascal.
        “Affine Equivariant Tyler’s M-Estimator Applied to Tail Parameter
        Learning of Elliptical Distributions.” arXiv, May 7, 2023.
@@ -1090,7 +1083,7 @@ def _cov_iter(data, weights_func, weights_args=None, cov_init=None,
     return res
 
 
-def _cov_starting(data, standardize=False, quantile=0.5, retransform=False):
+def _cov_starting(data, standardize=False, quantile=0.5):
     """compute some robust starting covariances
 
     The returned covariance matrices are intended as starting values
@@ -1134,7 +1127,7 @@ def _cov_starting(data, standardize=False, quantile=0.5, retransform=False):
 
     cov_all = []
     d = mahalanobis(xs, cov=None, cov_inv=np.eye(k_vars))
-    percentiles = [(k_vars+2) / nobs * 100 * 2, 25, 50, 85]
+    percentiles = [(k_vars+2) / nobs * 100 * 2, 25, 50]
     cutoffs = np.percentile(d, percentiles)
     for p, cutoff in zip(percentiles, cutoffs):
         xsp = xs[d < cutoff]
@@ -1156,929 +1149,15 @@ def _cov_starting(data, standardize=False, quantile=0.5, retransform=False):
         c03 = _cov_iter(xs, weights_quantile, weights_args=(quantile,),
                         rescale="med", cov_init=c02.cov, maxiter=100)
 
-        if not standardize or not retransform:
+        if standardize:
             cov_all.extend([c0, c01, c02, c03])
         else:
             # compensate for initial rescaling
-            # TODO: this does not return list of Holder anymore
             s = np.outer(std, std)
             cov_all.extend([r.cov * s for r in [c0, c01, c02, c03]])
 
     c2 = cov_ogk(xs)
     cov_all.append(c2)
 
-    c2raw = Holder(
-            cov=c2.cov_raw,
-            mean=c2.loc_raw * std + center,
-            method="ogk_raw",
-            )
-    cov_all.append(c2raw)
-
-    z_tanh = np.tanh(xs)
-    c_th = Holder(
-            cov=np.corrcoef(z_tanh.T),  # not consistently scaled for cov
-            mean=center,  # TODO: do we add inverted mean z_tanh ?
-            method="tanh",
-            )
-    cov_all.append(c_th)
-
-    x_spatial = xs / np.sqrt(np.sum(xs**2, axis=1))[:, None]
-    c_th = Holder(
-            cov=np.cov(x_spatial.T),
-            mean=center,
-            method="spatial",
-            )
-    cov_all.append(c_th)
-
-    c_th = Holder(
-            # not consistently scaled for cov
-            # cov=stats.spearmanr(xs)[0], # not correct shape if k=1 or 2
-            cov=corr_rank(xs),  # always returns matrix, np.corrcoef result
-            mean=center,
-            method="spearman",
-            )
-    cov_all.append(c_th)
-
-    c_ns = Holder(
-            cov=corr_normal_scores(xs),  # not consistently scaled for cov
-            mean=center,  # TODO: do we add inverted mean z_tanh ?
-            method="normal-scores",
-            )
-    cov_all.append(c_ns)
-
     # TODO: rescale back to original space using center and std
     return cov_all
-
-
-# ####### Det, CovDet and helper functions, might be moved to separate module
-
-
-class _Standardize():
-    """Robust standardization of random variable
-
-    """
-
-    def __init__(self, x, func_center=None, func_scale=None):
-        # naming mean or center
-        # maybe also allow str func_scale for robust.scale, e.g. for not
-        #    vectorized Qn
-        if func_center is None:
-            center = np.median(x, axis=0)
-        else:
-            center = func_center(x)
-        xdm = x - center
-        if func_scale is None:
-            scale = mad(xdm, center=0)
-        else:
-            # assumes vectorized
-            scale = func_scale(x)
-
-        self.x_stand = xdm / scale
-
-        self.center = center
-        self.scale = scale
-
-    def transform(self, x):
-        return (x - self.center) / self.scale
-
-    def untransform_mom(self, m, c):
-        mean = self.center + m
-        cov = c * np.outer(self.scale, self.scale)
-        return mean, cov
-
-
-def _orthogonalize_det(x, corr, loc_func, scale_func):
-    """Orthogonalize
-
-    This is a simplified version of the OGK method.
-    version from DetMCD works on zscored data
-    (does not return mean and cov of original data)
-    so we drop the compensation for scaling in zscoring
-
-    z is the data here, zscored with robust estimators,
-    e.g. median and Qn in DetMCD
-    """
-    evals, evecs = np.linalg.eigh(corr)  # noqa: F841
-    z = x.dot(evecs)
-    transf0 = evecs
-
-    scale_z = scale_func(z)  # scale of principal components
-    cov = (transf0 * scale_z**2).dot(transf0.T)
-    # extra step in DetMCD, sphering data with new cov to compute center
-    # I think this is equivalent to scaling z
-    # loc_z = loc_func(z / scale_z) * scale_z  # center of principal components
-    # loc = (transf0 * scale_z).dot(loc_z)
-    transf1 = (transf0 * scale_z).dot(transf0.T)
-    # transf1inv = (transf0 * scale_z**(-1)).dot(transf0.T)
-
-    # loc = loc_func(x @ transf1inv) @ transf1
-    loc = loc_func((z / scale_z).dot(transf0.T)) @ transf1
-
-    return loc, cov
-
-
-def _get_detcov_startidx(z, h, options_start=None, methods_cov="all"):
-    """Starting sets for deterministic robust covariance estimators.
-
-    These are intended as starting sets for DetMCD, DetS and DetMM.
-    """
-
-    if options_start is None:
-        options_start = {}
-
-    loc_func = options_start.get("loc_func", median)
-    scale_func = options_start.get("scale_func", mad)
-    z = (z - loc_func(z)) / scale_func(z)
-
-    if np.squeeze(z).ndim == 1:
-        # only one random variable
-        z = np.squeeze(z)
-        nobs = z.shape[0]
-        idx_sel = np.argpartition(np.abs(z), h)[:h]
-        idx_all = [(idx_sel, "abs-resid")]
-        # next uses symmetric equal-tail trimming
-        idx_sorted = np.argsort(z)
-        h_tail = (nobs - h) // 2
-        idx_all.append((idx_sorted[h_tail : h_tail + h], "trimmed-tail"))
-        return idx_all
-
-    # continue if more than 1 random variable
-    cov_all = _cov_starting(z, standardize=False, quantile=0.5)
-
-    # orthogonalization step
-    idx_all = []
-    for c in cov_all:
-        if not hasattr(c, "method"):
-            continue
-        method = c.method
-        mean, cov = _orthogonalize_det(z, c.cov, loc_func, scale_func)
-        d = mahalanobis(z, mean, cov)
-        idx_sel = np.argpartition(d, h)[:h]
-        idx_all.append((idx_sel, method))
-
-    return idx_all
-
-
-class CovM:
-    """M-estimator for multivariate Mean and Scatter.
-
-    Interface incomplete and experimental.
-
-    Parameters
-    ----------
-    data : array-like
-        Multivariate data set with observation in rows and variables in
-        columns.
-    norm_mean : norm instance
-        If None, then TukeyBiweight norm is used.
-        (Currently no other norms are supported for calling the initial
-        S-estimator)
-    norm_scatter : None or norm instance
-        If norm_scatter is None, then the norm_mean will be used.
-    breakdown_point : float in (0, 0.5]
-        Breakdown point for first stage S-estimator.
-    scale_bias : None or float
-        Must currently be provided if norm_mean is not None.
-    method : str
-        Currently only S-estimator has automatic selection of scale function.
-    """
-
-    def __init__(self, data, norm_mean=None, norm_scatter=None,
-                 scale_bias=None, method="S"):
-        # todo: method defines how norm_mean and norm_scatter are linked
-        #       currently I try for S-estimator
-
-        if method.lower() not in ["s"]:
-            msg = f"method {method} option not recognize or implemented"
-            raise ValueError(msg)
-
-        self.data = np.asarray(data)
-        self.k_vars = k_vars = self.data.shape[1]
-
-        # Todo: check interface for scale bias
-        self.scale_bias = scale_bias
-        if norm_mean is None:
-            norm_mean = rnorms.TukeyBiweight()
-            c = rtools.tuning_s_cov(norm_mean, k_vars, breakdown_point=0.5)
-            norm_mean._set_tuning_param(c, inplace=True)
-            self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
-
-        self.norm_mean = norm_mean
-        if norm_scatter is None:
-            self.norm_scatter = self.norm_mean
-        else:
-            self.norm_scatter = norm_scatter
-
-        self.weights_mean = self.norm_mean.weights
-        self.weights_scatter = self.weights_mean
-        # self.weights_scatter = lambda d: self.norm_mean.rho(d) / d**2
-        # this is for S-estimator, M-scale
-        self.rho = self.norm_scatter.rho
-
-    def _fit_mean_shape(self, mean, shape, scale):
-        """Estimate mean and shape in iteration step.
-
-        This does only one step.
-
-        Parameters
-        ----------
-        mean : ndarray
-            Starting value for mean
-        shape : ndarray
-            Starting value for shape matrix.
-        scale : float
-            Starting value for scale.
-
-        Returns
-        -------
-        Holder instance with updated estimates.
-        """
-        d = mahalanobis(self.data - mean, shape, sqrt=True) / scale
-        weights_mean = self.weights_mean(d)
-        weights_cov = self.weights_scatter(d)
-
-        res = cov_weighted(
-            self.data,
-            weights=weights_mean,
-            center=None,
-            weights_cov=weights_cov,
-            weights_cov_denom="det",
-            ddof=1,
-            )
-        return res
-
-    def _fit_scale(self, maha, start_scale=None, maxiter=100, rtol=1e-5,
-                   atol=1e-5):
-        """Estimate iterated M-scale.
-
-        Parameters
-        ----------
-        maha : ndarray
-        start_scale : None or float
-            Starting scale. If it is None, the mad of maha wi
-        maxiter : int
-            Maximum iterations to compute M-scale
-        rtol, atol : float
-            Relative and absolute convergence criteria for scale used with
-            allclose.
-
-        Returns
-        -------
-        float : scale estimate
-        """
-        if start_scale is None:
-            # TODO: this does not really make sense
-            # better scale to median of maha and chi or chi2
-            start_scale = mad(maha)
-
-        scale = rscale._scale_iter(
-            maha,
-            scale0=start_scale,
-            maxiter=maxiter,
-            rtol=rtol,
-            atol=atol,
-            meef_scale=self.rho,
-            scale_bias=self.scale_bias,
-            )
-        return scale
-
-    def fit(self, start_mean=None, start_shape=None, start_scale=None,
-            maxiter=100, update_scale=True):
-        """Estimate mean, shape and scale parameters with MM-estimator.
-
-        Parameters
-        ----------
-        start_mean : None or float
-            Starting value for mean, center.
-            If None, then median is used.
-        start_shape : None or 2-dim ndarray
-            Starting value of shape matrix, i.e. scatter matrix normalized
-            to det(scatter) = 1.
-            If None, then scaled covariance matrix of data is used.
-        start_scale : None or float.
-            Starting value of scale.
-        maxiter : int
-            Maximum number of iterations.
-        update_scale : bool
-            If update_scale is False, then
-
-        Returns
-        -------
-        results instance with mean, shape, scale, cov and other attributes.
-
-        Notes
-        -----
-        If start_scale is provided and update_scale is False, then this is
-        an M-estimator with a predetermined scale as used in the second
-        stage of an MM-estimator.
-
-        """
-
-        converged = False
-
-        if start_scale is not None:
-            scale_old = start_scale
-        else:
-            scale_old = 1
-            # will be reset if start_shape is also None.
-        if start_mean is not None:
-            mean_old = start_mean
-        else:
-            mean_old = np.median(self.data, axis=0)
-        if start_shape is not None:
-            shape_old = start_shape
-        else:
-            shape_old = np.cov(self.data.T)
-            scale = np.linalg.det(shape_old) ** (1 / self.k_vars)
-            shape_old /= scale
-            if start_scale is not None:
-                scale_old = scale
-
-        if update_scale is False:
-            scale = start_scale
-
-        for i in range(maxiter):
-            shape, mean = self._fit_mean_shape(mean_old, shape_old, scale_old)
-            d = mahalanobis(self.data - mean, shape, sqrt=True)
-            if update_scale:
-                scale = self._fit_scale(d, start_scale=scale_old, maxiter=10)
-
-            if (np.allclose(scale, scale_old, rtol=1e-5) and
-                    np.allclose(mean, mean_old, rtol=1e-5) and
-                    np.allclose(shape, shape_old, rtol=1e-5)
-                    ):  # noqa E124
-                converged = True
-                break
-            scale_old = scale
-            mean_old = mean
-            shape_old = shape
-
-        maha = mahalanobis(self.data - mean, shape / scale, sqrt=True)
-
-        res = Holder(
-            mean=mean,
-            shape=shape,
-            scale=scale,
-            cov=shape * scale**2,
-            converged=converged,
-            n_iter=i,
-            mahalanobis=maha,
-            )
-        return res
-
-
-class CovDetMCD:
-    """Minimum covariance determinant estimator with deterministic starts.
-
-    preliminary version
-
-    reproducability:
-
-    This uses deterministic starting sets and there is no randomness in the
-    estimator.
-    However, this will not be reprodusible across statsmodels versions
-    when the methods for starting sets or tuning parameters for the
-    optimization change.
-
-    Parameters
-    ----------
-    data : array-like
-        Multivariate data set with observation in rows and variables in
-        columns.
-
-    Notes
-    -----
-    The correction to the scale to take account of trimming in the reweighting
-    estimator is based on the chisquare tail probability.
-    This differs from CovMcd in R which uses the observed fraction of
-    observations above the metric trimming threshold.
-
-    References
-    ----------
-    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
-       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
-       Computational Statistics & Data Analysis 81 (January): 64–75.
-       https://doi.org/10.1016/j.csda.2014.07.013.
-    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
-       Deterministic Algorithm for Robust Location and Scatter.” Journal of
-       Computational and Graphical Statistics 21 (3): 618–37.
-       https://doi.org/10.1080/10618600.2012.672100.
-    """
-
-    def __init__(self, data):
-        # no options yet, methods were written as functions
-        self.data = np.asarray(data)
-
-    def _cstep(self, x, mean, cov, h, maxiter=2, tol=1e-8):
-        """C-step for mcd iteration
-
-        x is data, perc is percentile h / nobs, don't need perc when we
-        use np.argpartition
-        requires starting mean and cov
-        """
-
-        converged = False
-
-        for _ in range(maxiter):
-            d = mahalanobis(x - mean, cov)
-            idx_sel = np.argpartition(d, h)[:h]
-            x_sel = x[idx_sel]
-            mean = x_sel.mean(0)
-            cov_new = np.cov(x_sel.T, ddof=1)
-
-            if ((cov - cov_new)**2).mean() < tol:
-                cov = cov_new
-                converged = True
-                break
-
-            cov = cov_new
-
-        return mean, cov, converged
-
-    def _fit_one(self, x, idx, h, maxiter=2, mean=None, cov=None):
-        """Compute mcd for one starting set of observations.
-
-        Parameters
-        ----------
-        x : ndarray
-            Data.
-        idx : ndarray
-            Indices or mask of observation in starting set, used as ``x[idx]``
-        h : int
-            Number of observations in evaluation set for cov.
-        maxiter : int
-            Maximum number of c-steps.
-
-        Returns
-        -------
-        mean : ndarray
-            Estimated mean.
-        cov : ndarray
-            Estimated covariance.
-        det : float
-            Determinant of estimated covariance matrix.
-
-        Notes
-        -----
-        This does not do any preprocessing of the data and returns the
-        empirical mean and covariance of evaluation set of the data ``x``.
-        """
-        if idx is not None:
-            x_sel = x[idx]
-        else:
-            x_sel = x
-        if mean is None:
-            mean = x_sel.mean(0)
-        if cov is None:
-            cov = np.cov(x_sel.T, ddof=1)
-
-        # updated with c-step
-        mean, cov, conv = self._cstep(x, mean, cov, h, maxiter=maxiter)
-        det = np.linalg.det(cov)
-
-        return mean, cov, det, conv
-
-    def fit(self, h, *, h_start=None, mean_func=None, scale_func=None,
-            maxiter=100, options_start=None, reweight=True,
-            trim_frac=0.975, maxiter_step=100):
-        """
-        Compute minimum covariance determinant estimate of mean and covariance.
-
-        x : array-like
-            Data with observation in rows and variables in columns.
-        h : int
-            Number of observations in evaluation set for minimimizing
-            determinant.
-        h_start : int
-            Number of observations used in starting mean and covariance.
-        mean_func, scale_func : callable or None.
-            Mean and scale function for initial standardization.
-            Current defaults, if they are None, are median and mad, but
-            default scale_func will likely change.
-        options_start : None or dict
-            Options for the starting estimators.
-            currently not used
-            TODO: which options? e.g. for OGK
-        reweight : bool
-            If reweight is true, then a reweighted estimator is returned. The
-            reweighting is based on a chisquare trimming of Mahalanobis
-            distances. The raw results are in the ``results_raw`` attribute.
-        trim_frac : float in (0, 1)
-            Trim fraction used if reweight is true. Used to compute quantile
-            of chisquare distribution with tail probability 1 - trim_frac.
-        maxiter_step : int
-            Number of iteration in the c-step.
-            In the current implementation a small maxiter in the c-step does
-            not find the optimal solution.
-
-        Returns
-        -------
-        Holder instance with results
-        """
-
-        x = self.data
-        nobs, k_vars = x.shape
-
-        if h is None:
-            h = (nobs + k_vars + 1) // 2  # check with literature
-        if mean_func is None:
-            mean_func = lambda x: np.median(x, axis=0)  # noqa
-        if scale_func is None:
-            scale_func = mad
-        if options_start is None:
-            options_start = {}
-        if h_start is None:
-            nobs, k_vars = x.shape
-            h_start = max(nobs // 2 + 1, k_vars + 1)
-
-        m = mean_func(x)
-        s = scale_func(x)
-        z = (x - m) / s
-        # get initial mean, cov of standardized data, we only need ranking
-        # of obs
-        starts = _get_detcov_startidx(z, h_start, options_start)
-
-        fac_trunc = coef_normalize_cov_truncated(h / nobs, k_vars)
-
-        res = {}
-        for ii, ini in enumerate(starts):
-            idx_sel, method = ini
-            mean, cov, det, _ = self._fit_one(x, idx_sel, h,
-                                              maxiter=maxiter_step)
-            res[ii] = Holder(
-                mean=mean,
-                cov=cov * fac_trunc,
-                det_subset=det,
-                method=method,
-                )
-
-        det_all = np.array([i.det_subset for i in res.values()])
-        idx_best = np.argmin(det_all)
-        best = res[idx_best]
-        # mean = best.mean
-        # cov = best.cov
-
-        # need to c-step to convergence for best,
-        # is with best 2 in original DetMCD
-        if maxiter_step < maxiter:
-            mean, cov, det, conv = self._fit_one(x, None, h, maxiter=maxiter,
-                                                 mean=best.mean, cov=best.cov)
-            best = Holder(
-                mean=mean,
-                cov=cov * fac_trunc,
-                det_subset=det,
-                method=method,
-                converged=conv,
-                )
-
-        # include extra info in returned Holder instance
-        best.det_all = det_all
-        best.idx_best = idx_best
-
-        best.tmean = m
-        best.tscale = s
-
-        if reweight:
-            cov, mean = _reweight(x, best.mean, best.cov, trim_frac=trim_frac,
-                                  ddof=1)
-            fac_trunc = coef_normalize_cov_truncated(trim_frac, k_vars)
-            best_w = Holder(
-                mean=mean,
-                cov=cov * fac_trunc,
-                # det_subset=det,
-                method=method,
-                results_raw=best,
-                )
-
-            return best_w
-        else:
-            return best  # is Holder instance already
-
-
-class CovDetS:
-    """S-estimator for mean and covariance with deterministic starts.
-
-    Parameters
-    ----------
-    data : array-like
-        Multivariate data set with observation in rows and variables in
-        columns.
-    norm : norm instance
-        If None, then TukeyBiweight norm is used.
-        (Currently no other norms are supported for calling the initial
-        S-estimator)
-    breakdown_point : float in (0, 0.5]
-        Breakdown point for first stage S-estimator.
-
-    Notes
-    -----
-    Reproducability:
-
-    This uses deterministic starting sets and there is no randomness in the
-    estimator.
-    However, the estimates may not be reproducable across statsmodels versions
-    when the methods for starting sets or default tuning parameters for the
-    optimization change. With different starting sets, the estimate can
-    converge to a different local optimum.
-
-    References
-    ----------
-    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
-       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
-       Computational Statistics & Data Analysis 81 (January): 64–75.
-       https://doi.org/10.1016/j.csda.2014.07.013.
-    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
-       Deterministic Algorithm for Robust Location and Scatter.” Journal of
-       Computational and Graphical Statistics 21 (3): 618–37.
-       https://doi.org/10.1080/10618600.2012.672100.
-    """
-
-    def __init__(self, data, norm=None, breakdown_point=0.5):
-        # no options yet, methods were written as functions
-        self.data = np.asarray(data)
-        self.nobs, k_vars = self.data.shape
-        self.k_vars = k_vars
-
-        # self.scale_bias = scale_bias
-        if norm is None:
-            norm = rnorms.TukeyBiweight()
-            c = rtools.tuning_s_cov(norm, k_vars, breakdown_point=0.5)
-            norm._set_tuning_param(c, inplace=True)
-            self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
-        else:
-            raise NotImplementedError("only Biweight norm is supported")
-
-        self.norm = norm
-
-        self.mod = CovM(data, norm_mean=norm, norm_scatter=norm,
-                        scale_bias=self.scale_bias, method="S")
-
-    def _get_start_params(self, idx):
-        """Starting parameters from a subsample given by index
-
-        Parameters
-        ----------
-        idx : ndarray
-            Index used to select observations from the data. The index is used
-            for numpy arrays, so it can be either a boolean mask or integers.
-
-        Returns
-        -------
-        mean : ndarray
-            Mean of subsample
-        shape : ndarray
-            The shape matrix of the subsample which is the covariance
-            normalized so that determinant of shape is one.
-        scale : float
-            Scale of subsample, computed so that cov = shape * scale.
-        """
-        x_sel = self.data[idx]
-        k = x_sel.shape[1]
-
-        mean = x_sel.mean(0)
-        cov = np.cov(x_sel.T)
-
-        scale2 = np.linalg.det(cov) ** (1 / k)
-        shape = cov / scale2
-        scale = np.sqrt(scale2)
-        return mean, shape, scale
-
-    def _fit_one(self, mean=None, shape=None, scale=None, maxiter=100):
-        """Compute local M-estimator for one starting set of observations.
-
-        Parameters
-        ----------
-        x : ndarray
-            Data.
-        idx : ndarray
-            Indices or mask of observation in starting set, used as ``x[idx]``
-        h : int
-            Number of observations in evaluation set for cov.
-        maxiter : int
-            Maximum number of c-steps.
-
-        Returns
-        -------
-        mean : ndarray
-            Estimated mean.
-        cov : ndarray
-            Estimated covariance.
-        det : float
-            Determinant of estimated covariance matrix.
-
-        Notes
-        -----
-        This uses CovM to solve for the local optimum for given starting
-        values.
-        """
-
-        res = self.mod.fit(
-            start_mean=mean,
-            start_shape=shape,
-            start_scale=scale,
-            maxiter=maxiter,
-            update_scale=True
-            )
-
-        return res
-
-    def fit(self, *, h_start=None, mean_func=None, scale_func=None,
-            maxiter=100, options_start=None, maxiter_step=5):
-        """Compute S-estimator of mean and covariance.
-
-        Parameters
-        ----------
-        h_start : int
-            Number of observations used in starting mean and covariance.
-        mean_func, scale_func : callable or None.
-            Mean and scale function for initial standardization.
-            Current defaults, if they are None, are median and mad, but
-            default scale_func will likely change.
-        options_start : None or dict
-           Options for the starting estimators.
-           TODO: which options? e.g. for OGK
-
-        Returns
-        -------
-        Holder instance with results
-        """
-
-        x = self.data
-        nobs, k_vars = x.shape
-
-        if mean_func is None:
-            mean_func = lambda x: np.median(x, axis=0)  # noqa
-        if scale_func is None:
-            scale_func = mad
-        if options_start is None:
-            options_start = {}
-        if h_start is None:
-            nobs, k_vars = x.shape
-            h_start = max(nobs // 2 + 1, k_vars + 1)
-
-        m = mean_func(x)
-        s = scale_func(x)
-        z = (x - m) / s
-        # get initial mean, cov of standardized data, we only need ranking
-        # of obs
-        starts = _get_detcov_startidx(z, h_start, options_start)
-
-        res = {}
-        for ii, ini in enumerate(starts):
-            idx_sel, method = ini
-            mean0, shape0, scale0 = self._get_start_params(idx_sel)
-            res_i = self._fit_one(
-                mean=mean0,
-                shape=shape0,
-                scale=scale0,
-                maxiter=maxiter_step,
-                )
-
-            res_i.method = method
-            res[ii] = res_i
-
-        scale_all = np.array([i.scale for i in res.values()])
-        idx_best = np.argmin(scale_all)
-        best = res[idx_best]
-        # mean = best.mean
-        # cov = best.cov
-
-        # need to c-step to convergence for best,
-        # is with best 2 in original DetMCD
-        if maxiter_step < maxiter:
-            best = self._fit_one(
-                mean=best.mean,
-                shape=best.shape,
-                scale=best.scale,
-                maxiter=maxiter,
-                )
-
-        # include extra info in returned Holder instance
-        best.scale_all = scale_all
-        best.idx_best = idx_best
-
-        best.tmean = m
-        best.tscale = s
-
-        return best  # is Holder instance already
-
-
-class CovDetMM():
-    """MM estimator using DetS as first stage estimator.
-
-    Note: The tuning parameter for second stage M estimator is currently only
-    available for a small number of variables and only three values of
-    efficiency. For other cases, the user has to provide the norm instance
-    with desirec tuning parameter.
-
-    Parameters
-    ----------
-    data : array-like
-        Multivariate data set with observation in rows and variables in
-        columns.
-    norm : norm instance
-        If None, then TukeyBiweight norm is used.
-        (Currently no other norms are supported for calling the initial
-        S-estimator)
-        If ``norm`` is an instance of TukeyBiweight, then it will be used in
-        the second stage M-estimation. The ``efficiency`` argument is ignored
-        and the tuning parameter of the user provided instance is not changed.
-    breakdown_point : float in (0, 0.5]
-        Breakdown point for first stage S-estimator.
-    efficiency : float
-        Asymptotic efficiency of second stage M estimator.
-
-    Notes
-    -----
-    Current limitation is that only TukeyBiweight is supported.
-
-    The tuning parameter for second stage M estimator uses a table of values
-    for number of variables up to 15 and efficiency in
-    [0.75, 0.8, 0.85, 0.9, 0.95, 0.975, 0.99]. The tuning parameter for other
-    cases needs to be computed by numerical integration and rootfinding.
-    Alternatively, the user can provide a norm instance with desired tuning
-    parameter.
-
-    References
-    ----------
-    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
-       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
-       Computational Statistics & Data Analysis 81 (January): 64–75.
-       https://doi.org/10.1016/j.csda.2014.07.013.
-    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
-       Deterministic Algorithm for Robust Location and Scatter.” Journal of
-       Computational and Graphical Statistics 21 (3): 618–37.
-       https://doi.org/10.1080/10618600.2012.672100.
-    ..[3] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators and
-       M-Estimators of Multivariate Location and Covariance.” The Annals of
-       Statistics 17 (4): 1662–83.
-    ..[4] Salibián-Barrera, Matías, Stefan Van Aelst, and Gert Willems. 2006.
-       “Principal Components Analysis Based on Multivariate MM Estimators with
-       Fast and Robust Bootstrap.” Journal of the American Statistical
-       Association 101 (475): 1198–1211.
-    ..[5] Tatsuoka, Kay S., and David E. Tyler. 2000. “On the Uniqueness of
-       S-Functionals and M-Functionals under Nonelliptical Distributions.” The
-       Annals of Statistics 28 (4): 1219–43.
-    """
-
-    def __init__(self, data, norm=None, breakdown_point=0.5, efficiency=0.95):
-        self.data = np.asarray(data)
-        self.nobs, k_vars = self.data.shape
-        self.k_vars = k_vars
-        self.breakdown_point = breakdown_point
-
-        # self.scale_bias = scale_bias
-        if norm is None:
-            norm = rnorms.TukeyBiweight()
-            c = rtools.tukeybiweight_mvmean_eff(k_vars, efficiency)
-            norm._set_tuning_param(c, inplace=True)
-            # scale_bias is not used for second stage MM norm
-            # self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
-        elif not isinstance(norm, rnorms.TukeyBiweight):
-            raise NotImplementedError("only Biweight norm is supported")
-        # We allow tukeybiweight norm instance with user provided c
-
-        self.norm = norm
-        # model for second stage M-estimator
-        self.mod = CovM(data, norm_mean=norm, norm_scatter=norm,
-                        scale_bias=None, method="S")
-
-    def fit(self, maxiter=100):
-        """Estimate model parameters.
-
-        Parameters
-        ----------
-        maxiter : int
-            Maximum number of iterations in the second stage M-estimation.
-        fit args : dict
-            currently missing
-
-        Returns
-        -------
-        Instance of a results or holder class.
-
-        Notes
-        -----
-        This uses CovDetS for the first stage estimation and CovM with fixed
-        scale in the second stage MM-estimation.
-
-        TODO: fit options are missing.
-
-        """
-        # first stage estimate
-        mod_s = CovDetS(
-            self.data,
-            norm=None,
-            breakdown_point=self.breakdown_point
-            )
-        res_s = mod_s.fit()
-
-        res = self.mod.fit(
-            start_mean=res_s.mean,
-            start_shape=res_s.shape,
-            start_scale=res_s.scale,
-            maxiter=maxiter,
-            update_scale=False,
-            )
-
-        return res
diff --git a/statsmodels/robust/norms.py b/statsmodels/robust/norms.py
index 1e7edac3b..b3993a68e 100644
--- a/statsmodels/robust/norms.py
+++ b/statsmodels/robust/norms.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from . import tools as rtools
+# TODO: add plots to weighting functions for online docs.
 
 
 def _cabs(x):
@@ -41,9 +41,6 @@ class RobustNorm:
 
     continuous = 1
 
-    def __repr__(self):
-        return self.__class__.__name__
-
     def rho(self, z):
         """
         The robust criterion estimator function.
@@ -202,16 +199,12 @@ class HuberT(RobustNorm):
     def __init__(self, t=1.345):
         self.t = t
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        if inplace:
-            self.t = c
-            return self
-        else:
-            return self.__class__(t=c)
+        self.t = c
 
     def max_rho(self):
         return np.inf
@@ -329,18 +322,6 @@ class RamsayE(RobustNorm):
     def __init__(self, a=.3):
         self.a = a
 
-    def _set_tuning_param(self, c, inplace=False):
-        """Set and change the tuning parameter of the Norm.
-
-        Warning: this needs to wipe cached attributes that depend on the param.
-        """
-        # todo : change default to inplace=False, when tools are fixed
-        if inplace:
-            self.a = c
-            return self
-        else:
-            return self.__class__(a=c)
-
     def max_rho(self):
         return np.inf
 
@@ -438,16 +419,12 @@ class AndrewWave(RobustNorm):
     def __init__(self, a=1.339):
         self.a = a
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, a):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        if inplace:
-            self.a = c
-            return self
-        else:
-            return self.__class__(a=c)
+        self.a = a
 
     def max_rho(self):
         return 2 * self.a**2
@@ -576,16 +553,12 @@ class TrimmedMean(RobustNorm):
     def __init__(self, c=2.):
         self.c = c
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        if inplace:
-            self.c = c
-            return self
-        else:
-            return self.__class__(c=c)
+        self.c = c
 
     def max_rho(self):
         return self.rho(self.c)
@@ -701,20 +674,14 @@ class Hampel(RobustNorm):
         self.b = b
         self.c = c
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        a = c / 4
-        b = c / 2
-        if inplace:
-            self.c = c
-            self.a = a
-            self.b = b
-            return self
-        else:
-            return self.__class__(a=a, b=b, c=c)
+        self.c = c
+        self.a = c / 4
+        self.b = c / 2
 
     def max_rho(self):
         return self.rho(self.c)
@@ -895,49 +862,12 @@ class TukeyBiweight(RobustNorm):
     def __init__(self, c=4.685):
         self.c = c
 
-    def __repr__(self):
-        return f"{self.__class__.__name__}(c={self.c})"
-
-    @classmethod
-    def get_tuning(cls, bp=None, eff=None):
-        """Tuning parameter for given breakdown point or efficiency.
-
-        This currently only return values from a table.
-
-        Parameters
-        ----------
-        bp : float in [0.05, 0.5] or None
-            Required breakdown point
-            Either bp or eff has to be specified, but not both.
-        eff : float or None
-            Required asymptotic efficiency.
-            Either bp or eff has to be specified, but not both.
-
-        Returns
-        -------
-        float : tuning parameter.
-
-        """
-        if ((bp is None and eff is None) or
-                (bp is not None and eff is not None)):
-            raise ValueError("exactly one of bp and eff needs to be provided")
-
-        if bp is not None:
-            return rtools.tukeybiweight_bp[bp]
-        elif eff is not None:
-            return rtools.tukeybiweight_eff[eff]
-
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        # todo : change default to inplace=False, when tools are fixed
-        if inplace:
-            self.c = c
-            return self
-        else:
-            return self.__class__(c=c)
+        self.c = c
 
     def max_rho(self):
         return self.rho(self.c)
@@ -961,7 +891,7 @@ class TukeyBiweight(RobustNorm):
         Returns
         -------
         rho : ndarray
-            rho(z) = -(1 - (z/c)**2)**3 * c**2/6 + c**2/6   for \|z\| <= R
+            rho(z) = -(1 - (z/c)**2)**3 * c**2/6.   for \|z\| <= R
 
             rho(z) = 0                              for \|z\| > R
         """
@@ -1053,16 +983,12 @@ class TukeyQuartic(RobustNorm):
         self.c = c
         self.k = k
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        if inplace:
-            self.c = c
-            return self
-        else:
-            return self.__class__(c=c, k=self.k)
+        self.c = c
 
     def max_rho(self):
         return self.rho(self.c)
@@ -1195,16 +1121,12 @@ class StudentT(RobustNorm):
         self.c = c
         self.df = df
 
-    def _set_tuning_param(self, c, inplace=False):
+    def _set_tuning_param(self, c):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        if inplace:
-            self.c = c
-            return self
-        else:
-            return self.__class__(c=c, df=self.df)
+        self.c = c
 
     def max_rho(self):
         return np.inf
diff --git a/statsmodels/robust/resistant_linear_model.py b/statsmodels/robust/resistant_linear_model.py
deleted file mode 100644
index 843b616ab..000000000
--- a/statsmodels/robust/resistant_linear_model.py
+++ /dev/null
@@ -1,248 +0,0 @@
-"""
-Created on Apr. 19, 2024 12:17:03 p.m.
-
-Author: Josef Perktold
-License: BSD-3
-"""
-
-import numpy as np
-from statsmodels.base.model import Model
-from statsmodels.regression.linear_model import OLS
-from statsmodels.tools.testing import Holder
-
-from statsmodels.robust.robust_linear_model import RLM
-import statsmodels.robust.norms as rnorms
-import statsmodels.robust.scale as rscale
-from statsmodels.robust.covariance import _get_detcov_startidx
-
-
-class RLMDetS(Model):
-    """S-estimator for linear model with deterministic starts.
-
-
-    Notes
-    -----
-    This estimator combines the method of Fast-S regression (Saliban-Barrera
-    et al 2006) using starting sets similar to the deterministic estimation of
-    multivariate location and scatter DetS and DetMM of Hubert et al (2012).
-
-
-    References
-    ----------
-
-    .. [1] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
-       Deterministic Algorithm for Robust Location and Scatter.” Journal of
-       Computational and Graphical Statistics 21 (3): 618–37.
-       https://doi.org/10.1080/10618600.2012.672100.
-
-    .. [2] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck.
-       2015. “The DetS and DetMM Estimators for Multivariate Location and
-       Scatter.” Computational Statistics & Data Analysis 81 (January): 64–75.
-       https://doi.org/10.1016/j.csda.2014.07.013.
-
-    .. [3] Rousseeuw, Peter J., Stefan Van Aelst, Katrien Van Driessen, and
-        Jose Agulló. 2004. “Robust Multivariate Regression.”
-        Technometrics 46 (3): 293–305.
-
-    .. [4] Salibian-Barrera, Matías, and Víctor J. Yohai. 2006. “A Fast
-       Algorithm for S-Regression Estimates.” Journal of Computational and
-       Graphical Statistics 15 (2): 414–27.
-
-
-    """
-
-    def __init__(self, endog, exog, norm=None, breakdown_point=0.5,
-                 col_indices=None, include_endog=False):
-        super().__init__(endog, exog)
-
-        if norm is None:
-            norm = rnorms.TukeyBiweight()
-
-        tune = norm.get_tuning(bp=breakdown_point)
-        c = tune[0]
-        scale_bias = tune[2]
-        norm = norm._set_tuning_param(c, inplace=False)
-        self.mscale = rscale.MScale(norm, scale_bias)
-
-        self.norm = norm
-        self.breakdown_point = breakdown_point
-
-        # TODO: detect constant
-        if col_indices is None:
-            exog_start = self.exog[:, 1:]
-        else:
-            exog_start = self.exog[:, col_indices]
-
-        # data for robust mahalanobis distance of starting sets
-        if include_endog:
-            self.data_start = np.column_stack((endog, exog_start))
-        else:
-            self.data_start = exog_start
-
-    def _get_start_params(self, h):
-        # I think we should use iterator with yield
-
-        if self.data_start.shape[1] == 0 and self.exog.shape[1] == 1:
-            quantiles = np.quantile(self.endog, [0.25, 0.5, 0.75])
-            start_params_all = [np.atleast_1d([q]) for q in quantiles]
-            return start_params_all
-
-        starts = _get_detcov_startidx(
-            self.data_start, h, options_start=None, methods_cov="all")
-
-        start_params_all = [
-            OLS(self.endog[idx], self.exog[idx]).fit().params
-            for (idx, method) in starts
-            ]
-        return start_params_all
-
-    def _fit_one(self, start_params, maxiter=100):
-        mod = RLM(self.endog, self.exog, M=self.norm)
-        res = mod.fit(start_params=start_params,
-                      scale_est=self.mscale,
-                      maxiter=maxiter)
-        return res
-
-    def fit(self, h, maxiter=100, maxiter_step=5, start_params_extra=None):
-
-        start_params_all = self._get_start_params(h)
-        if start_params_extra:
-            start_params_all.extend(start_params_extra)
-        res = {}
-        for ii, sp in enumerate(start_params_all):
-            res_ii = self._fit_one(sp, maxiter=maxiter_step)
-            res[ii] = Holder(
-                scale=res_ii.scale,
-                params=res_ii.params,
-                method=ii,  # method  # TODO need start set method
-                )
-
-        scale_all = np.array([i.scale for i in res.values()])
-        scale_sorted = np.argsort(scale_all)
-        best_idx = scale_sorted[0]
-
-        # TODO: iterate until convergence if start fits are not converged
-        res_best = self._fit_one(res[best_idx].params, maxiter=maxiter)
-
-        # TODO: add extra start and convergence info
-        res_best._results.results_iter = res
-        # results instance of _fit_once has RLM as `model`
-        res_best.model_dets = self
-        return res_best
-
-
-class RLMDetSMM(RLMDetS):
-    """MM-estimator with S-estimator starting values.
-
-    Parameters
-    ----------
-    endog : array-like, 1-dim
-        Dependent, endogenous variable.
-    exog array-like, 1-dim
-        Inependent, exogenous regressor variables.
-    norm : robust norm
-        Redescending robust norm used for S- and MM-estimation.
-        Default is TukeyBiweight.
-    efficiency : float in (0, 1)
-        Asymptotic efficiency of the MM-estimator (used in second stage).
-    breakdown_point : float in (0, 0.5)
-        Breakdown point of the preliminary S-estimator.
-    col_indices : None or array-like of ints
-        Index of columns of exog to use in the mahalanobis distance computation
-        for the starting sets of the S-estimator.
-        Default is all exog except first column (constant). Todo: will change
-        when we autodetect the constant column
-    include_endog : bool
-        If true, then the endog variable is combined with the exog variables
-        to compute the the mahalanobis distances for the starting sets of the
-        S-estimator.
-
-    """
-    def __init__(self, endog, exog, norm=None, efficiency=0.95,
-                 breakdown_point=0.5, col_indices=None, include_endog=False):
-        super().__init__(
-            endog,
-            exog,
-            norm=norm,
-            breakdown_point=breakdown_point,
-            col_indices=col_indices,
-            include_endog=include_endog
-            )
-
-        self.efficiency = efficiency
-        if norm is None:
-            norm = rnorms.TukeyBiweight()
-
-        c = norm.get_tuning(eff=efficiency)[0]
-        norm = norm._set_tuning_param(c, inplace=False)
-        self.norm_mean = norm
-
-    def fit(self, h=None, scale_binding=False, start=None):
-        """Estimate the model
-
-        Parameters
-        ----------
-        h : int
-            The size of the initial sets for the S-estimator.
-            Default is ....  (todo)
-        scale_binding : bool
-            If true, then the scale is fixed in the second stage M-estimation,
-            i.e. this is the MM-estimator.
-            If false, then the high breakdown point M-scale is used also in the
-            second stage M-estimation if that estimated scale is smaller than
-            the scale of the preliminary, first stage S-estimato.
-        start : tuple or None
-            If None, then the starting parameters and scale for the second
-            stage M-estimation are taken from the fist stage S-estimator.
-            Alternatively, the starting parameters and starting scale can be
-            provided by the user as tuple (start_params, start_scale). In this
-            case the first stage S-estimation in skipped.
-        maxiter, other optimization parameters are still missing (todo)
-
-        Returns
-        -------
-        results instance
-
-        Notes
-        -----
-        If scale_binding is false, then the estimator is a standard
-        MM-estimator with fixed scale in the second stage M-estimation.
-        If scale_binding is true, then the estimator will try to find an
-        estimate with lower M-scale using the same scale-norm rho as in the
-        first stage S-estimator. If the estimated scale, is not smaller than
-        then the scale estimated in the first stage S-estimator, then the
-        fixed scale MM-estimator is returned.
-
-
-        """
-        norm_m = self.norm_mean
-        if start is None:
-            res_s = super().fit(h)
-            start_params = np.asarray(res_s.params)
-            start_scale = res_s.scale
-        else:
-            start_params, start_scale = start
-            res_s = None
-
-        mod_m = RLM(self.endog, self.exog, M=norm_m)
-        res_mm = mod_m.fit(
-            start_params=start_params,
-            start_scale=start_scale,
-            update_scale=False
-            )
-
-        if not scale_binding:
-            # we can compute this first and skip MM if scale decrease
-            mod_sm = RLM(self.endog, self.exog, M=norm_m)
-            res_sm = mod_sm.fit(
-                start_params=start_params,
-                scale_est=self.mscale
-                )
-
-        if not scale_binding and res_sm.scale < res_mm.scale:
-            res = res_sm
-        else:
-            res = res_mm
-
-        res._results.results_dets = res_s
-        return res
diff --git a/statsmodels/robust/robust_linear_model.py b/statsmodels/robust/robust_linear_model.py
index 2bd52909b..d3c105731 100644
--- a/statsmodels/robust/robust_linear_model.py
+++ b/statsmodels/robust/robust_linear_model.py
@@ -263,7 +263,6 @@ class RLM(base.LikelihoodModel):
             wls_results = lm.WLS(self.endog, self.exog).fit()
         else:
             start_params = np.asarray(start_params, dtype=np.double).squeeze()
-            start_params = np.atleast_1d(start_params)
             if (start_params.shape[0] != self.exog.shape[1] or
                     start_params.ndim != 1):
                 raise ValueError('start_params must by a 1-d array with {} '
@@ -277,8 +276,6 @@ class RLM(base.LikelihoodModel):
             self.scale = self._estimate_scale(wls_results.resid)
         elif start_scale:
             self.scale = start_scale
-            if not update_scale:
-                self.scale_est = scale_est = "fixed"
 
         history = dict(params=[np.inf], scale=[])
         if conv == 'coefs':
diff --git a/statsmodels/robust/scale.py b/statsmodels/robust/scale.py
index b9049a6f6..d866d52f1 100644
--- a/statsmodels/robust/scale.py
+++ b/statsmodels/robust/scale.py
@@ -426,9 +426,6 @@ class MScale:
         self.chi_func = chi_func
         self.scale_bias = scale_bias
 
-    def __repr__(self):
-        return repr(self.chi_func)
-
     def __call__(self, data, **kwds):
         return self.fit(data, **kwds)
 
diff --git a/statsmodels/robust/tools.py b/statsmodels/robust/tools.py
index d923b3776..a424d40f6 100644
--- a/statsmodels/robust/tools.py
+++ b/statsmodels/robust/tools.py
@@ -4,7 +4,7 @@ Created on Mar. 11, 2024 10:41:37 p.m.
 Author: Josef Perktold
 License: BSD-3
 """
-# flake8: noqa E731
+
 import numpy as np
 from scipy import stats, integrate, optimize
 
@@ -57,8 +57,8 @@ def _var_normal(norm):
 
 
     """
-    num = stats.norm.expect(lambda x: norm.psi(x)**2)  #noqa E731
-    denom = stats.norm.expect(lambda x: norm.psi_deriv(x))**2  #noqa E731
+    num = stats.norm.expect(lambda x: norm.psi(x)**2)
+    denom = stats.norm.expect(lambda x: norm.psi_deriv(x))**2
     return num / denom
 
 
@@ -100,7 +100,7 @@ def _var_normal_jump(norm):
 
 
     """
-    num = stats.norm.expect(lambda x: norm.psi(x)**2)  #noqa E731
+    num = stats.norm.expect(lambda x: norm.psi(x)**2)
 
     def func(x):
         # derivative normal pdf
@@ -142,6 +142,7 @@ def _get_tuning_param(norm, eff, kwd="c", kwargs=None, use_jump=False,
         efficiency.
 
     """
+    kwds = {} if kwargs is None else kwargs
     if bracket is None:
         bracket = [0.1, 10]
 
@@ -149,12 +150,12 @@ def _get_tuning_param(norm, eff, kwd="c", kwargs=None, use_jump=False,
         def func(c):
             # kwds.update({kwd: c})
             # return _var_normal(norm(**kwds)) - 1 / eff
-            norm._set_tuning_param(c, inplace=True)
+            norm._set_tuning_param(c)
             return _var_normal(norm) - 1 / eff
     else:
         def func(c):
-            norm._set_tuning_param(c, inplace=True)
-            return _var_normal_jump(norm) - 1 / eff
+            norm._set_tuning_param(c)
+            return _var_normal_jump(norm(**kwds) - 1 / eff)
 
     res = optimize.brentq(func, *bracket)
     return res
@@ -214,15 +215,14 @@ def tuning_s_estimator_mean(norm, breakdown=None):
 
     def func(c):
         norm_ = norm
-        norm_._set_tuning_param(c, inplace=True)
-        bp = (stats.norm.expect(lambda x : norm_.rho(x)) /  #noqa E731
-              norm_.max_rho())
+        norm_._set_tuning_param(c)
+        bp = stats.norm.expect(lambda x : norm_.rho(x)) / norm_.max_rho()
         return bp
 
     res = []
     for bp in bps:
-        c_bp = optimize.brentq(lambda c0: func(c0) - bp, 0.1, 10)  #noqa E731
-        norm._set_tuning_param(c_bp, inplace=True)  # inplace modification
+        c_bp = optimize.brentq(lambda c0: func(c0) - bp, 0.1, 10)
+        norm._set_tuning_param(c_bp)  # inplace modification
         eff = 1 / _var_normal(norm)
         b = stats.norm.expect(lambda x : norm.rho(x))
         res.append([bp, eff, c_bp, b])
@@ -242,321 +242,3 @@ def tuning_s_estimator_mean(norm, breakdown=None):
         )
 
     return res2
-
-
-def scale_bias_cov_biw(c, k_vars):
-    """Multivariate scale bias correction for TukeyBiweight norm.
-
-    This uses the chisquare distribution as reference distribution for the
-    squared Mahalanobis distance.
-    """
-    p = k_vars  # alias for formula
-    chip, chip2, chip4, chip6 = stats.chi2.cdf(c**2, [p, p + 2, p + 4, p + 6])
-    b = p / 2 * chip2 - p * (p + 2) / (2 * c**2) * chip4
-    b += p * (p + 2) * (p + 4) / (6 * c**4) * chip6 + c**2 / 6 * (1 - chip)
-    return b, b / (c**2 / 6)
-
-
-def scale_bias_cov(norm, k_vars):
-    """Multivariate scale bias correction.
-
-
-    Parameter
-    ---------
-    norm : norm instance
-        The rho function of the norm is used in the moment condition for
-        estimating scale.
-    k_vars : int
-        Number of random variables in the multivariate data.
-
-    Returns
-    -------
-    scale_bias: float
-    breakdown_point : float
-        Breakdown point computed as scale bias divided by max rho.
-    """
-
-    rho = lambda x: (norm.rho(np.sqrt(x)))  # noqa
-    scale_bias = stats.chi2.expect(rho, args=(k_vars,))
-    return scale_bias, scale_bias / norm.max_rho()
-
-
-def tuning_s_cov(norm, k_vars, breakdown_point=0.5, limits=()):
-    """Tuning parameter for multivariate S-estimator given breakdown point.
-    """
-    from .norms import TukeyBiweight  # avoid circular import
-
-    if not limits:
-        limits = (0.5, 30)
-
-    if isinstance(norm, TukeyBiweight):
-        def func(c):
-            return scale_bias_cov_biw(c, k_vars)[1] - breakdown_point
-    else:
-        norm = norm._set_tuning_param(2., inplace=False)  # create copy
-
-        def func(c):
-            norm._set_tuning_param(c, inplace=True)
-            return scale_bias_cov(norm, k_vars)[1] - breakdown_point
-
-    p_tune = optimize.brentq(func, limits[0], limits[1])
-    return p_tune
-
-
-def eff_mvmean(norm, k_vars):
-    """Efficiency for M-estimator of multivariate mean at normal distribution.
-
-    This also applies to estimators that are locally equivalent to an
-    M-estimator such as S- and MM-estimators.
-
-    Parameters
-    ----------
-    norm : instance of norm class
-    k_vars : int
-        Number of variables in multivariate random variable, i.e. dimension.
-
-    Returns
-    -------
-    eff : float
-        Asymptotic relative efficiency of mean at normal distribution.
-    alpha : float
-        Numerical integral. Efficiency is beta**2 / alpha
-    beta : float
-        Numerical integral.
-
-    Notes
-    -----
-    This implements equ. (5.3) p. 1671 in Lopuhaä 1989
-
-    References
-    ----------
-
-    .. [1] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators
-       and M-Estimators of Multivariate Location and Covariance.”
-       The Annals of Statistics 17 (4): 1662–83.
-
-    """
-    k = k_vars  # shortcut
-    f_alpha = lambda d: norm.psi(d)**2 / k
-    f_beta = lambda d: ((1 - 1 / k) * norm.weights(d) +
-                        1 / k * norm.psi_deriv(d))
-    alpha = stats.chi(k).expect(f_alpha)
-    beta = stats.chi(k).expect(f_beta)
-    return beta**2 / alpha, alpha, beta
-
-
-def eff_mvshape(norm, k_vars):
-    """Efficiency of M-estimator of multivariate shape at normal distribution.
-
-    This also applies to estimators that are locally equivalent to an
-    M-estimator such as S- and MM-estimators.
-
-    Parameters
-    ----------
-    norm : instance of norm class
-    k_vars : int
-        Number of variables in multivariate random variable, i.e. dimension.
-
-    Returns
-    -------
-    eff : float
-        Asymptotic relative efficiency of mean at normal distribution.
-    alpha : float
-        Numerical integral. Efficiency is beta**2 / alpha
-    beta : float
-        Numerical integral.
-
-    Notes
-    -----
-    This implements sigma_1 in equ. (5.5) p. 1671 in Lopuhaä 1989.
-    Efficiency of shape is approximately 1 / sigma1.
-
-    References
-    ----------
-
-    .. [1] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators
-       and M-Estimators of Multivariate Location and Covariance.”
-       The Annals of Statistics 17 (4): 1662–83.
-
-    """
-
-    k = k_vars  # shortcut
-    f_a = lambda d: k * (k + 2) * norm.psi(d)**2 * d**2
-    f_b = lambda d: norm.psi_deriv(d) * d**2 + (k + 1) * norm.psi(d) * d
-    a = stats.chi(k).expect(f_a)
-    b = stats.chi(k).expect(f_b)
-    return b**2 / a, a, b
-
-
-def tuning_m_cov_eff(norm, k_vars, efficiency=0.95, eff_mean=True, limits=()):
-    """Tuning parameter for multivariate M-estimator given efficiency.
-
-    This also applies to estimators that are locally equivalent to an
-    M-estimator such as S- and MM-estimators.
-
-    Parameters
-    ----------
-    norm : instance of norm class
-    k_vars : int
-        Number of variables in multivariate random variable, i.e. dimension.
-    efficiency : float < 1
-        Desired asymptotic relative efficiency of mean estimator.
-        Default is 0.95.
-    eff_mean : bool
-        If eff_mean is true (default), then tuning parameter is to achieve
-        efficiency of mean estimate.
-        If eff_mean is fale, then tuning parameter is to achieve efficiency
-        of shape estimate.
-    limits : tuple
-        Limits for rootfinding with scipy.optimize.brentq.
-        In some cases the interval limits for rootfinding can be too small
-        and not cover the root. Current default limits are (0.5, 30).
-
-    Returns
-    -------
-    float : Tuning parameter for the norm to achieve desired efficiency.
-        Asymptotic relative efficiency of mean at normal distribution.
-
-    Notes
-    -----
-    This uses numerical integration and rootfinding and will be
-    relatively slow.
-    """
-    if not limits:
-        limits = (0.5, 30)
-
-    # make copy of norm
-    norm = norm._set_tuning_param(1, inplace=False)
-
-    if eff_mean:
-        def func(c):
-            norm._set_tuning_param(c, inplace=True)
-            return eff_mvmean(norm, k_vars)[0] - efficiency
-    else:
-        def func(c):
-            norm._set_tuning_param(c, inplace=True)
-            return eff_mvshape(norm, k_vars)[0] - efficiency
-
-    p_tune = optimize.brentq(func, limits[0], limits[1])
-    return p_tune
-
-
-#  ##### tables
-
-tukeybiweight_bp = {
-    # breakdown point : (tuning parameter, efficiency, scale bias)
-    0.50: (1.547645, 0.286826, 0.199600),
-    0.45: (1.756059, 0.369761, 0.231281),
-    0.40: (1.987965, 0.461886, 0.263467),
-    0.35: (2.251831, 0.560447, 0.295793),
-    0.30: (2.560843, 0.661350, 0.327896),
-    0.25: (2.937015, 0.759040, 0.359419),
-    0.20: (3.420681, 0.846734, 0.390035),
-    0.15: (4.096255, 0.917435, 0.419483),
-    0.10: (5.182361, 0.966162, 0.447614),
-    0.05: (7.545252, 0.992424, 0.474424),
-    }
-
-tukeybiweight_eff = {
-    # efficiency : (tuning parameter, breakdown point)
-    0.65: (2.523102, 0.305646),
-    0.70: (2.697221, 0.280593),
-    0.75: (2.897166, 0.254790),
-    0.80: (3.136909, 0.227597),
-    0.85: (3.443690, 0.197957),
-    0.90: (3.882662, 0.163779),
-    0.95: (4.685065, 0.119414),
-    0.97: (5.596823, 0.087088),
-    0.98: (5.920719, 0.078604),
-    0.99: (7.041392, 0.056969),
-    }
-
-# relative efficiency for M and MM estimator of multivariate location
-# Table 2 from Kudraszow and Maronna JMA 2011
-# k in [1, 2, 3, 4, 5, 10], eff in [0.8, 0.9, 0.95
-# TODO: need to replace with more larger table and more digits.
-# (4, 0.95): 5.76 -. 5.81 to better match R rrcov and numerical integration
-tukeybiweight_mvmean_eff_km = {
-        (1, 0.8): 3.14, (2, 0.8): 3.51, (3, 0.8): 3.82, (4, 0.8): 4.1,
-        (5, 0.8): 4.34, (10, 0.8): 5.39,
-        (1, 0.9): 3.88, (2, 0.9): 4.28, (3, 0.9): 4.62, (4, 0.9): 4.91,
-        (5, 0.9): 5.18, (10, 0.9): 6.38,
-        (1, 0.95): 4.68, (2, 0.95): 5.12, (3, 0.95): 5.48, (4, 0.95): 5.81,
-        (5, 0.95): 6.1, (10, 0.95): 7.67,
-        }
-
-_table_biweight_mvmean_eff = np.array([
-    [2.89717, 3.13691, 3.44369, 3.88266, 4.68506, 5.59682, 7.04139],
-    [3.26396, 3.51006, 3.82643, 4.2821, 5.12299, 6.0869, 7.62344],
-    [3.5721, 3.82354, 4.14794, 4.61754, 5.49025, 6.49697, 8.10889],
-    [3.84155, 4.09757, 4.42889, 4.91044, 5.81032, 6.85346, 8.52956],
-    [4.08323, 4.34327, 4.68065, 5.17267, 6.09627, 7.17117, 8.90335],
-    [4.30385, 4.56746, 4.91023, 5.41157, 6.35622, 7.45933, 9.24141],
-    [4.50783, 4.77466, 5.12228, 5.63199, 6.59558, 7.72408, 9.5512],
-    [4.69828, 4.96802, 5.32006, 5.83738, 6.81817, 7.96977, 9.83797],
-    [4.87744, 5.14986, 5.506, 6.03022, 7.02679, 8.19958, 10.10558],
-    [5.04704, 5.32191, 5.68171, 6.21243, 7.22354, 8.41594, 10.35696],
-    [5.20839, 5.48554, 5.84879, 6.38547, 7.41008, 8.62071, 10.59439],
-    [5.36254, 5.64181, 6.00827, 6.55051, 7.58772, 8.81538, 10.81968],
-    [5.51034, 5.7916, 6.16106, 6.7085, 7.75752, 9.00118, 11.03428],
-    [5.65249, 5.9356, 6.3079, 6.86021, 7.92034, 9.17908, 11.23939],
-    ])
-
-_table_biweight_mvshape_eff = np.array([
-    [3.57210, 3.82354, 4.14794, 4.61754, 5.49025, 6.49697, 8.10889],
-    [3.84155, 4.09757, 4.42889, 4.91044, 5.81032, 6.85346, 8.52956],
-    [4.08323, 4.34327, 4.68065, 5.17267, 6.09627, 7.17117, 8.90335],
-    [4.30385, 4.56746, 4.91023, 5.41157, 6.35622, 7.45933, 9.24141],
-    [4.50783, 4.77466, 5.12228, 5.63199, 6.59558, 7.72408, 9.55120],
-    [4.69828, 4.96802, 5.32006, 5.83738, 6.81817, 7.96977, 9.83797],
-    [4.87744, 5.14986, 5.50600, 6.03022, 7.02679, 8.19958, 10.10558],
-    [5.04704, 5.32191, 5.68171, 6.21243, 7.22354, 8.41594, 10.35696],
-    [5.20839, 5.48554, 5.84879, 6.38547, 7.41008, 8.62071, 10.59439],
-    [5.36254, 5.64181, 6.00827, 6.55051, 7.58772, 8.81538, 10.81968],
-    [5.51034, 5.79160, 6.16106, 6.70849, 7.75752, 9.00118, 11.03428],
-    [5.65249, 5.93560, 6.30790, 6.86021, 7.92034, 9.17908, 11.23939],
-    [5.78957, 6.07443, 6.44938, 7.00630, 8.07692, 9.34991, 11.43603],
-    [5.92206, 6.20858, 6.58604, 7.14731, 8.22785, 9.51437, 11.62502],
-    ])
-
-
-def _convert_to_dict_mvmean_effs(eff_mean=True):
-    effs_mvmean = [0.75, 0.8, 0.85, 0.9, 0.95, 0.975, 0.99]
-    ks = list(range(1, 15))
-    if eff_mean:
-        table = _table_biweight_mvmean_eff
-    else:
-        table = _table_biweight_mvshape_eff
-    tp = {}
-    for i, k in enumerate(ks):
-        for j, eff in enumerate(effs_mvmean):
-            tp[(k, eff)] = table[i, j]
-
-    return tp
-
-
-tukeybiweight_mvmean_eff_d = _convert_to_dict_mvmean_effs(eff_mean=True)
-tukeybiweight_mvshape_eff_d = _convert_to_dict_mvmean_effs(eff_mean=False)
-
-
-def tukeybiweight_mvmean_eff(k, eff, eff_mean=True):
-    """tuning parameter for biweight norm to achieve efficiency for mv-mean.
-
-    Uses values from precomputed table if available, otherwise computes it
-    numerically and adds it to the module global dict.
-
-    """
-    if eff_mean:
-        table_dict = tukeybiweight_mvmean_eff_d
-    else:
-        table_dict = tukeybiweight_mvshape_eff_d
-
-    try:
-        tp = table_dict[(k, eff)]
-    except KeyError:
-        # compute and cache
-        from .norms import TukeyBiweight  # avoid circular import
-        norm = TukeyBiweight(c=1)
-        tp = tuning_m_cov_eff(norm, k, efficiency=eff, eff_mean=eff_mean)
-        table_dict[(k, eff)] = tp
-    return tp
diff --git a/statsmodels/stats/covariance.py b/statsmodels/stats/covariance.py
index 3346912c3..341d952c8 100644
--- a/statsmodels/stats/covariance.py
+++ b/statsmodels/stats/covariance.py
@@ -149,18 +149,6 @@ def transform_corr_normal(corr, method, return_var=False, possdef=True):
         return corr_n
 
 
-def corr_rank(data):
-    """Spearman rank correlation
-
-    simplified version of scipy.stats.spearmanr
-    """
-    x = np.asarray(data)
-    axisout = 0
-    ar = np.apply_along_axis(stats.rankdata, axisout, x)
-    corr = np.corrcoef(ar, rowvar=False)
-    return corr
-
-
 def corr_normal_scores(data):
     """Gaussian rank (normal scores) correlation
 
diff --git a/statsmodels/stats/rates.py b/statsmodels/stats/rates.py
index dd6654c7f..8b226c9be 100644
--- a/statsmodels/stats/rates.py
+++ b/statsmodels/stats/rates.py
@@ -186,15 +186,15 @@ def test_poisson(count, nobs, value, method=None, alternative="two-sided",
     return res
 
 
-def confint_poisson(count, exposure, method=None, alpha=0.05):
+def confint_poisson(count, exposure, method=None, alpha=0.05,
+                    alternative="two-sided"):
     """Confidence interval for a Poisson mean or rate
 
     The function is vectorized for all methods except "midp-c", which uses
     an iterative method to invert the hypothesis test function.
 
     All current methods are central, that is the probability of each tail is
-    smaller or equal to alpha / 2. The one-sided interval limits can be
-    obtained by doubling alpha.
+    smaller or equal to alpha / 2.
 
     Parameters
     ----------
@@ -209,10 +209,16 @@ def confint_poisson(count, exposure, method=None, alpha=0.05):
     alpha : float in (0, 1)
         Significance level, nominal coverage of the confidence interval is
         1 - alpha.
+    alternative : {"two-sider", "larger", "smaller")
+        default: "two-sided"
+        Specifies whether to calculate a two-sided or one-sided confidence
+        interval.
 
     Returns
     -------
     tuple (low, upp) : confidence limits.
+        When alternative is not "two-sided", lower or upper bound is set to
+        0 or inf respectively.
 
     Notes
     -----
@@ -270,7 +276,12 @@ def confint_poisson(count, exposure, method=None, alpha=0.05):
     """
     n = exposure  # short hand
     rate = count / exposure
-    alpha = alpha / 2  # two-sided
+
+    if alternative == 'two-sided':
+        alpha = alpha / 2
+    elif alternative not in ['larger', 'smaller']:
+        raise NotImplementedError(
+            f"alternative {alternative} is not available")
 
     if method is None:
         msg = "method needs to be specified, currently no default method"
@@ -365,6 +376,11 @@ def confint_poisson(count, exposure, method=None, alpha=0.05):
     else:
         raise ValueError("unknown method %s" % method)
 
+    if alternative == "larger":
+        ci = (0, ci[1])
+    elif alternative == "smaller":
+        ci = (ci[0], np.inf)
+
     ci = (np.maximum(ci[0], 0), ci[1])
     return ci
 
diff --git a/statsmodels/stats/weightstats.py b/statsmodels/stats/weightstats.py
index 0e098e574..57f136c9e 100644
--- a/statsmodels/stats/weightstats.py
+++ b/statsmodels/stats/weightstats.py
@@ -498,7 +498,7 @@ class DescrStatsW:
         tstat : float
             test statistic
         pvalue : float
-            pvalue of the z-test
+            pvalue of the t-test
 
         Notes
         -----
@@ -1524,7 +1524,7 @@ def ztest(
     tstat : float
         test statistic
     pvalue : float
-        pvalue of the z-test
+        pvalue of the t-test
 
     Notes
     -----
diff --git a/statsmodels/tools/data.py b/statsmodels/tools/data.py
index 9e507030c..115bc1699 100644
--- a/statsmodels/tools/data.py
+++ b/statsmodels/tools/data.py
@@ -19,8 +19,7 @@ def _check_period_index(x, freq="M"):
     if not inferred_freq.startswith(freq):
         raise ValueError("Expected frequency {}. Got {}".format(freq,
                                                                 inferred_freq))
-def is_series(obj):
-    return isinstance(obj, pd.Series)
+
 
 def is_data_frame(obj):
     return isinstance(obj, pd.DataFrame)
@@ -122,24 +121,3 @@ def _is_recarray(data):
         return isinstance(data, np.core.recarray)
     else:
         return isinstance(data, np.rec.recarray)
-
-def _as_array_with_name(obj, default_name):
-    """
-    Call np.asarray() on obj and attempt to get the name if its a Series.
-
-    Parameters
-    ----------
-    obj: pd.Series
-        Series to convert to an array
-    default_name: str
-        The default name to return in case the object isn't a pd.Series or has
-        no name attribute.
-
-    Returns
-    -------
-    array_and_name: tuple[np.ndarray, str]
-        The data casted to np.ndarra and the series name or None
-    """
-    if is_series(obj):
-        return (np.asarray(obj), obj.name)
-    return (np.asarray(obj), default_name)


Test patch:
diff --git a/statsmodels/genmod/tests/test_glm.py b/statsmodels/genmod/tests/test_glm.py
index b3b6876a4..fd169d36f 100644
--- a/statsmodels/genmod/tests/test_glm.py
+++ b/statsmodels/genmod/tests/test_glm.py
@@ -2661,62 +2661,3 @@ def test_tweedie_score():
             nhess = approx_hess_cs(pa, lambda x: model.loglike(x, scale=1))
             ahess = model.hessian(pa, scale=1)
             assert_allclose(nhess, ahess, atol=5e-8, rtol=5e-8)
-
-def test_names():
-    """Test the name properties if using a pandas series.
-
-    They should not be the defaults if the series has a name.
-
-    Don't care about the data here, only testing the name properties.
-    """
-    y = pd.Series([0, 1], name="endog_not_default")
-    x = pd.DataFrame({"a": [1, 1], "b": [1, 0]})
-    exposure = pd.Series([0, 0], name="exposure_not_default")
-    freq_weights = pd.Series([0, 0], name="freq_weights_not_default")
-    offset = pd.Series([0, 0], name="offset_not_default")
-    var_weights = pd.Series([0, 0], name="var_weights_not_default")
-
-    model = GLM(
-        endog=y,
-        exog=x,
-        exposure=exposure,
-        freq_weights=freq_weights,
-        offset=offset,
-        var_weights=var_weights,
-        family=sm.families.Tweedie(),
-    )
-    assert model.offset_name == "offset_not_default"
-    assert model.exposure_name == "exposure_not_default"
-    assert model.freq_weights_name == "freq_weights_not_default"
-    assert model.var_weights_name == "var_weights_not_default"
-    assert model.endog_names == "endog_not_default"
-    assert model.exog_names == ["a", "b"]
-
-
-def test_names_default():
-    """Test the name properties if using a numpy arrays.
-
-    Don't care about the data here, only testing the name properties.
-    """
-    y = np.array([0, 1])
-    x = np.array([[1, 1,], [1, 0]])
-    exposure = np.array([0, 0])
-    freq_weights = np.array([0, 0])
-    offset = np.array([0, 0])
-    var_weights = np.array([0, 0])
-
-    model = GLM(
-        endog=y,
-        exog=x,
-        exposure=exposure,
-        freq_weights=freq_weights,
-        offset=offset,
-        var_weights=var_weights,
-        family=sm.families.Tweedie(),
-    )
-    assert model.offset_name == "offset"
-    assert model.exposure_name == "exposure"
-    assert model.freq_weights_name == "freq_weights"
-    assert model.var_weights_name == "var_weights"
-    assert model.endog_names == "y"
-    assert model.exog_names == ["const", "x1"]
diff --git a/statsmodels/robust/tests/test_covariance.py b/statsmodels/robust/tests/test_covariance.py
index 9aa60c5a0..ac1eb3c64 100644
--- a/statsmodels/robust/tests/test_covariance.py
+++ b/statsmodels/robust/tests/test_covariance.py
@@ -6,7 +6,6 @@ from numpy.testing import assert_allclose, assert_equal
 import pandas as pd
 
 from statsmodels import robust
-import statsmodels.robust.norms as robnorms
 import statsmodels.robust.covariance as robcov
 import statsmodels.robust.scale as robscale
 
@@ -178,121 +177,6 @@ def test_tyler():
     assert res1.n_iter == 55
 
 
-def test_cov_ms():
-    # use CovM as local CovS
-
-    # > CovSest(x)  using package rrcov
-    # same result with > CovSest(x, method="sdet")
-    # but scale difers with method="biweight"
-    mean_r = np.array([1.53420879676, 1.82865741024, 1.65565146981])
-    cov_r = np.array([
-        [1.8090846049573, 0.0479283121828, 0.2446369025717],
-        [0.0479283121828, 1.8189886310494, 0.2513025527579],
-        [0.2446369025717, 0.2513025527579, 1.7287983150484],
-        ])
-
-    scale2_r = np.linalg.det(cov_r) ** (1/3)
-    shape_r = cov_r / scale2_r
-    scale_r = np.sqrt(scale2_r)
-
-    exog_df = dta_hbk[["X1", "X2", "X3"]]
-    mod = robcov.CovM(exog_df)
-    # with default start, default start could get wrong local optimum
-    res = mod.fit()
-    assert_allclose(res.mean, mean_r, rtol=1e-5)
-    assert_allclose(res.shape, shape_r, rtol=1e-5)
-    assert_allclose(res.cov, cov_r, rtol=1e-5)
-    assert_allclose(res.scale, scale_r, rtol=1e-5)
-
-    # with results start
-    res = mod.fit(start_mean=mean_r, start_shape=shape_r, start_scale=scale_r)
-    assert_allclose(res.mean, mean_r, rtol=1e-5)
-    assert_allclose(res.shape, shape_r, rtol=1e-5)
-    assert_allclose(res.cov, cov_r, rtol=1e-5)
-    assert_allclose(res.scale, scale_r, rtol=1e-5)
-
-    mod_s = robcov.CovDetS(exog_df)
-    res = mod_s.fit()
-    assert_allclose(res.mean, mean_r, rtol=1e-5)
-    assert_allclose(res.shape, shape_r, rtol=1e-5)
-    assert_allclose(res.cov, cov_r, rtol=1e-5)
-    assert_allclose(res.scale, scale_r, rtol=1e-5)
-
-
-def test_covdetmcd():
-
-    # results from rrcov
-    # > cdet = CovMcd(x = hbk, raw.only = TRUE, nsamp = "deterministic",
-    #                 use.correction=FALSE)
-    cov_dmcd_r = np.array("""
-    2.2059619213639   0.0223939863695   0.7898958050933   0.4060613360808
-    0.0223939863695   1.1384166802155   0.4315534571891  -0.2344041030201
-    0.7898958050933   0.4315534571891   1.8930117467493  -0.3292893001459
-    0.4060613360808  -0.2344041030201  -0.3292893001459   0.6179686100845
-    """.split(), float).reshape(4, 4)
-
-    mean_dmcd_r = np.array([1.7725, 2.2050, 1.5375, -0.0575])
-
-    mod = robcov.CovDetMCD(dta_hbk)
-    res = mod.fit(40, maxiter_step=100, reweight=False)
-    assert_allclose(res.mean, mean_dmcd_r, rtol=1e-5)
-    assert_allclose(res.cov, cov_dmcd_r, rtol=1e-5)
-
-    # with reweighting
-    # covMcd(x = hbk, nsamp = "deterministic", use.correction = FALSE)
-    # iBest: 5; C-step iterations: 7, 7, 7, 4, 6, 6
-    # Log(Det.):  -2.42931967153
-
-    mean_dmcdw_r = np.array([1.5338983050847, 1.8322033898305, 1.6745762711864,
-                            -0.0728813559322])
-    cov_dmcdw_r = np.array("""
-    1.5677744869295   0.09285770205078   0.252076010128   0.13873444408300
-    0.0928577020508   1.56769177397171   0.224929617385  -0.00516128856542
-    0.2520760101278   0.22492961738467   1.483829106079  -0.20275013775619
-    0.1387344440830  -0.00516128856542  -0.202750137756   0.43326701543885
-    """.split(), float).reshape(4, 4)
-
-    mod = robcov.CovDetMCD(dta_hbk)
-    res = mod.fit(40, maxiter_step=100)  # default is reweight=True
-    assert_allclose(res.mean, mean_dmcdw_r, rtol=1e-5)
-    # R uses different trimming correction
-    # compare only shape (using trace for simplicity)
-    shape = res.cov / np.trace(res.cov)
-    shape_r = cov_dmcdw_r / np.trace(cov_dmcdw_r)
-    assert_allclose(shape, shape_r, rtol=1e-5)
-
-
-def test_covdetmm():
-
-    # results from rrcov
-    # CovMMest(x = hbk, eff.shape=FALSE,
-    #          control=CovControlMMest(sest=CovControlSest(method="sdet")))
-    cov_dmm_r = np.array("""
-        1.72174266670826 0.06925842715939 0.20781848922667 0.10749343153015
-        0.06925842715939 1.74566218886362 0.22161135221404 -0.00517023660647
-        0.20781848922667 0.22161135221404 1.63937749762534 -0.17217102475913
-        0.10749343153015 -0.00517023660647 -0.17217102475913 0.48174480967136
-        """.split(), float).reshape(4, 4)
-
-    mean_dmm_r = np.array([1.5388643420460, 1.8027582110408, 1.6811517253521,
-                           -0.0755069488908])
-
-    # using same c as rrcov
-    c = 5.81031555752526
-    mod = robcov.CovDetMM(dta_hbk, norm=robnorms.TukeyBiweight(c=c))
-    res = mod.fit()
-
-    assert_allclose(res.mean, mean_dmm_r, rtol=1e-5)
-    assert_allclose(res.cov, cov_dmm_r, rtol=1e-5, atol=1e-5)
-
-    # using c from table,
-    mod = robcov.CovDetMM(dta_hbk)
-    res = mod.fit()
-
-    assert_allclose(res.mean, mean_dmm_r, rtol=1e-3)
-    assert_allclose(res.cov, cov_dmm_r, rtol=1e-3, atol=1e-3)
-
-
 def test_robcov_SMOKE():
     # currently only smoke test or very loose comparisons to dgp
     nobs, k_vars = 100, 3
diff --git a/statsmodels/robust/tests/test_tools.py b/statsmodels/robust/tests/test_tools.py
index 9b10b23a8..a3adbb708 100644
--- a/statsmodels/robust/tests/test_tools.py
+++ b/statsmodels/robust/tests/test_tools.py
@@ -49,26 +49,24 @@ results_menenez = [
 @pytest.mark.parametrize("case", results_menenez)
 def test_eff(case):
     norm, res2 = case
-    use_jump = False
 
     if norm.continuous == 2:
         var_func = _var_normal
     else:
         var_func = _var_normal_jump
-        use_jump = True
 
     res_eff = []
     for c in res2:
-        norm._set_tuning_param(c, inplace=True)
+        norm._set_tuning_param(c)
         res_eff.append(1 / var_func(norm))
 
     assert_allclose(res_eff, effs, atol=0.0005)
 
     for c in res2:
         # bp = stats.norm.expect(lambda x : norm.rho(x)) / norm.rho(norm.c)
-        norm._set_tuning_param(c, inplace=True)
-        eff = 1 / var_func(norm)
-        tune = _get_tuning_param(norm, eff, use_jump=use_jump)
+        norm._set_tuning_param(c)
+        eff = 1 / _var_normal(norm)
+        tune = _get_tuning_param(norm, eff)
         assert_allclose(tune, c, rtol=1e-6, atol=5e-4)
 
 
diff --git a/statsmodels/stats/tests/test_rates_poisson.py b/statsmodels/stats/tests/test_rates_poisson.py
index 28baefb97..ff57cef6b 100644
--- a/statsmodels/stats/tests/test_rates_poisson.py
+++ b/statsmodels/stats/tests/test_rates_poisson.py
@@ -1146,3 +1146,24 @@ def test_power_negbin():
 
     assert_allclose(pow_p, pow1, atol=5e-2)
     assert_allclose(pow_p, pow_, rtol=1e-13)
+
+
+@pytest.mark.parametrize('count',    [0, 1, 10, 100])
+@pytest.mark.parametrize('exposure', [1, 10, 100, 1000])
+@pytest.mark.parametrize('alpha',    [0.01, 0.05, 0.1])
+@pytest.mark.parametrize('method',
+                         method_names_poisson_1samp['confint'])
+@pytest.mark.parametrize('alternative', ['larger', 'smaller'])
+def test_confint_poisson_alternative(count, exposure, method, alpha,
+                                     alternative):
+    # regression test
+    two_sided_ci = confint_poisson(count, exposure, method=method,
+                                   alpha=alpha * 2, alternative='two-sided')
+    one_sided_ci = confint_poisson(count, exposure, method=method,
+                                   alpha=alpha, alternative=alternative)
+    if alternative == 'larger':
+        two_sided_ci = (0, two_sided_ci[1])
+        assert_allclose(one_sided_ci, two_sided_ci, rtol=1e-12)
+    elif alternative == 'smaller':
+        two_sided_ci = (two_sided_ci[0], np.inf)
+        assert_allclose(one_sided_ci, two_sided_ci, rtol=1e-12)
diff --git a/statsmodels/tools/tests/test_data.py b/statsmodels/tools/tests/test_data.py
index 21a4f7f8e..178004aa1 100644
--- a/statsmodels/tools/tests/test_data.py
+++ b/statsmodels/tools/tests/test_data.py
@@ -33,16 +33,3 @@ def test_patsy_577():
     np.testing.assert_(data._is_using_patsy(endog, None))
     exog = dmatrix("var2 - 1", df)
     np.testing.assert_(data._is_using_patsy(endog, exog))
-
-
-def test_as_array_with_name_series():
-    s = pandas.Series([1], name="hello")
-    arr, name = data._as_array_with_name(s, "not_used")
-    np.testing.assert_array_equal(np.array([1]), arr)
-    assert name == "hello"
-
-
-def test_as_array_with_name_array():
-    arr, name = data._as_array_with_name(np.array([1]), "default")
-    np.testing.assert_array_equal(np.array([1]), arr)
-    assert name == "default"